RK4 Orbit Integration

Integrate a circular two-body orbit with the fourth-order Runge-Kutta method and check energy drift over one period.

Problem

The classical Runge-Kutta method of order four (RK4) advances an ordinary differential equation by combining four slope evaluations per step, achieving error that shrinks as the fourth power of the step size. We integrate a body in a central gravitational field and track how well it conserves energy.

Equations

State is position and velocity in 2D. The acceleration points inward as -r / |r|^3 (unit gravitational parameter). Starting on a circular orbit, an ideal integrator returns the body to its start after one period.

import numpy as np
def deriv(s):
    x,y,vx,vy=s; r=(x*x+y*y)**1.5
    return np.array([vx,vy,-x/r,-y/r])
s=np.array([1.,0.,0.,1.])   # circular orbit, period 2*pi
dt=0.01
def energy(s): x,y,vx,vy=s; return 0.5*(vx*vx+vy*vy)-1/np.hypot(x,y)
E0=energy(s)
for _ in range(int(2*np.pi/dt)):
    k1=deriv(s); k2=deriv(s+dt/2*k1); k3=deriv(s+dt/2*k2); k4=deriv(s+dt*k3)
    s=s+dt/6*(k1+2*k2+2*k3+k4)
print('pos',np.round(s[:2],4),'energy drift',round(energy(s)-E0,3e-1*0+8))

Result

After one full period the body returns close to its starting point and the energy drift is tiny, on the order of 1e-8 for this step size. RK4's fourth-order accuracy means cutting the step in half reduces error by a factor of sixteen. It is not symplectic, so over very many orbits energy slowly drifts; long-term celestial integrations use symplectic methods instead.

• RK4 balances accuracy and simplicity, the default for smooth non-stiff problems.
• Adaptive step-size variants (RK45) adjust dt to hold a target error automatically.
• Particle-orbit and gyrocenter integrators in Kronos plasma tools use RK-family or symplectic pushers depending on the time scale.