1

我的闲散代码（应该）解决两个物体的运动方程，但结果是粒子运行方式，我无法找到错误在哪里

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中描述的算法

4

1 回答 1

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
于 2016-07-06T14:32:51.123 回答