Part 2 of 6. A time step has no reason to conserve energy. One that keeps the symplectic structure of Hamilton's equations keeps its energy error bounded and some momenta exact, though not every symmetry.
Part 1 discretised space; this part discretises time. The exact flow of a Hamiltonian system preserves phase-space area, runs backwards, and conserves the momentum of each symmetry of the Hamiltonian. The figures below measure which of these four common integrators keep.
For a Hamiltonian H(q, p) the equations of motion are
and for H = p2/2 + V(q) they are Newton's law. The Jacobian M of the flow over any time satisfies
In one degree of freedom this is det M = 1: the flow preserves area in the (q, p) plane, which is Liouville's theorem. With more degrees of freedom the condition is stronger than det = 1, since it preserves the oriented area of every two-dimensional projection, not just total volume. A map that satisfies it is called symplectic.
Take the pendulum, H = p2/2 − cos q, and three ways to advance it by a step h:
q ← q + h · p p ← p − h · sin(q_old)
k1 = f(y), k2 = f(y + h·k1/2) k3 = f(y + h·k2/2), k4 = f(y + h·k3) y ← y + h/6·(k1+2k2+2k3+k4)
p ← p − (h/2) · sin(q) q ← q + h · p p ← p − (h/2) · sin(q)
RK4 costs four force evaluations per step, Verlet one (its closing half-kick reuses the next step's force).
Figure 1. One start (θ₀, 0), stepped in parallel; the purple curve is the true level set of H. Below each portrait is H − H(0) against step number, each on its own labelled scale. Euler spirals out. RK4's energy drifts steadily down. Verlet's oscillates in a band that does not widen. The line above runs RK4 on until its drift leaves Verlet's band.
RK4 is more accurate for a while, and then it is not. The Jacobian of a single step shows why.
Euler's step has Jacobian
which stretches area near the stable equilibrium (cos q > 0) and shrinks it near the inverted position. Verlet's step is a product of three shears, each with determinant exactly 1, so the product is symplectic for every h:
RK4's Jacobian has no such factorisation. Figure 2 computes det J − 1 for one step at every point of the pendulum's phase plane, differentiating each update exactly.
Figure 2. Colour is det J − 1 for one step on a log scale: brown where the step grows area, teal where it shrinks it, blank where it keeps area to within 10−14. Each panel spans −π ≤ q ≤ π, −2.5 ≤ p ≤ 2.5. Euler shows the stripes of h2 cos q. RK4's error is about five orders of magnitude smaller at the default step, takes both signs, and is zero only on isolated curves. Verlet is blank at every h. Click a panel, or focus it and use the arrow keys, to move the probe.
A per-step error in area compounds. On the harmonic oscillator every one of these maps is linear, so a square of initial conditions stays a parallelogram and its area can be measured exactly.
Figure 3. 32 points on the edge of a small square, carried by the harmonic oscillator H = (p2 + q2)/2 with h = 0.05; dashed is the start. The true flow rotates the square rigidly. Euler's determinant here is exactly 1 + h2, and its readout compares the measured shoelace area with (1 + h2)n. Verlet shears the square and keeps its area.
Verlet is only second-order accurate, so it cannot conserve H. Backward error analysis shows that its trajectory follows, very closely, the exact flow of a nearby shadow Hamiltonian
where, for H = p2/2 + V(q) and this kick-drift-kick ordering, the first correction is
The series diverges in general. Benettin and Giorgilli (1994) and Hairer and Lubich (1997) showed that truncating it at the optimal order gives a function the method conserves up to an error exponentially small in 1/h, over times exponentially long in 1/h. So H itself should stay in a band of width O(h2), and H + h2H2 in one of width O(h4). Figure 4 measures both, and the correction with its sign flipped.
Figure 4. One Verlet run of the pendulum from θ₀ = 1.3 to t = 400. Left: H − H(0) and Ĥ − Ĥ(0) with Ĥ = H + h2H2, drawn as min–max envelopes. Right: the width of each band against h on log–log axes, with H − h2H2 in grey. The fitted slopes are the orders of the error: 2 for H, 4 for the correct Ĥ, and 2 again for the wrong sign, whose band is twice as wide as H's. Click a point to pick its h.
Euler and RK4 have no such function: neither preserves area, and no Hamiltonian flow changes it.
The pendulum's flow is reversible: flip the sign of p, run for the same time again, and you are back at the start with p reversed. Verlet's step is symmetric (its inverse is itself with h → −h), so it retraces exactly. Symplectic Euler, a full kick and then a drift, preserves area but is not symmetric. RK4 is neither.
Figure 5. From (1.3, 0) with h = 0.2: run N steps, flip p, run N more, and measure the distance from (1.3, 0). Log scale; the dashed line is double-precision rounding, 2.2 × 10−16. Verlet misses by accumulated rounding only. Symplectic Euler comes back near the right orbit but at the wrong phase, so its miss is on the scale of the orbit itself. RK4 and Euler miss by amounts that grow with N.
Hamiltonian Monte Carlo relies on this: its acceptance rule is correct only if the integrator preserves volume and is reversible, which is why it uses Verlet under the name leapfrog (Neal 2011).
Each continuous symmetry of H comes with a conserved momentum. Verlet's kick is the exact flow of V alone and its drift the exact flow of p2/2 alone. If both are rotation invariant, each substep conserves angular momentum L exactly, and so does the step. Symplectic Euler is built from the same two flows. RK4 treats a rotated orbit exactly as it treats the original, yet loses L anyway. The planar Kepler problem, H = |p|2/2 − 1/|q|, has the rotation symmetry to test this.
Figure 6. Kepler orbits from perihelion at distance 1 to t = 1000. Top: the first ten orbits under each method, true ellipse dashed, with the largest |ΔL| over the run. Bottom: the largest relative energy error so far, on log–log axes. The two symplectic methods hold L to rounding and their energy error in flat bands. RK4 holds neither exactly; its energy error starts far lower and keeps growing, in jumps at the close passes. The readout fits that growth as a power of t and extrapolates to where it would meet Verlet's band.
On any run this figure can show, RK4 has the smaller energy error, while losing angular momentum. The bounded band wins only on runs of the length the readout reports, which is the scale of long molecular-dynamics and solar-system integrations. Higher-order symplectic methods, such as Yoshida's (1990) compositions of Verlet steps, keep the band and the exact momenta with a smaller band width.
Kepler's problem has a hidden symmetry. The Runge–Lenz vector A = p × L − q/|q| points at perihelion and is conserved, which is why the orbits close. Verlet's shadow Hamiltonian is still rotation invariant, so L survives, but its potential is no longer exactly 1/r, so A does not.
Figure 7. Direction of A, averaged over each orbit, for 16 orbits at eccentricity 0.5. The slope is the precession per orbit, fitted by least squares; the readout reruns each method at h/2. Verlet and symplectic Euler precess at nearly the same rate, and halving h divides it by about 4. Symplectic Euler's line starts offset because its orbit is shifted at first order in h. RK4 precesses too, but at a rate that falls as h4 and is flat at this scale.
A lattice gas conserves mass and momentum at every collision, yet whether its flow is isotropic depends on a rank-4 tensor the square lattice cannot make isotropic (Part 4); Part 5 repairs the square lattice with weighted diagonals. First, Part 3 follows the symmetry of Conway's Life as it runs.