A step that must not be adapted

The rung below set two integrators against each other on the same orbit and found the ranking inverted from what accuracy would suggest. Runge–Kutta is far the more accurate method and its energy error grows steadily; leapfrog is cruder and its energy error oscillates inside a band and stays there. Over five billion years only one of those properties survives contact with the problem.

The band has a condition attached, and it is not the one anybody expects. It is not a condition on how small the step is. It is a condition that the step be the same at every step, and it is broken by the single most natural thing a person writing an integrator does: refining the step where the orbit is moving fast.

What the guarantee actually is #

A symplectic integrator does not conserve energy. It is worth saying that plainly, because the bounded band in the figure looks exactly like conservation and is not.

What it does is conserve something else. For a system with Hamiltonian $H$ and a fixed step $h$, the map the integrator applies is — to all orders that matter — the exact time-$h$ flow of a different Hamiltonian,

$$\tilde H = H + h^2 H_2 + h^4 H_4 + \dots,$$

called the shadow Hamiltonian. The numerical trajectory is not an approximation to the true one; it is an exact solution of a nearby problem. And exact solutions conserve their own energy exactly.

So the measured energy error is the difference between $H$ and $\tilde H$ evaluated along the trajectory. That difference is a smooth function of position and momentum, of size $h^2$, and it therefore oscillates as the body goes round its orbit and returns to the same value when the body returns. It does not grow, because there is nothing for it to grow out of.

Why changing the step destroys it #

The shadow Hamiltonian depends on $h$. Change the step and a different nearby Hamiltonian is being solved exactly.

While the step is constant, the numerical trajectory lies on a level set of $\tilde H$ and stays there for ever. Change the step at some moment, and the state — which is on a level set of $\tilde H_1$ — is now being propagated on a level set of $\tilde H_2$ passing through that same point. The value of the true energy $H$ at that point is whatever it was; the subsequent oscillation is now about a different centre.

Do that once and nothing much happens. Do it at every step, with the step responding to the body’s position, and the offsets accumulate. They accumulate as a random walk: each change contributes a small displacement whose sign depends on where in the orbit it happened, and a sum of many such displacements grows as the square root of their number rather than cancelling.

The result is the rising curve in the first figure. It is not a bug in the implementation and it is not a matter of the control being too aggressive. It is what a symplectic integrator does when it stops being symplectic, and the moment it stops is the moment the step changes.

Why the simplest scheme is the one with the property #

It is worth saying where the symplectic property comes from, because it explains why the method that has it is cruder than the one that does not.

The Hamiltonian of a gravitating system splits into two pieces, each of which is exactly solvable on its own. The kinetic part alone advances positions at constant velocity — a drift. The potential part alone changes velocities at constant position — a kick. Neither is the full problem and each is integrable in closed form.

A symplectic integrator is built by alternating exact solutions of the pieces. Half a kick, a full drift, half a kick: that is leapfrog, and every step of it is an exact flow of something, so the composition is an exact flow of something too. The map is therefore area-preserving in phase space — symplectic — by construction rather than by accident, and the shadow Hamiltonian is the thing it is the exact flow of.

The symmetry of the arrangement is what makes the error second order. Kick–drift–kick is time-reversible: run it backwards and the same sequence appears. A reversible map cannot have odd-order error terms, because those would change sign under reversal and the map does not — so the leading error is $h^2$ and the expansion of the shadow Hamiltonian contains only even powers.

That is the whole construction, and its crudeness is the point. A higher-order method built by taking a weighted combination of Taylor terms is more accurate and is not the exact flow of anything, so it has no shadow Hamiltonian and no bound. Higher-order symplectic methods exist and they are built the other way: by composing more drifts and kicks with carefully chosen coefficients, so that the composition remains an exact flow while the low-order error terms cancel.

The property belongs to the composition rather than to the accuracy, and any modification that breaks the composition breaks it — which is what a varying step does.

The arithmetic of the accumulation #

The random walk is worth putting numbers to, because the scaling is what decides whether the effect is a curiosity or a disaster.

A fixed-step method’s energy error is bounded by some $\epsilon \sim h^2$ and stays there. A varying-step method’s error is a sum of $N$ independent displacements of typical size $\delta$, where $\delta$ is how much the shadow Hamiltonian moves when the step changes by the amount it changes in one step. That sum grows as $\delta\sqrt{N}$.

Two consequences follow, and they pull in opposite directions. Because the growth is a square root rather than a linear drift, an adaptive symplectic method is still much better than a non-symplectic one of the same order, whose error grows as $N$. And because it grows at all, there is always a run length past which the fixed step wins — and that length is short. In the first figure it is a few hundred revolutions. For a solar-system integration a few hundred revolutions is a few thousand years, against a target of five billion.

The crossing point moves with the aggressiveness of the control. A step that varies by ten per cent around the orbit crosses later than one that varies by a factor of five. It always crosses.

What the step control was for #

None of this would matter if a constant step were adequate everywhere. It is not, and the reason is eccentricity.

An orbit at $e = 0.9$ moves nineteen times faster at perihelion than at aphelion, and its radius of curvature there is a hundredth of what it is at the far end. A step chosen to resolve perihelion is absurdly fine for the other ninety-five per cent of the orbit; a step chosen for the mean is nowhere near fine enough for perihelion, and the error committed in the few steps near closest approach dominates the whole revolution.

The natural response is to make the step follow the free-fall time, $h \propto r^{3/2}$ — more steps where the physics is fast and fewer where it is slow. Judged one step at a time it is unambiguously better. That is the trap: it is better at every step and ruinous over a million of them.

What is done instead #

Three things, and they are all ways of avoiding the choice rather than making it.

Take a fixed step short enough for the worst encounter that will ever happen. This is what most planetary integrations do, and it is affordable only because the planets’ orbits are nearly circular and their mutual encounters are distant. It fails the moment a small body on a crossing orbit is included.

Change the independent variable rather than the step. If the equations are rewritten so that the integration proceeds in a variable that is not time — an eccentric anomaly, or a fictitious time defined by $dt = r\,ds$ — then a constant step in the new variable is a variable step in time, automatically fine near perihelion, and the map remains symplectic in the transformed system. This is the same manoeuvre that turns the singularity of a collision orbit into a change of variable, and it is close kin to the solve that does not ask which conic it is. It is not free: the transformation has to be built into the splitting, and a time-transformed symplectic method for a general $N$-body problem is considerably harder to construct than the two-line leapfrog it replaces. Or switch integrators for the duration of the encounter. This is what long-term solar-system codes actually do. The system is split into a dominant Keplerian part, solved exactly, and a perturbation, applied as periodic kicks; that splitting is the Wisdom–Holman map and it is what makes a step of a twentieth of the shortest orbital period sufficient for the planets, whose mutual perturbations are the small terms the elements respond to. When two bodies come within a few Hill radii of each other the splitting stops being a good one, because the “perturbation” is no longer small, and the code hands that piece of the trajectory to a conventional adaptive method — a Bulirsch–Stoer or a high-order Runge–Kutta — for as long as the encounter lasts, then hands it back. What counts as an encounter is set by the radius within which the planet rather than the Sun is in charge.

The hybrid is not symplectic during the encounter and does not pretend to be. What it relies on is that encounters are rare and brief: a bounded error committed over a few hundred years out of five billion is a bounded error, and the drift-free behaviour is restored for the whole of the remaining time.

What it costs, and how it is checked #

The bookkeeping is the awkward part. Handing a subsystem to a different integrator means deciding when the encounter begins and ends, and a switch that happens at slightly the wrong moment is itself a step-size change of the kind this essay is about. Good implementations blend the two smoothly, with a changeover function that goes from zero to one over a few steps, precisely so that nothing discontinuous happens to the map.

And the check that any of it worked is the same check the rung below used, run over a longer baseline: the energy error is monitored, and a run whose energy has drifted is a run whose answer is not trusted whatever else it says.

What the long integrations report #

The five-billion-year integrations that produce the stability statistics for the solar system are all built this way, and their published diagnostics are the ones this essay is about.

They quote a relative energy error, typically a few parts in $10^{11}$ over the whole run, and they quote it as a bound rather than as an endpoint — because the whole point of the method is that the endpoint is not special. They quote the angular momentum error separately, because the splitting conserves it to machine precision by construction and a departure would indicate a coding error rather than a discretisation one. And they report how many close encounters were handed to the non-symplectic branch and for how long, because that is the part of the run whose error is not bounded by anything structural.

What none of them claims is to know where Mercury will be. A five-billion-year integration is a sample of what the system might do, not a prediction of what it will, and the reason is one the collection has already made: an exponential divergence turns a centimetre of initial uncertainty into an entire orbit in a hundred million years, so the trajectory computed is one of many consistent with the data.

The error that is not the method’s #

There is a floor under everything above and it is arithmetic rather than analysis.

A five-billion-year integration at a step of a fiftieth of Mercury’s year is about $10^{12}$ steps. Each one adds two numbers in floating point, and each addition rounds. Rounding errors are not systematic — they are as often up as down — so they accumulate as a random walk, and $10^{12}$ of them at a relative size of $10^{-16}$ each gives $10^{-16}\sqrt{10^{12}} = 10^{-10}$.

That is comparable with the energy errors the published integrations quote, which means the arithmetic and the discretisation are contributing at the same level and neither can be reduced by improving the other.

The response is to change how the addition is done. In compensated summation, the part of each addition lost to rounding is computed explicitly and carried forward into the next step, so the error no longer accumulates as a random walk but stays bounded. It costs a few extra operations per step and it removes a term that no amount of care about the physics would touch.

There is a second and subtler arithmetic effect. A position accumulated as a running sum loses precision as it grows, so a body whose coordinates are large in the chosen units is integrated less accurately than one whose coordinates are small — for no physical reason whatever. The fix is to carry the drift as an increment rather than as an absolute position, which is a change of variables made purely to keep the numbers near unity.

Neither of those is physics and both change the answer, which is the usual position for a computation whose interesting property is a bound on an accumulated quantity: the bound has to hold against the arithmetic as well as against the equations.

The general shape of the mistake #

It is worth naming the pattern, because it recurs and it is not about integrators.

A method is chosen for a structural property — a conserved quantity, a bounded error, an exact symmetry. The property depends on a condition that is not the obvious one. Somebody then improves the method along the obvious axis, the condition is violated as a side effect, and the improvement is measurably better on every short test and worse on the long one the method was chosen for.

Here the obvious axis is accuracy and the condition is uniformity. The tell is the first figure: over its first three revolutions the adaptive run is nearly six times better, which is exactly the evidence a short test would produce. The property the essay is about is a property of a fixed step, so it is worth drawing at a second eccentricity and a second range of step sizes.

Where this ladder goes next #

This rung has taken a property that a previous rung established and found its actual precondition: the bounded energy error of a symplectic method belongs to a fixed step, and every natural improvement to the step destroys it.

The rung above is the splitting itself. Separating a Hamiltonian into an exactly solvable Kepler part and a small perturbation is what makes a twenty-step-per-orbit integration of the solar system accurate at all, and the choice of how to split it — and of the order in which the pieces are applied — is where the modern work is.

Beside it lies the question of what “accurate” should mean for a chaotic system. If two integrations of the same initial condition diverge exponentially, then no integration is following the true trajectory, and what is being claimed instead is that the computed trajectory is the exact trajectory of a nearby initial condition — a shadowing statement, which is a different and much weaker claim than the one people usually think they are making.

And below it, the habit: check what a guarantee is a guarantee about. The energy band belongs to the step being constant, not to the step being small, and the difference between those two readings is the difference between an integration that means something after five billion years and one that does not.