Why your physics sim drifts, and when RK4 is the wrong fix

Put a satellite in a perfect circular orbit, step it forward with the obvious method, and watch it slowly spiral outward and escape. You did not change the physics. The gravity is exact, the initial velocity is exact. The orbit decays because of how you advanced time, and fixing it properly is more interesting than "use a better integrator."

The obvious method, and why it leaks

A simulation tracks a state s (here, position and velocity) and a function f(s) that gives its rate of change. The most obvious step is Euler's method: multiply the current rate by a small time step and add it on.

def euler_step(s, f, dt):
    return s + dt * f(s)

The problem is that f(s) is the slope at the start of the step, and you use it for the whole step. For anything curved, that single slope is wrong before you finish using it, and the errors do not cancel, they accumulate. For an orbit the bias is systematic: Euler consistently overshoots, the satellite ends each step slightly too far out and slightly too fast, and over thousands of steps the orbital energy climbs without bound. The orbit spirals open. Halve dt and the error only halves with it, because Euler's global error is proportional to the step size. To get a clean orbit you would need an absurdly tiny step.

RK4: sample the slope four times

The standard fix is fourth-order Runge-Kutta. Instead of trusting one slope, it samples four across the step, a beginning, two estimates of the middle, and an end, and takes a weighted average:

def rk4_step(s, f, dt):
    k1 = f(s)
    k2 = f(s + 0.5 * dt * k1)
    k3 = f(s + 0.5 * dt * k2)
    k4 = f(s + dt * k3)
    return s + dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6.0

The payoff is dramatic. RK4's global error is proportional to the step size to the fourth power, so halving dt cuts the error by sixteen, not two. With RK4 the spiral that wrecked the Euler orbit vanishes at a step size where Euler is still hopeless. For most simulations, this is the answer, and it is why RK4 is the workhorse of every physics engine. If you take one thing from this, it is: do not ship Euler.

The twist: RK4 is accurate, but it leaks too

Here is what the usual tutorial does not tell you. RK4 is accurate, but it is not symplectic, and for a system that conserves energy, like an orbit, accuracy is not the same as conservation.

Run an RK4 orbit for a few thousand steps and it looks perfect. Run it for ten million steps, a long planetary integration, and the energy slowly, monotonically drifts. The error per step is tiny, but it has a consistent sign, so it integrates into a trend over enough steps. RK4 does not blow up like Euler, but it does not hold the orbit's energy steady forever either. For short flights this never matters. For long-term celestial mechanics it is the whole ballgame.

The fix is not a higher-order RK. It is a different kind of method. A symplectic integrator preserves the geometric structure of Hamiltonian systems, and the simplest one is symplectic Euler: update the velocity using the old position, then update the position using the new velocity.

def symplectic_step(pos, vel, accel, dt):
    vel = vel + dt * accel(pos)   # use the old position
    pos = pos + dt * vel          # use the just-updated velocity
    return pos, vel

That one reordering changes everything. Symplectic Euler is only first-order accurate, technically worse than RK4 at any given step. But it does not conserve the true energy and instead conserves a nearby "shadow" energy, which means its energy error oscillates within a fixed band and never trends. Run it for a billion steps and the orbit is still an orbit. The velocity Verlet method (leapfrog) is the second-order version everyone actually uses, same idea.

So the honest ranking has no single winner:

• For a general ODE over a short horizon, RK4. Accurate, robust, easy.
• For a Hamiltonian system over a long horizon, a symplectic method. Lower order, but bounded energy error forever.
• For stiff systems or wildly varying dynamics, neither fixed-step method; you want adaptive step size or an implicit solver.

The reason this is worth knowing is that "use RK4" is the advice everyone gives, and it is right exactly until the moment it is quietly wrong. The way to internalize it is to integrate the same circular orbit three ways, Euler, RK4, and symplectic Euler, and plot the energy over a long run. Euler's climbs, RK4's drifts, and the first-order symplectic method, the supposedly worst one, is the only line that stays flat. That plot rearranges your intuition about what "accuracy" even means.