Реализация алгоритма Верле в Python

Question

У меня проблема с реализацией алгоритма верлета в Python. Я попробовал этот код:

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

При обращении за помощью в 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 значения в списке, так как они будут рассчитываться по другим параметрам с использованием соответствующего уравнения силы.