Astronomy & Celestial Mechanics
the sky, computed
Series

Numerical integration — the series

5 essays on one idea, from the one that introduces it to the one that assumes the rest.
  1. Energy error over 240 revolutions, at one step size. The relative error in total energy against revolution number, for three integrators run on the same Kepler orbit at e = 0.5 with the same step of 200 per revolution. The exact energy is a constant, so every curve here is the method rather than the problem. Euler climbs steadily: its energy at the end is 106.5% wrong, and the orbit it draws has spiralled outwards. Runge–Kutta 4 begins 4.7e+2 times more accurate than leapfrog and ends at 2.22e-4, having grown by a factor of 10 across the run: the error is SECULAR. Leapfrog oscillates inside a band and stays there — worst error 2.62e-3, and the second half of the run is no worse than the first, which is measured here rather than claimed. That is the property that decides whether a five-billion-year integration means anything, and it is not accuracy: a symplectic method is the exact solution of a Hamiltonian a step-size away from the intended one, so its energy cannot wander, while a more accurate non-symplectic method has no such constraint and eventually wanders further.

    Wrong about where, and right about how much

    Runge–Kutta is the more accurate method and loses energy steadily; leapfrog is cruder and its energy error never leaves a band. Over five billion years only one of those properties survives — and neither method knows where the planet is.

    part 1 · gravitation
  2. The adaptive step starts 5.7× more accurate and ends 7.5× worse. The envelope of the relative energy error against elapsed revolutions — the worst error within each plotted interval rather than the error at one sample in it — for the same second-order symplectic integrator run two ways on a Kepler orbit at e = 0.5. The flat band is a fixed step of one 200th of a period: its error oscillates once per revolution and the envelope does not grow, because a symplectic method at constant h is the exact solution of a nearby Hamiltonian and is conserving that one. The rising curve is the same method with the step refined where the orbit is fast — h ∝ r^3/2, the free-fall time, varying by a factor of 5 around the orbit — which is the first thing anybody reaches for at a close encounter and which is strictly more accurate step for step: over its first three revolutions it stays a factor of 5.7 below the fixed run. By 1500 revolutions it is a factor of 7.5 above it and still climbing. Changing the step changes which Hamiltonian is being conserved, the errors from successive steps stop cancelling, and what is left is a random walk with no bound at all. The practical consequence is that a solar-system integration cannot adapt its step: it either keeps a step short enough for the closest encounter it will ever meet, or it detects the encounter and hands that piece of the trajectory to an entirely different, non-symplectic method for the duration — which is what every long-term integration of the planets actually does.

    A step that must not be adapted

    A symplectic integrator's bounded energy error is a property of a fixed step. It is conserving a Hamiltonian a step-size away from the intended one, and changing the step changes which Hamiltonian — so refining the step at a close encounter, which is the first thing anybody does, destroys the only property the method was chosen for.

    part 2 · gravitation
  3. A step size below which a smaller step is worse. The error left in a long integration against the step size, for methods of three different orders, with both contributions drawn. The falling lines are truncation error, whose slope on these axes is exactly the order of the method. The rising line is round-off, identical for all three because it is a property of the arithmetic and not of the algorithm: every operation loses a few bits, the losses are independent, and they accumulate as the square root of the number of steps — which is why its slope is −1/2 and why it rises as the step shrinks. Each method's total has a minimum, at a step of 1.0e-6, 5.0e-6, 9.7e-4 for orders 1, 2, 4. Below that minimum every halving of the step costs time and makes the answer worse. That is the practical reason a solar-system integration is not run at an arbitrarily fine step, and it is a reason with nothing to do with computer time.

    An error that grows like a random walk

    A long integration accumulates two errors with opposite habits. One falls when the step is made smaller and grows in proportion to the time; the other grows when the step is made smaller and accumulates as a square root. Which of the two dominates decides whether a billion-year integration means anything.

    part 3 · gravitation
  4. Sixty futures that began a trillionth apart. The momentum of the standard map at K = 3, a kicked rotor whose phase space is almost entirely chaotic, for 60 members started 10⁻¹² apart in angle, against iteration. For 41 iterations they are one trajectory; then the separation, growing by a factor of about 2.0 each step, reaches the size of the system and they become independent random walks. Three of the 60 cross the threshold at |p| = 90 within 400 iterations (dark), the stand-in for an orbit that goes unstable. Which members cross is decided by differences far below anything that could be measured; how many cross is a property of the dynamics, and it is the only thing the ensemble can be asked for.

    A probability for a system with one future

    The Solar System will do one thing, and an integration of it is wrong about which thing within a hundred million years. What long integrations report instead is a probability — one per cent that Mercury's orbit goes unstable before the Sun dies — and a probability of a single event is a strange object. It is a count over an ensemble of futures that differ by less than anything measurable, which makes sense only because chaos makes those futures independent; and how many of them are needed is a binomial question with an unforgiving answer for rare events.

    part 4 · gravitation
  5. A rare crossing reached in stages. Multilevel splitting on a kicked rotor (the standard map at K = 10), aimed at the event |p| > 450 within 400 iterations. 16 members are run until each reaches the first level or runs out of time; those that reach it are copied back up to 16, each copy restarted from its parent's crossing point with its angle nudged by a thousandth of a radian, and run on toward the next. The levels (dashed) are not fixed values of |p| but curves that rise as time runs out, so a member counts as closer to the event only if it has the time left to finish. 4 levels; the fractions passing them are 0.13, 0.06, 0.63, 0.88, and their product, 0.0043, is the estimate. Each stage is a common event; the rare one is never waited for.

    A rare event reached by cloning the near misses

    Counting how often an ensemble of integrations goes unstable works for an event of one in a hundred and becomes hopeless for one in a million, because a ten per cent answer needs a hundred events and so a hundred million futures. Splitting gets there by a different route: wait for the members that come closest, copy them, and run the copies on. Each stage measures a common event and the rare one is a product of common ones — provided the copies forget their parents and the stages know how much time is left.

    part 5 · gravitation

All series