Astronomy & Celestial Mechanics
the sky, computed
Series

Numerical integration — the series

3 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

All series