python - Python 中的 Verlet 集成导致粒子逃跑

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)

我使用了 velocity Verlet 中描述的算法

Answer 1

加速度不会像速度和位置那样随着时间的推移而累积:它应该在过程的每个步骤中从头开始计算。所以，

1. 从 acc 和 verlet 的参数列表中删除 a
2. 在 acc 开头用零初始化 a。
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