Wrong about where, and right about how much

The two-body problem is solved, so nobody integrates it. That is exactly what makes it the right thing to integrate: the answer is a conic and the energy is a constant, so every departure a numerical method produces is a property of the method rather than of the problem. Hand the same orbit to three integrators, run it for a few hundred revolutions, and what comes back is a portrait of what each one does when nobody is watching.

The result is not the one the accuracy of each method would suggest.

Three methods, one orbit #

Euler’s method is the definition of a derivative used as an instruction: take the current velocity, multiply by the step, add it to the position; take the current acceleration, multiply by the step, add it to the velocity. First order, and its failure on an orbit is not subtle. Each step’s position update uses the velocity from the beginning of the step, which on a curving path always points slightly outward from where the body will actually be. The error has a direction, the same direction every step, and the orbit spirals out. Enough steps and the body leaves.

Runge–Kutta 4 evaluates the acceleration four times per step at cleverly chosen intermediate points and combines them so that the leading four terms of the Taylor expansion are matched. Fourth order: halve the step and the error falls by sixteen. It is the default numerical integrator of most of science for good reasons.

Leapfrog — kick, drift, kick — is the construction Newton drew to derive the area law turned into arithmetic. It updates the velocity by half a step, then the position by a full step using that half-updated velocity, then the velocity by the remaining half using the acceleration at the new position. Second order. Two force evaluations per step against Runge–Kutta’s four, and by every ordinary measure the worse method.

The property that is not accuracy #

What the energy figure shows is that “how large is the error” and “what does the error do over time” are different questions with different answers.

Runge–Kutta’s error at any given moment is small. Its error at the end of a long run is not, because it accumulates in one direction. Every step deposits a tiny amount of energy in the same sign, and a tiny amount deposited a hundred million times is not tiny. Nothing in the construction of Runge–Kutta forbids this: the method matches Taylor terms, and a Taylor expansion has no opinion about conserved quantities.

Leapfrog’s error at any given moment is larger, by more than two orders of magnitude at the same step size. Its error at the end of the run is the same size as its error at the beginning. The curve wobbles up and down within a band whose width is set by the step size, and the band does not widen.

The reason is structural rather than fortunate. Leapfrog is the exact solution of a different problem. There exists a Hamiltonian $\tilde H$, differing from the true one by terms of order $h^2$, whose exact trajectories are precisely the points leapfrog produces. So the sequence of states it generates conserves $\tilde H$ exactly — not approximately, exactly — and since $\tilde H$ differs from $H$ by a bounded amount, $H$ is bounded along the numerical trajectory too. The method has not been made accurate. It has been made to solve the wrong problem consistently, and consistency is what a conservation law is.

That is the content of the word symplectic. A Hamiltonian flow preserves a particular structure on phase space — areas, in the two-dimensional case — and a symplectic integrator is one whose single step is itself a map preserving that structure. Composing structure-preserving maps gives a structure-preserving map, however many of them there are, which is why the guarantee does not decay with the length of the run.

And neither of them knows where the planet is #

Here is the half of the story that a summary of the form “symplectic integrators conserve energy” throws away.

A symplectic integrator’s semi-major axis stays where it should. Its phase does not. The body goes round very slightly too fast or too slow, and after enough revolutions it is on the far side of the orbit from where it should be — at which point the position error saturates at about the size of the orbit and stops growing, having exhausted the available places to be wrong.

So a five-billion-year integration of the solar system produces two kinds of statement, from one run, and they have to be sorted:

• About the shape of the system — do the planets stay on orbits of roughly these sizes and eccentricities, does anything cross anything else, does anything leave — the run is evidence, because those questions are about the conserved structure and the structure is what the method preserves — the same structure that makes six numbers a complete description of an orbit in the first place.
• About where anything is — the run is worthless, because the phase error is unbounded and, as it happens, the physics makes it worse than the numerics do.

What was actually measured #

Every number in the figures above came from arithmetic, and it is worth being exact about which parts of this subject are observational and which are not.

The integrator comparison is not. It is a controlled experiment on a problem with a known answer, run at a stated step size, and its findings are facts about the algorithms.

What is observational is the initial conditions, and they are the reason the whole exercise has a horizon. A long-term solar-system integration starts from a planetary ephemeris: a fit to radar ranging to the inner planets, spacecraft tracking at Mars, Jupiter and Saturn, and two centuries of optical astrometry. The Earth–Mars distance is known to a few metres by ranging, so the positions that begin the integration carry errors of that order — the best-determined initial conditions in dynamics, and not remotely good enough, because the errors do not stay that size.

The second observational input is the mass of everything, and it enters as $GM$ per body rather than as a mass, for the reason that the gravitational constant is the worst-measured constant in physics. The asteroids are the practical limit: their combined $GM$ perturbs Mars measurably, the individual masses of most of them are unknown, and a long integration has to model them as a smoothed ring or accept an error it cannot bound.

The measurement that settles whether any of this works is a backwards one. An integration run backwards from modern initial conditions predicts the circumstances of eclipses recorded in Babylon, in China and in the Arab world — dates, and more usefully the geographic tracks of totality. Those tracks are reproduced, once the Earth’s own rotation is allowed to slow, across 2,700 years. That is about $10^6$ orbits of the Moon, which is a real test of a numerical scheme and a small fraction of the age of the solar system.

What the integrator is actually solving #

The bounded energy error looks like a happy accident and it is a theorem, and the statement of it explains the behaviour exactly.

A symplectic integrator applied to a Hamiltonian system does not approximately solve that system’s equations. It exactly solves the equations of a slightly different Hamiltonian — a shadow Hamiltonian, differing from the true one by terms of the order of the step size raised to the method’s order.

That single fact accounts for everything observed. The shadow system’s energy is conserved exactly, because it is a Hamiltonian system being solved exactly; and the true energy differs from the shadow energy by a fixed small amount that depends on the position in the orbit but does not accumulate. So the true energy oscillates within a band and never drifts.

A non-symplectic method has no such shadow. Its error at each step is a genuine error rather than a faithful solution of a nearby problem, and successive errors accumulate in whatever direction the truncation happens to favour — which for the energy is generally one direction.

The comparison is therefore not between a more accurate method and a less accurate one. It is between a method that answers a slightly wrong question correctly and one that answers the right question approximately, and over a billion orbits the first is what a long integration needs.

The step that must not change #

The theorem has a condition attached, and it is the practical constraint on every use of these methods: the shadow Hamiltonian depends on the step size, so changing the step changes which problem is being solved exactly.

An adaptive integrator that shortens its step near periapsis and lengthens it afterwards is therefore hopping between shadow systems, and the energy conservation is destroyed — the error resumes accumulating, and the method loses the property it was chosen for.

That is awkward, because the situations where a fixed step is expensive are the ones where a variable step is most wanted: a close encounter demands a step hundreds of times shorter than the rest of the orbit, and paying that everywhere is unaffordable.

The standard resolution is a hybrid. Integrate with the fixed-step symplectic method while all separations are large, and switch to a conventional adaptive method for the duration of any close encounter, switching back afterwards. The switch introduces an energy error each time, but encounters are rare, so the errors are occasional discrete jumps rather than a continuous drift.

The design is therefore an admission: the good long-term behaviour is available only while nothing interesting is happening, and the interesting events have to be paid for one at a time.

It is also why the published long-term integrations of the solar system state their step size and their encounter handling in the abstract: those two choices decide what the result is worth, and neither is a detail of the implementation.

Where the model stops #

The comparison above is on the two-body problem, and the two-body problem has no chaos in it. Everything the figures show about phase error is a statement about accumulated arithmetic. In the real solar system there is a second source of the same symptom and it is far larger.

There are also failure modes that belong to the methods rather than to the problem.

Close approaches. A fixed-step symplectic integrator handles a nearly-circular orbit beautifully and a close encounter appallingly, because the timescale collapses and the step no longer resolves it. Adapting the step destroys the symplectic property — the composition is no longer of one fixed map — so the standard remedies are elaborate: a symplectic corrector, a change of time variable that is itself canonical, or a hybrid scheme that switches to a high-order non-symplectic method for the duration of the encounter and switches back. Each of these is a compromise and the literature about which compromise is a substantial one.

Round-off. Over $10^{11}$ steps, double-precision rounding is not negligible; it is a random walk in the energy of amplitude $\sqrt{N}$ times the machine epsilon, which for $10^{11}$ steps is about $10^{-11}$ — comparable with the truncation error the scheme was designed around. Long integrations use compensated summation for this reason, which is a bookkeeping trick rather than a numerical method and recovers most of it.

The model itself. An integration reproduces whatever forces it is given. General relativity contributes Mercury’s 43 arcseconds a century, which over five billion years is a complete extra turn of the orbit twenty times over; leaving it out changes the answer about Mercury’s stability entirely, because it shifts the precession rate out of a secular resonance with Jupiter. So do the tides, the solar mass loss, and the Moon receding from the Earth — each a small term, each integrated for a very long time. Each is a small term, and the whole subject is about what small terms do when integrated for long enough.

What the picture cannot show #

Any real system. Every figure here is one orbit or three point masses. The solar system is a nine-body problem at least, with several hundred thousand small bodies whose combined effect is not negligible, and nothing about a two-body test says how a scheme behaves on it.

The cost. The energy figure holds the step size fixed and compares methods, which is the wrong comparison for a practitioner: Runge–Kutta uses four force evaluations per step to leapfrog’s two, so at equal cost the comparison is leapfrog at half the step, which halves its error band. The right axis for that question is wall-clock time, and it is not on any of these plots.

Whether the answer is right. These figures compare integrators against an exact solution, which exists only because the test problem was chosen for it. On the problem anybody cares about there is no exact solution to compare against, and the working substitute — run it again at half the step and see whether the answer moves — is a test of convergence rather than of correctness. That is the same distinction the zero-velocity curve draws between what an integration can establish and what a conserved quantity can: one is evidence and the other is proof, and the restricted three-body problem is the last place in this subject where proof is available.

The result that changed what the question was #

Before 1988 the question “is the solar system stable” was understood as having a yes-or-no answer that a long enough integration would deliver.

Jacques Laskar’s integration of the secular equations, and then the direct integrations that followed, changed the question rather than answering it. What came out was a Lyapunov time for the inner planets of about five million years — so an error of a centimetre in Mercury’s position grows to the size of its orbit in a hundred million years, and no initial conditions that could ever be measured will determine where Mercury is at the age of the Sun.

What survived was a different kind of statement. Run the integration many times from initial conditions scattered within the observational errors, and count: in about one per cent of the runs Mercury’s eccentricity grows enough over five billion years for it to cross Venus. That is the modern answer, and it is a probability over an ensemble rather than a prediction. It could not have been reached without integrators whose energy error is bounded, because a run whose semi-major axes are slowly drifting for numerical reasons cannot be distinguished from one whose semi-major axes are drifting for physical ones, and the whole result is about which of those is happening.

The separation between the two errors is the essay’s subject, and both of them are worth reading at a second eccentricity.

Where the ladder goes next #

Later rungs on this anchor: the Wisdom–Holman mapping, which splits the Hamiltonian into a Kepler part solved exactly and a perturbation applied as a kick, and buys a factor of a thousand in speed for solar-system work. Symplectic correctors and the higher-order compositions that reach fourth and sixth order while staying symplectic. Regularisation of close encounters, and the Kustaanheimo–Stiefel transformation that removes the collision singularity outright. Backward error analysis as a subject, of which the shadow Hamiltonian above is one instance. And the tree and multipole methods that make an $N$-body problem cost less than $N^2$, which is what turns a planetary integrator into a galactic one.

Leapfrog is older than the analysis that justifies it: it appears in Newton’s Principia as the geometrical construction by which he derives the area law, a century and a half before Hamilton and three before anyone said the word symplectic. Newton drew it because it was the obvious way to turn a continuous force into a sequence of impulses. It turns out to be the right way, for a reason he had no vocabulary to state.