Интеграция Verlet в Python приводит к тому, что частицы убегают

Question

Мой следующий код (должен был) решить уравнение движения для двух тел, но в результате частицы бегут, а я не смог найти, где ошибка

import numpy as np
import matplotlib.pyplot as plt

plt.style.use('ggplot')

DIM = 2
N = 2
ITER = 1000

def acc(r, v, a):
    for i in range(N - 1):
        for j in range(i+1, N):
            r2 = 0.0
            rij = [0.0, 0.0]
            for k in range(DIM):
                rij[k] = r[i][k] - r[j][k]
                r2 += rij[k] * rij[k]
            fg = -1.0 /np.power(np.sqrt(r2), 3)
            for k in range(DIM):
                a[i][k] += fg * rij[k]
                a[j][k] -= fg * rij[k]
    return a

def verlet(r, v, a, dt):
    for i in range(N):
        for k in range(DIM):
            v[i][k] += 0.5 * a[i][k] * dt 
            r[i][k] += v[i][k] * dt
    a = acc(r, v, a)
    for i in range(N):
        for k in range(DIM):
            v[i][k] += 0.5 * a[i][k] * dt
    return [r,v]


r = np.zeros((N, DIM))
v = np.zeros((N ,DIM))
a = np.zeros((N, DIM))  

r[0] = [0.5,0.0]
v[0] = [0.0,1.0]

r[1] = [-0.5,0.0]
v[1] = [0.0,-1.0]

dt = 0.01
plt.ion()
for i in range(ITER):
    r,v = verlet(r, v, a, dt)
    plt.scatter(r[0][0], r[0][1])
    plt.scatter(r[1][0], r[1][1],color='yellow')
    plt.pause(0.00005)

И я использовал алгоритм, описанный в скорости Verlet

Answer 1

Ускорение не накапливается с течением времени, как скорость и положение: его нужно вычислять с нуля на каждом этапе процесса. Так,

1. Удалить a из списка аргументов обоих acc а также verlet
2. инициализировать a нулями в начале acc,
3. Переместить звонок a = acc(r, v) к началу verlet,

Вы увидите, что статьи вращаются вокруг друг друга, как и ожидалось.

Не имеет отношения к вашей проблеме, но важно для эффективного использования NumPy: избавьтесь от петель над векторными координатами, оставьте для NumPy сложение и вычитание векторов. Я переписал acc метод таким образом:

def acc(r, v):
    a = np.zeros((N, DIM))  
    for i in range(N - 1):
        for j in range(i+1, N):
            rij = r[i, :] - r[j, :]
            r2 = np.dot(rij, rij)
            fg = -1.0 /np.power(np.sqrt(r2), 3)
            a[i, :] += fg * rij
            a[j, :] -= fg * rij
    return a