Verlet algorithm implementation in Python
I have problem with implementation of verlet algorithm in Python. I tried this code:
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
But it is not working properly. What is wrong? I am looking for general code for Verlet algorithm.
Answer
Your question isn't very clear, but there are a few sources of error in your code.
Eg, for i > 0
x[i+1] = x[i]-v[i]*dt+(a[i]*(dt**2)*0.5)
tries to use the value of v[i], but that element doesn't exist yet in the v list.
To give a concrete example, when i = 1, you need v[1], but the only thing in the v list at that stage is v[0]; v[1] isn't computed until the following line.
That error should cause the Python interpreter to display an error message:
IndexError: list index out of range
When asking for help on Stack Overflow please post the error messages with your question, preferably starting with the line:
Traceback (most recent call last):
That makes it a lot easier for people reading your question to debug your code.
FWIW, on the first iteration of your for loop, when i == 0, x[i+1] and x[i-1] both refer to the same element, since there are two elements in the x list at that stage, so x[-1] is x[1].
Also, it's strange that you are storing your t values in a list. You don't need to do that. Just store t as a simple float value and increment it by dt each time through the loop; note that t itself isn't used in any of the calculations, but you may want to print it.
Your formulas don't appear to match those given on the Wikipedia page, either for the Basic Störmer–Verlet, or the Velocity Verlet. For example, in the code line I quoted earlier you are subtracting v[i]*dt but you should be adding it.
Perhaps you should consider implementing the related Leapfrog integration method. The synchronised version of Leapfrog is easy to code, and quite effective, IME.
From the Wikipedia link:
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
Generally, it's not necessary to actually store the a values in a list, since they will be calculated from the other parameters using the relevant force equation.