Velocity Verlet integration producing massive results python

Question

As part of a vector class/ "physics engine" I am making, I am trying to implement a velocity verlet algorithm to solve Newton's equations of motion. However, the algorithm that I have implemented is producing massive values, as you can see below. i have tried what I believe to be every option on my own, and was wondering if someone could help. Thanks.

Xpos = 0 OldPos = 1 dt = 0.02 Xaccel = 9.81 def doVerletX(currentPos, oldPos, accel, dt): newPos = (currentPos + (currentPos-oldPos) + (accel*dt*dt)) return newPos for a in range (1,10,1): newPos = doVerletX(Xpos,OldPos,dt,Xaccel) print(newPos) OldPos = Xpos dt+=dt

The output to the above code was:

0.9247220000000003 3.8494440000000005 7.698888000000001 15.397776000000002 30.795552000000004 61.59110400000001 123.18220800000002 246.36441600000003 492.72883200000007

I know that this is obviously wrong and was wondering what to do to fix it?

dt is delta t? why are you doing dt+=dt? Also I think your parameter are in wrong order when calling the funciton. In function definition dt is 4th argument, during call you are passing dt as 3rd argument. — Anonta
– Anonta, Commented Sep 29, 2017 at 11:10

Lutz Lehmann · Accepted Answer · 2017-10-03 15:12:17Z

Already mentioned was the dt+=dt problem that should be t+=dt.

In the time propagation, you also need to shift the newPos value to Xpos.

The interface definition for doVerlet has dt as last argument, you have to use it that way also in the function call.

Xpos = 0 OldPos = 1 t=0 dt = 0.02 Xaccel = 9.81 def doVerletX(currentPos, oldPos, accel, dt): newPos = (currentPos + (currentPos-oldPos) + (accel*dt*dt)) return newPos for a in range (1,10,1): newPos = doVerletX(Xpos,OldPos,Xaccel,dt) t+=dt print(t,newPos) OldPos, Xpos = Xpos, newPos

Using that code gives the result

(0.02, -0.996076) (0.04, -1.988228) (0.06, -2.976456) (0.08, -3.960760) (0.10, -4.941140) (0.12, -5.917596) (0.14, -6.890128) (0.16, -7.858736) (0.18, -8.823420)

which is believable as the initial velocity is -50 m/s and the acceleration towards the positive.

To get the correct initial velocity and thus exact solution one would have to solve

x(t)=x0+v0*t+0.5*g*t^2

for v0 at t=-0.02 using x(0)=x0=0 and x(-0.02)=1 giving

v0=-(1-0.02^2/2*9.81)/0.02=-49.9019 and x(0.18)=-8.82342

which is exactly the value computed above. This is as it should be with an order 2 method.

In real life, that object would have traveled -9 units but the verlet algorithm here says-8.82 and they is an error of approximately 2% and so I was wondering if this error would cumulatively grow, and if so make the simulation horrifically inaccurate as the error of 2% is after only 10 iterations.
The formula for the position is x(t)=-50*t+0.5*9.81*t^2 which at t=0.18 gives -8.841078 which is not that far off.
How does the velocity equal -50, because when you do currentPos - oldPos you are not dividing by dt. Also, what does the OldPos, Xpos = Xpos, newPos part mean (I know you are equating things, but I am not sure which things).
(x(0)-x(-0.02))/0.02 = -50. Of course that is only an estimate by the mean value. I added the completely exact computation of the exact solution to the answer.

Collectives™ on Stack Overflow

Velocity Verlet integration producing massive results python

1 Answer 1

4 Comments

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

4 Comments

Related