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来求解基本的运动微分方程。

Answer 2

您的问题不是很清楚，但是您的代码中有一些错误源。

例如，对于i>；0

x[i+1] = x[i]-v[i]*dt+(a[i]*(dt**2)*0.5)

尝试使用v[i]的值，但该元素尚不存在于v列表中。

举一个具体的例子，当i=1时，您需要v[1]，但该阶段v列表中唯一的东西是v[0]；v[1]直到下一行才计算。

该错误应导致Python解释器显示错误消息：

IndexError: list index out of range

在寻求有关堆栈溢出的帮助时，请将错误消息与您的问题一起发布，最好从以下行开始：

Traceback (most recent call last):

这使得人们更容易阅读你的问题来调试你的代码。

FWIW，在for循环的第一次迭代中，当i==0，x[i+1]和x[i-1]都引用同一个元素时，因为在该阶段的x列表中有两个元素，所以x[-1]是x[1]。

另外，奇怪的是，您将t值存储在列表中。你不需要那么做。只需将t存储为一个简单的float值，并在每次循环中递增dt；请注意，t本身不用于任何计算，但您可能需要打印它。

您的公式似乎与维基百科页面上给出的公式不匹配，无论是针对Basic Störmer–Verlet，还是针对Velocity Verlet。例如，在我前面引用的代码行中，您要减去v[i]*dt，但您应该添加它。

也许您应该考虑实现相关的Leapfrog integration方法。同步版本的leapprog很容易编写代码，而且非常有效。

从维基百科链接：

x[i+1] = x[i] + v[i] * dt + 0.5 * a[i] * dt * dt
v[i+1] = v[i] + 0.5 * (a[i] + a[i+1]) * dt

通常，不需要将a值实际存储在列表中，因为它们将使用相关的力方程从其他参数计算。