Verlet, Leapfrog & RK4: Choosing a Numerical Integrator

What an integrator actually does

Nearly every simulation on this site solves the same problem: you know the forces acting on a body right now, and you want its position a moment later. Newton gives you a second-order ODE, a = F(x)/m, and a numerical integrator turns it into a rule for stepping (x, v) forward by a finite time h. Which rule you choose decides whether a planet stays in orbit for an hour of wall-clock time or spirals into the sun by the third minute.

Two properties matter and they are not the same thing. Accuracy is how small the error is after one step (local truncation error) or after a fixed time (global error). Stability is whether the error stays bounded as you keep stepping. A method can be accurate and unstable, or inaccurate and beautifully stable — Velocity Verlet is the second kind, and for physics that is usually what you want.

Explicit Euler: the one that lies

The cheapest scheme evaluates the force once and moves everything with it:

// explicit (forward) Euler — do not use for orbits
v += a(x) * h;
x += v * h;        // v here is the OLD v if you write it in this order

Its local error is O(h²) and its global error O(h) — first order. Worse, for an oscillator it is unconditionally unstable: the amplitude grows every single step, no matter how small h is. Simulate a spring or a circular orbit with explicit Euler and the energy climbs monotonically; the planet spirals outward. Shrinking the step only slows the disaster, it does not prevent it.

The nearly free fix is semi-implicit (symplectic) Euler: update the velocity first, then move the position with the new velocity. Same cost, one line reordered — and now the method is symplectic and the energy stops drifting. This is the difference between the two Euler variants that most tutorials never mention.

// semi-implicit Euler — stable, 1st order, one force evaluation
v += a(x) * h;
x += v * h;        // v is the NEW v — this ordering is the whole trick

Symplectic, and why it beats raw accuracy

Hamiltonian systems — gravity, springs, pendulums, molecular dynamics — conserve phase-space volume. A symplectic integrator conserves a discrete analogue of it, and as a consequence it exactly conserves a slightly perturbed "shadow" energy that stays within O(h²) of the true energy forever. Practically: the energy of a symplectic run oscillates around the correct value instead of drifting away from it. A non-symplectic method of higher order will have a smaller error over one orbit but a secular drift over ten thousand.

Velocity Verlet and Leapfrog

Velocity Verlet is the workhorse of molecular dynamics and of every game engine's rigid-body solver. It is symplectic, time-reversible, second order, and needs exactly one force evaluation per step if you cache the acceleration:

// Velocity Verlet — 2nd order, symplectic, 1 force eval/step
x += v * h + 0.5 * aOld * h * h;
const aNew = a(x);                 // the only force evaluation
v += 0.5 * (aOld + aNew) * h;
aOld = aNew;

Leapfrog is the same method wearing different clothes: positions live on integer steps, velocities on half-steps, and the two "leap" over each other. It is algebraically equivalent to Velocity Verlet for constant h (identical trajectories, up to rounding) and is the standard choice in astrophysical N-body codes. Its one annoyance is that positions and velocities are never known at the same instant, so kinetic energy needs an extra half-kick to be reported honestly — and, because the equivalence relies on a constant step, naively changing h mid-run breaks the symplectic property.

There is also position Verlet ("Störmer"), which stores the previous position instead of a velocity: x_next = 2·x − x_prev + a·h². It is what cloth and rope solvers use, because a positional constraint is applied by simply moving a point — the velocity is implied by the difference of positions and updates itself for free.

RK4: accurate, not symplectic

Classical fourth-order Runge–Kutta samples the derivative four times per step and blends them. Its global error is O(h⁴), so it is dramatically more accurate per step than Verlet — but it costs four force evaluations, and it is not symplectic. Run RK4 on a closed orbit and the energy decays slowly and monotonically: the planet spirals in. For a chaotic system integrated over a short horizon (the double pendulum, the Lorenz attractor) that is fine and the accuracy is worth it; for a solar system integrated over a million years it is a disaster.

const k1 = f(t, y);
const k2 = f(t + h/2, y + h/2 * k1);
const k3 = f(t + h/2, y + h/2 * k2);
const k4 = f(t + h,   y + h   * k3);
y += h/6 * (k1 + 2*k2 + 2*k3 + k4);   // 4 evaluations, O(h⁴), not symplectic

Choosing the step size

The step is bounded by the fastest thing in your system. For an oscillator of angular frequency ω, an explicit second-order method needs roughly h·ω < 2 to remain stable, and you want far below that — a rule of thumb is 20 to 50 steps per period of the stiffest mode before the motion looks right. For particle systems the analogous bound is the CFL condition: no particle may cross more than a fraction of the interaction radius in one step.

Two practical habits. First, decouple the physics step from the frame rate: run a fixed h (say 1/240 s) inside an accumulator loop and interpolate for rendering, otherwise a dropped frame silently changes your integrator's behaviour. Second, plot the total energy. It is the cheapest possible bug detector — a monotone climb means explicit Euler somewhere, a monotone decay means a non-symplectic method or too much damping, and a bounded wobble means you are fine.

What this site uses where

orbits, N-body, springs, molecular dynamics  → Velocity Verlet / Leapfrog
cloth, ropes, soft bodies                    → position Verlet + constraints
SPH fluids                                   → Leapfrog with a CFL-limited step
double pendulum, Lorenz, chaotic systems     → RK4 (short horizon, high accuracy)

Frequently asked questions

Try it live

Everything above runs in your browser — open Double Pendulum and change the parameters while it is running. Nothing is installed, nothing is uploaded, the whole model lives in one tab.