Электри c сила между частицами с использованием numpy массивов

Question

Я пытаюсь смоделировать частицу, летящую на другой частице во время электрического отталкивания (или притяжения), называемого рассеянием Резерфорда. Мне удалось смоделировать (несколько) частиц с помощью циклов и списков python. Однако теперь я хочу использовать массивы numpy. Модель будет использовать следующие шаги:

1. Для всех частиц:
   1. Рассчитать радиальное расстояние со всеми другими частицами
   2. Рассчитать угол со всеми другими частицами
   3. Рассчитать силу netto в направлении x и y
2. Создать матрицу с netto xForce и yForce для каждой частицы
3. Создать ускорение (также компоненты x и y ) матрица a = F / масса
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. (Это поможет с расширяемостью, если вы перейдете к трехмерным точкам позже)
• определите промежуточные переменные distance, angle, force как NxN массивы для представления парных взаимодействий

Numpy о чем нужно знать:

• Вы можете вызывать большинство числительных c функций для массивов, если массивы имеют одинаковую форму (или соответствующие формы, что является нетривиальным topi c ...)
• meshgrid помогает генерировать индексы массивов, необходимые для сдвига массива Nx2 массивов для вычисления NxN результатов
• и тангенциальной ноты (ха-ха ) arctan2() вычисляет угол со знаком, так что вы можете обойти сложный «какой квадрант» logi c

Например, вы можете сделать что-то подобное. Обратите внимание, что в get_dist и get_angle арифметические операции c между точками выполняются в самом нижнем измерении:

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.        ]])

Answer 2

Вот один go только с al oop для каждого временного шага, и он должен работать для любого количества измерений, я тоже проверил с 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 в результирующей силе, ie. Сила между + / + должна быть отрицательной, и я суммировал неправильную ось, извиняюсь за это. Теперь в случаях, когда masses[0] = 5, система развивается правильно :