1

我在 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 算法的通用代码。

4

2 回答 2

3

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

例如，对于i> 0

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

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

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

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

IndexError: list index out of range

在 Stack Overflow 上寻求帮助时，请在您的问题中发布错误消息，最好从以下行开始：

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 集成方法。Leapfrog 的同步版本很容易编码，而且非常有效，IME。

从维基百科链接：

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值实际存储在列表中，因为它们将使用相关的力方程从其他参数计算出来。

于 2015-03-12T13:11:10.363 回答
2

您的实施存在几个问题，

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

于 2015-03-12T13:19:05.790 回答