python - 使用 numpy 阵列的粒子之间的电力

Question

我试图模拟一个粒子在经历电排斥(或吸引力)的同时向另一个粒子飞行，称为卢瑟福散射。我已经成功地使用 for 循环和 python 列表模拟了(一些)粒子。但是，现在我想改用 numpy 数组。该模型将使用以下步骤:

1. 对于所有粒子:
   1. 计算与所有其他粒子的径向距离
   2. 计算与所有其他粒子的角度
   3. 计算 x 方向和 y 方向的净力
2. 使用 netto xForce 和 yForce 为每个粒子创建矩阵
3. 通过 a = F/mass 创建加速度(也是 x 和 y 分量)矩阵
4. 更新速度矩阵
5. 更新位置矩阵

我的问题是 我不知道如何使用 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, yforce

def 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 parameters
t = 0
N = 3
masses = 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

Answer 1

用数组模拟这个问题的一种方法可能是:

• 将点坐标定义为 Nx2 数组。 (如果您稍后进入 3-D 点，这将有助于扩展)
• 将中间变量distance、angle、force定义为NxN数组来表示成对交互<

Numpy 事情要知道:

• 如果数组具有相同的形状(或一致的形状，这是一个不平凡的话题...)，您可以在数组上调用大多数数值函数
• 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.        ]])