python - Verlet算法的Python实现

Question

我在 Python 中实现 Verlet 算法时遇到问题。我尝试了这段代码:

x[0] = 1
v[0] = 0
t[0] = 0
a[0] = 1
for i in range (0, 1000):
    x[i+1] = x[i] - v[i] * dt + (a[i] * (dt**2) * 0.5)
    v[i] = (x[i+1] - x[i-1]) * 0.5 * dt
    t[i+1] = t[i] + dt

但它无法正常工作。怎么了？ 我正在寻找 Verlet 算法的通用代码。

Answer 1

您的实现存在几个问题，

1. 在步骤 0 中，您尝试访问不存在的索引 -1，并且索引 i+1 也不存在最后一个索引。
2. 加速度仅在第一步定义，其他步骤不定义。
3. 计算下一个位置时存在符号错误。

这样的事情应该可以解决它，

import numpy as np
N = 1000
dt  = 0.1 
x = np.zeros(N)
v = np.zeros(N)
t = np.arange(0,(N+0.5)*dt, dt)
a = np.ones(N)*1.0  # initial condition
x[[0,1]] = 1
v[1] = v[0]+a[0]*dt

for i in range (1,N-1):
    x[i+1] = x[i]+v[i]*dt+(a[i]*(dt**2)*0.5)
    v[i+1] = v[i] + a[i]*dt

但这不会非常有效，最好使用scipy.integrate.ode来求解底层的运动微分方程。