我試圖模擬一個粒子在經歷電排斥（或吸引力）的同時向另一個粒子飛行，稱為盧瑟福散射。我已經成功地使用 for 循環和 python 列表模擬了（一些）粒子。但是，現在我想改用 numpy 數組。該模型將使用以下步驟：對于所有粒子：計算所有其他粒子的徑向距離計算與所有其他粒子的角度計算 x 方向和 y 方向的凈力使用 netto xForce 和 yForce 為每個粒子創建矩陣通過 a = F/mass 創建加速度（也是 x 和 y 分量）矩陣更新速度矩陣更新位置矩陣我的問題是我不知道如何使用 numpy 數組來計算力分量。 下面是我無法運行的代碼。import numpy as np# I used this function to calculate the force while using for-loops.def force(x1, y1, x2, x2): angle = math.atan((y2 - y1)/(x2 - x1)) dr = ((x1-x2)**2 + (y1-y2)**2)**0.5 force = charge2 * charge2 / dr**2 xforce = math.cos(angle) * force yforce = math.sin(angle) * force # The direction of force depends on relative location if x1 > x2 and y1<y2: xforce = xforce yforce = yforce elif x1< x2 and y1< y2: xforce = -1 * xforce yforce = -1 * yforce elif x1 > x2 and y1 > y2: xforce = xforce yforce = yforce else: xforce = -1 * xforce yforce = -1* yforce return xforce, yforcedef update(array): # this for loop defeats the entire use of numpy arrays for particle in range(len(array[0])): # find distance of all particles pov from 1 particle # find all x-forces and y-forces on that particle xforce = # sum of all x-forces from all particles yforce = # sum of all y-forces from all particles force_arr[0, particle] = xforce force_arr[1, particle] = yforce return force# begin parameterst = 0N = 3masses = np.ones(N)charges = np.ones(N)loc_arr = np.random.rand(2, N)speed_arr = np.random.rand(2, N)acc_arr = np.random.rand(2, N)force = np.random.rand(2, N)while t < 0.5: force_arr = update(loc_arry) acc_arr = force_arr / masses speed_arr += acc_array loc_arr += speed_arr t += dt # plot animation
3 回答
汪汪一只貓
TA貢獻1898條經驗 獲得超8個贊
用數組對這個問題建模的一種方法可能是：
將點坐標定義為
Nx2數組。（如果您稍后進入 3-D 點，這將有助于擴展性）將中間變量
distance,angle,定義force為NxN表示成對交互的數組
麻木的事情要知道：
如果數組具有相同的形狀（或一致的形狀，這是一個重要的話題......），您可以在數組上調用大多數數值函數
meshgrid幫助您生成對數組進行變形Nx2以計算NxN結果所需的數組索引和一個切線音符（哈哈）
arctan2()計算一個有符號的角度，所以你可以繞過復雜的“哪個象限”邏輯
例如，你可以做這樣的事情。注意點get_dist之間get_angle的算術運算發生在最底部的維度：
import numpy as np
# 2-D locations of particles
points = np.array([[1,0],[2,1],[2,2]])
N = len(points) # 3
def get_dist(p1, p2):
r = p2 - p1
return np.sqrt(np.sum(r*r, axis=2))
def get_angle(p1, p2):
r = p2 - p1
return np.arctan2(r[:,:,1], r[:,:,0])
ii = np.arange(N)
ix, iy = np.meshgrid(ii, ii)
dist = get_dist(points[ix], points[iy])
angle = get_angle(points[ix], points[iy])
# ... compute force
# ... apply the force, etc.
對于上面顯示的示例 3 點向量：
In [246]: dist
Out[246]:
array([[0. , 1.41421356, 2.23606798],
[1.41421356, 0. , 1. ],
[2.23606798, 1. , 0. ]])
In [247]: angle / np.pi # divide by Pi to make the numbers recognizable
Out[247]:
array([[ 0. , -0.75 , -0.64758362],
[ 0.25 , 0. , -0.5 ],
[ 0.35241638, 0.5 , 0. ]])
月關寶盒
TA貢獻1772條經驗 獲得超5個贊
這是每個時間步只有一個循環的一次嘗試，它應該適用于任意數量的維度，我也用 3 進行了測試：
from matplotlib import pyplot as plt
import numpy as np
fig, ax = plt.subplots()
N = 4
ndim = 2
masses = np.ones(N)
charges = np.array([-1, 1, -1, 1]) * 2
# loc_arr = np.random.rand(N, ndim)
loc_arr = np.array(((-1,0), (1,0), (0,-1), (0,1)), dtype=float)
speed_arr = np.zeros((N, ndim))
# compute charge matrix, ie c1 * c2
charge_matrix = -1 * np.outer(charges, charges)
time = np.linspace(0, 0.5)
dt = np.ediff1d(time).mean()
for i, t in enumerate(time):
# get (dx, dy) for every point
delta = (loc_arr.T[..., np.newaxis] - loc_arr.T[:, np.newaxis]).T
# calculate Euclidean distance
distances = np.linalg.norm(delta, axis=-1)
# and normalised unit vector
unit_vector = (delta.T / distances).T
unit_vector[np.isnan(unit_vector)] = 0 # replace NaN values with 0
# calculate force
force = charge_matrix / distances**2 # norm gives length of delta vector
force[np.isinf(force)] = 0 # NaN forces are 0
# calculate acceleration in all dimensions
acc = (unit_vector.T * force / masses).T.sum(axis=1)
# v = a * dt
speed_arr += acc * dt
# increment position, xyz = v * dt
loc_arr += speed_arr * dt
# plotting
if not i:
color = 'k'
zorder = 3
ms = 3
for i, pt in enumerate(loc_arr):
ax.text(*pt + 0.1, s='{}q {}m'.format(charges[i], masses[i]))
elif i == len(time)-1:
color = 'b'
zroder = 3
ms = 3
else:
color = 'r'
zorder = 1
ms = 1
ax.plot(loc_arr[:,0], loc_arr[:,1], '.', color=color, ms=ms, zorder=zorder)
ax.set_aspect('equal')
上面的示例生成，其中黑色和藍色點分別表示開始和結束位置：
當電荷相等時charges = np.ones(N) * 2，系統對稱性被保留并且電荷排斥：
最后是一些隨機的初始速度speed_arr = np.random.rand(N, 2)：
編輯
對上面的代碼做了一些小改動，以確保它是正確的。（我在合力上遺漏了 -1，即 +/+ 之間的力應該是負數，并且我總結了錯誤的軸，為此道歉。現在在 的情況下masses[0] = 5，系統正確發展：
人到中年有點甜
TA貢獻1895條經驗 獲得超7個贊
經典的方法是計算系統中所有粒子的電場。假設您有 3 個帶正電荷的帶電粒子：
particles = np.array([[1,0,0],[2,1,0],[2,2,0]]) # location of each particle
q = np.array([1,1,1]) # charge of each particle
計算每個粒子位置的電場的最簡單方法是 for 循環：
def for_method(pos,q):
"""Computes electric field vectors for all particles using for-loop."""
Evect = np.zeros( (len(pos),len(pos[0])) ) # define output electric field vector
k = 1 / (4 * np.pi * const.epsilon_0) * np.ones((len(pos),len(pos[0]))) * 1.602e-19 # make this into matrix as matrix addition is faster
# alternatively you can get rid of np.ones and just define this as a number
for i, v0 in enumerate(pos): # s_p - selected particle | iterate over all particles | v0 reference particle
for v, qc in zip(pos,q): # loop over all particles and calculate electric force sum | v particle being calculated for
if all((v0 == v)): # do not compute for the same particle
continue
else:
r = v0 - v #
Evect[i] += r / np.linalg.norm(r) ** 3 * qc #! multiply by charge
return Evect * k
# to find electric field at each particle`s location call
for_method(particles, q)
此函數返回與輸入粒子數組具有相同形狀的向量數組。要找到每個上的力，您只需將此向量乘以q電荷數組。從那里開始，您可以使用您最喜歡的 ODE 求解器輕松找到您的加速并集成系統。
性能優化和準確性
因為方法是最慢的方法?？梢詢H使用線性代數來計算該場，從而顯著提高速度。以下代碼對這個問題非常有效的 Numpy 矩陣“單線”（幾乎是單線）：
def CPU_matrix_method(pos,q):
"""Classic vectorization of for Coulomb law using numpy arrays."""
k = 1 / (4 * np.pi * const.epsilon_0) * np.ones((len(pos),3)) * 1.602e-19 # define electric constant
dist = distance.cdist(pos,pos) # compute distances
return k * np.sum( (( np.tile(pos,len(pos)).reshape((len(pos),len(pos),3)) - np.tile(pos,(len(pos),1,1))) * q.reshape(len(q),1)).T * np.power(dist,-3, where = dist != 0),axis = 1).T
請注意，此代碼和以下代碼還返回每個粒子的電場矢量。
如果您使用 Cupy 庫將其卸載到 GPU 上，您可以獲得更高的性能。以下代碼幾乎與CPU_matrix_method相同，我只是稍微擴展了單行代碼，以便您可以更好地看到發生了什么：
def GPU_matrix_method(pos,q):
"""GPU Coulomb law vectorization.
Takes in numpy arrays, performs computations and returns cupy array"""
# compute distance matrix between each particle
k_cp = 1 / (4 * cp.pi * const.epsilon_0) * cp.ones((len(pos),3)) * 1.602e-19 # define electric constant, runs faster if this is matrix
dist = cp.array(distance.cdist(pos,pos)) # could speed this up with cupy cdist function! use this: cupyx.scipy.spatial.distance.cdist
pos, q = cp.array(pos), cp.array(q) # load inputs to GPU memory
dist_mod = cp.power(dist,-3) # compute inverse cube of distance
dist_mod[dist_mod == cp.inf] = 0 # set all infinity entries to 0 (i.e. diagonal elements/ same particle-particle pairs)
# compute by magic
return k_cp * cp.sum((( cp.tile(pos,len(pos)).reshape((len(pos),len(pos),3)) - cp.tile(pos,(len(pos),1,1))) * q.reshape(len(q),1)).T * dist_mod, axis = 1).T
關于上述算法的準確性，如果你計算粒子陣列上的 3 種方法，你會得到相同的結果：
[[-6.37828367e-10 -7.66608512e-10 0.00000000e+00]
[ 5.09048221e-10 -9.30757576e-10 0.00000000e+00]
[ 1.28780145e-10 1.69736609e-09 0.00000000e+00]]
關于性能，我在 2 到 5000 個帶電粒子的系統上計算了每種算法。此外，我還包括了 for_method 的 Numba 預編譯版本，以使 for-loop 方法具有競爭力：
我們看到 for-loop 執行非常需要超過 400 秒來計算具有 5000 個粒子的系統。放大到底部：
這表明解決這個問題的矩陣方法要好幾個數量級。準確地說，Numba for-loop 對 5000 個粒子的評估需要 18.5 秒，CPU 矩陣需要 4 秒（比 Numba 快 5 倍），GPU 矩陣*需要 0.8 秒（比 Numba 快 23 倍）。較大的陣列顯示出顯著差異。
* 使用的 GPU 是 Nvidia K100。
添加回答
舉報
0/150
提交
取消