Why a worse method gives a better orbit
RK4 is fourth order. Velocity Verlet is second. Over a long orbit Verlet wins comprehensively, and the reason has nothing to do with accuracy.
Everything so far has ranked methods by order of accuracy. Higher order, smaller error, better method. That ranking is about to fail badly.
Take a harmonic oscillator — the cleanest conservative system there is — and integrate it for a long time with RK4 (fourth order) and velocity Verlet (second order). Then look not at the trajectory but at the energy, which the true dynamics conserves exactly.
Integrating a harmonic oscillator for 4000 time units at h = 0.1. What happens to the energy under RK4 (order 4) versus velocity Verlet (order 2)?
Switch the view to energy and lengthen the run. The vertical axis is the relative change in total energy; the exact answer is a flat line at zero.
Two qualitatively different behaviours. RK4's energy error is smaller at first and then keeps growing. Verlet's is larger at first and then stops.
That distinction — bounded versus secular — matters more than order for any long integration, and it is the difference between a solar-system model that stays a solar system and one where the planets slowly spiral into the sun.
One line of code
Here is what makes Verlet different. Forward Euler, on a system with position and velocity :
const vNext = v + h * a(q); // both updates use the OLD position
const qNext = q + h * v; // and the OLD velocitySymplectic Euler:
const vNext = v + h * a(q); // update velocity first
const qNext = q + h * vNext; // then use the NEW velocity for positionThat is the entire difference. One variable is updated before the other instead of both being updated from the old state. The order of accuracy is unchanged — both are first order — and yet one has bounded energy error forever and the other does not.
Velocity Verlet is the second-order member of the same family:
const a = accel(q);
const qNext = q + h * v + 0.5 * h * h * a; // half-step-aware position
const aNext = accel(qNext);
const vNext = v + 0.5 * h * (a + aNext); // average the accelerationsForward Euler computes both updates from the old state. Symplectic Euler updates the velocity first, then uses that already-updated velocity to advance the position.
// forward Euler // symplectic Euler
v1 = v0 + h*a(q0); v1 = v0 + h*a(q0);
q1 = q0 + h*v0; q1 = q0 + h*v1; // ← uses v1, not v0It matters because the reordered version is a symplectic map: it exactly preserves the phase-space area element that Hamiltonian flow preserves. Forward Euler expands that area every step, which shows up as energy growing without bound. Both methods are first order, so accuracy does not distinguish them at all — the structure does.
Why bounded, not just small
The explanation is the deepest idea in this path, and it reframes what a numerical method is.
A symplectic integrator does not approximately solve your problem. It exactly solves a nearby one. There exists a modified Hamiltonian — the shadow Hamiltonian —
such that the numerical trajectory lies (to exponentially small error, over exponentially long times) on exact orbits of .
Everything follows from that. is conserved exactly because it is a Hamiltonian and the method is its exact flow. And differs from the true by . So the true energy along the numerical trajectory stays within of its initial value — forever. It oscillates as the trajectory moves around the shadow orbit, but it cannot wander off, because it is pinned to a conserved quantity.
RK4 has no shadow Hamiltonian. It is not a symplectic map, so there is nothing pinning its energy, and the per-step error — though far smaller — accumulates in one direction.
That is the whole story:
- RK4: small error, no constraint → error accumulates linearly in time.
- Verlet: larger error, but bounded by a conserved quantity → error never grows.
Over ten periods RK4 wins easily. Over ten million, it is not close.
The orbit you can see it in
The harmonic oscillator is clean but forgiving. A Kepler orbit with real eccentricity is where the difference becomes visceral: the perihelion passage is fast and tightly curved, and a fixed-step method has to survive it.
Verlet's ellipse precesses slightly — a phase error, which is second order and expected — but it stays an ellipse of the right size. RK4's ellipse shrinks. Given long enough, RK4's planet falls into the star, for no physical reason whatsoever.
Hamiltonian flow has a property stronger than energy conservation: it preserves the symplectic 2-form . Geometrically, it preserves oriented area in each conjugate plane — Liouville's theorem, and the reason phase-space volume is conserved.
A numerical method is symplectic when its update map preserves exactly, not approximately. For a one-degree-of-freedom system that reduces to the Jacobian of the step map having determinant exactly 1. Check it for symplectic Euler on the oscillator and you get exactly 1, with no remainder. It is an algebraic identity, not an approximation.
Backward error analysis then supplies the shadow Hamiltonian, and the bounded energy error follows.
The constraints are real, and worth stating plainly:
- Fixed step size. Varying breaks the conjugacy to a single shadow Hamiltonian and reintroduces drift. Adaptive symplectic integration is genuinely hard, and the usual answer is to change time variable instead (regularisation) rather than to change .
- Separable Hamiltonians. gives explicit methods. Non-separable systems need implicit ones, which cost far more per step.
- No dissipation. These methods encode conservative structure. Add friction and the premise is gone.
- Trajectory accuracy is still only order . Verlet's phase error grows linearly in time. If you need to know where the planet is on a specific date, high order still matters; if you need to know the orbit is still an orbit, structure matters.
Time-reversibility is a related and often sufficient property: velocity Verlet is reversible, and reversibility alone rules out the monotone drift that dissipation-free non-symplectic methods exhibit.
Symplectic integration is one instance of geometric numerical integration: build the discretisation so that structure is preserved exactly, and let accuracy be secondary.
- Variational integrators discretise the action and derive the update from a discrete Euler–Lagrange equation. Symplecticity and a discrete Noether theorem come out automatically — conserved momenta stay conserved to machine precision. This is arguably the cleaner derivation, and it extends to constrained and non-smooth mechanics.
- Lie-group integrators keep the state on its manifold — a rotation stays a rotation, rather than drifting off and needing re-orthonormalisation.
- Structure-preserving PDE discretisations (mimetic finite differences, discrete exterior calculus) preserve or discrete vorticity exactly, which is why long-run magnetohydrodynamics and climate models care.
- Structure-preserving machine learning — Hamiltonian and Lagrangian neural networks — builds the same constraint into learned dynamics, so a learned model conserves energy by construction instead of by hoping the training data covered it.
The through-line: a discretisation that respects the structure of the problem beats a more accurate one that does not, whenever you integrate long enough.
What to carry forward
- Order of accuracy is not a total ordering of methods. It measures short-term error and says nothing about long-term qualitative behaviour.
- Ask what the system conserves, and whether your method knows about it.
- Bounded error and small error are different goals. Long integrations want the first.
- The next chapter takes this further: methods that abandon time-stepping altogether.