Verlet Integration of an Orbit

Simulate a planet around a star with velocity Verlet and see why symplectic integrators conserve energy over long runs.

Why not RK4 here

RK4 is accurate short-term but not symplectic: over millions of orbital steps its energy slowly drifts. Verlet integration is symplectic - it conserves a nearby shadow energy exactly - so orbits stay bounded for very long simulations.

Velocity Verlet

• x_new = x + v dt + 0.5 a dt^2
• compute a_new from the new position
• v_new = v + 0.5 (a + a_new) dt

Gravity

For a unit-mass planet around a unit-mass star at the origin, acceleration is a = -G M r / |r|^3. Start with a position and a perpendicular velocity to get an ellipse.

import numpy as np
def accel(r): return -r/np.linalg.norm(r)**3
r=np.array([1.0,0.0]); v=np.array([0.0,1.0]); dt=0.01
a=accel(r); E0=0.5*v@v-1/np.linalg.norm(r)
for n in range(100000):
    r=r+v*dt+0.5*a*dt*dt
    anew=accel(r); v=v+0.5*(a+anew)*dt; a=anew
E=0.5*v@v-1/np.linalg.norm(r)
print('energy drift:',round(abs(E-E0),6))  # stays tiny

What to observe

The orbit closes cleanly and the total energy oscillates within a tiny band instead of drifting, even after 100000 steps. An RK4 run of the same length would show a visible spiral in or out. This bounded-energy property is why molecular dynamics and n-body codes almost always use Verlet-family integrators.