CODEWITHSUNDEEP
pythonphysics
Bob
Bob

Reputation: 51

Velocity Verlet integration producing massive results python

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?

Upvotes: 0

Views: 1263

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 26040

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.

Upvotes: 1

Related Questions