How to Lose a Planet by Choosing the Wrong Integration Method

One Equation, Applied Over and Over

Gravity between two bodies is not a complicated law. One mass pulls the other straight toward itself with a strength that falls off as the square of the distance. Write the equation down and a computer can step it forward in time as easily as it adds two numbers, because stepping it forward in time is, quite literally, all that "integrating" an orbit means: take the current position and velocity, apply a small correction, and repeat, tens of thousands of times, until what comes out the other end looks like a planet going around a star.

The catch is that "apply a small correction" is not one algorithm. It is a family of them, and family members disagree, sometimes wildly, about how much energy the orbit is allowed to keep by the end. Choose the wrong member of the family and the planet you just simulated does something no real planet does: it gains energy for free, or loses it for free, for reasons that live entirely in the arithmetic and nowhere in the physics. Watch that happen once, on your own laptop, on an orbit you set up yourself, and it changes how you read every simulated trajectory you did not personally check.

Real orbital-dynamics research takes this seriously for reasons well beyond a school club's curiosity. The solar system itself is chaotic on timescales of tens of millions of years, meaning two integrations that start an immeasurably tiny distance apart, including two runs of the same code with different roundoff, diverge exponentially rather than staying close [15]; serious long-term integrations of planetary orbits are correspondingly run at extremely fine resolution and checked against independent codes for exactly this reason [6]. Nothing in our own two-body problem is chaotic. The failures below are not sensitivity to initial conditions; they are a wrong answer, repeatably, from a fixed method at a fixed resolution, which is in some ways a plainer thing to demonstrate.

This is a methods study, not a discovery about gravity. The two-body problem is one of the few gravity problems with a closed-form answer, and we chose it for exactly that reason: every number below is checked against the exact Kepler ellipse, not against another simulation. Nothing here was observed. The computation is the whole experiment, and every figure, every table entry and every claim in the abstract traces back to a line printed by the program listed under "Analysis code" above.

Three Ways to Take a Step

Work in units where $GM=1$ and the semi-major axis $a=1$, which makes the orbital period exactly $T = 2\pi$ and turns every quantity below into a clean, dimensionless number. The state of the planet is its position and velocity, $(x, y, v_x, v_y)$, and the acceleration is always $\mathbf{a} = -\mathbf{r}/r^3$. Three rules turn one state at time $t$ into the next state at time $t+\Delta t$:

Forward Euler. The oldest and simplest rule there is: $x_{n+1} = x_n + v_n\Delta t$, $v_{n+1} = v_n + a(x_n)\Delta t$. One force evaluation per step. First order, meaning its error per fixed time interval should shrink in proportion to $\Delta t$ itself.

Velocity Verlet. Split each step into a half-kick, a full drift, and a second half-kick: $v_{n+1/2} = v_n + \tfrac12 a(x_n)\Delta t$, $x_{n+1} = x_n + v_{n+1/2}\Delta t$, $v_{n+1} = v_{n+1/2} + \tfrac12 a(x_{n+1})\Delta t$. Also one force evaluation per step, since the acceleration computed at the end of one step is reused as the start of the next. Second order, and symplectic, meaning the discrete map it defines preserves phase-space volume exactly [7], a geometric property Euler and RK4 both lack. Higher-order symplectic methods exist too, built by composing several Verlet-like substeps with carefully chosen weights [4]; we did not test one, and §7 explains why a fourth-order symplectic method would have been a fairer opponent for RK4 than plain second-order Verlet.

Fourth-order Runge-Kutta (RK4). The workhorse of general-purpose numerical integration: four intermediate force evaluations per step, combined with weights $1,2,2,1$, give fourth-order accuracy, meaning its error per fixed time interval should shrink as the fourth power of $\Delta t$. Four times the cost of Euler or Verlet, step for step, and, as the theory section below shows plainly, not symplectic despite the higher order, a point the broader symplectic-integration literature makes at length [8] and one that genuinely surprises people meeting it for the first time.

Cost is charged the honest way, in force evaluations, since that is what actually eats compute time in a real simulation: one per step for Euler and Verlet, four per step for RK4. That accounting matters later, in §9, once the question stops being "which method is more accurate at the same timestep" and becomes "which method is more accurate for the same money."

Checking the Machine Before Trusting It

Two validations ran before any orbit was trusted.

A: does the Kepler-equation solver actually solve the equation? Every comparison below depends on an exact reference orbit, obtained by solving Kepler's equation $M = E - e\sin E$ for the eccentric anomaly $E$ by Newton's method [12]. Checked over 740 combinations of mean anomaly and eccentricity spanning everything used in this study, the worst residual was $8.882\times10^{-16}$, at machine precision, needing at most 9 iterations to converge. Club value: worst residual $8.882\times10^{-16}$. Accepted value: exactly 0. Difference: $8.882\times10^{-16}$, the floor of double-precision arithmetic, not a real disagreement.

B: does the reference orbit actually conserve energy and angular momentum? Sampling the closed-form solution at 4,000 points across 7 orbits at $e=0.6$ and recomputing energy and angular momentum independently from the sampled position and velocity, the worst relative deviation from the theoretical constants ($E_0=-GM/2a=-0.5$ exactly, $L_0=\sqrt{GMa(1-e^2)}=0.8$ exactly) was $4.441\times10^{-15}$ for energy and $5.551\times10^{-16}$ for angular momentum, both at the double-precision floor.

Both checks passed before a single stepping rule was run. The rest of this article compares three approximate methods against this verified exact one, never against each other directly.

The Headline Number

No fitting in this section. Just what the three integrators actually did, in plain arithmetic, at the one configuration used everywhere in this article: eccentricity 0.6, one hundred steps per orbit, $\Delta t = T/100 \approx 0.0628$.

Start energy: $E_0 = -0.500000000000000$, exact by construction. After a single orbit, Euler's energy is $E_0 + 1.166583 \times |E_0| = -0.5 + 0.583 = +0.083$. Positive. The two-body problem with positive energy is not an ellipse; it is unbound. Euler did this inside one orbit, at a resolution, one hundred steps to go around once, that would strike most people as more than enough. It does not recover: 300 orbits later the energy is $+0.0828$, a number within a rounding error of where it landed after orbit one. Once Euler pumps that much energy in, there is no elliptical orbit left to drift away from.

Verlet, same 300 orbits: energy after orbit 1 is $-0.4865$; after orbit 300 it is $-0.4855$. Those two numbers are close not because Verlet stopped moving, but because it keeps returning to roughly the same place. In between, every single orbit, it swings from close to $0\%$ relative error up to about $2.99\%$ and back (Figure 2, lower panel; see §5). RK4, same 300 orbits: $-0.500105$ after orbit 1, creeping to $-0.532241$ by orbit 300, a small, steady contraction, more bound with each orbit rather than less.

"Permanent" means the trace never falls back under half the threshold again for the rest of the run; "transient" means it does. Verlet crosses 1% once, briefly, two hundredths of an orbit in, then falls back and spends the remaining 299.98 orbits oscillating inside a band it never permanently leaves.

How Bad Does It Get

Figure 1 is the honest picture, not a flattering one. Euler's trace is not merely off the ellipse, it dwarfs it: by $x=-9.45$ the true orbit's entire diameter, 2 units, would fit inside the gap between Euler's excursion and the sun nearly five times over. Verlet and RK4, by contrast, are close enough to the dashed exact curve that separating them at this scale takes the next figure, not this one.

Two shapes, not one. Verlet's error is a band with a fixed width, set once, near the start, and never revisited. RK4's error is a line with a slope, small at first, unmistakable by orbit 150. The bound on Verlet's band is not a coincidence of this one run: backward error analysis shows a symplectic integrator's numerical trajectory exactly conserves a nearby "shadow" Hamiltonian, close to the true one by an amount set by the timestep, so the true energy can only ever wander within a fixed distance of its start rather than drift away from it [10], a result laid out in full in the standard reference on the subject [11]. The ordinary-least-squares slope of RK4's relative energy error against orbit number is $-2.146\times10^{-4}$ per orbit with a standard error of $9.79\times10^{-9}$, a slope 21,915 times its own standard error: about as far from "consistent with zero" as a number in this article gets. Verlet's equivalent slope is $-1.26\times10^{-8}$ per orbit against a standard error of $3.36\times10^{-7}$, zero to zero decimal places that matter. Euler's slope, for what it is worth after the orbit is already gone, sits at 2.6 standard errors from zero, a small residual wobble inside an already-ruined trajectory, interesting only for showing that the chaotic wandering after collapse is itself not perfectly static.

Validation C: angular momentum. Symplectic theory makes a second, sharper prediction for velocity Verlet specifically. Over the same 300 orbits, Verlet's angular momentum drifted by at most $2.387\times10^{-14}$ relative to its starting value, 107 machine epsilons, against a naive roundoff random walk estimate of $3.85\times10^{-14}$ over 30,000 steps: a ratio of 0.62, consistent with pure floating-point noise and nothing else. Forward Euler's angular momentum, over the identical run, drifted $1.48\times10^{16}$ times further; RK4's drifted $3.06\times10^{11}$ times further. Club value: Verlet max $|\Delta L/L_0| = 2.387\times 10^{-14}$. Accepted value: exactly zero, up to roundoff, as the theory in the margin note requires. Difference: consistent with roundoff, not a real disagreement.

Does the Error Actually Shrink the Way Theory Says

A stepping rule's "order" is a specific, falsifiable claim: halve the timestep and the global error after a fixed amount of simulated time should shrink by a fixed factor, $2^1$ for a first-order method, $2^2$ for second order, $2^4$ for fourth. We tested this directly, at $e=0.6$, by integrating exactly one orbit at seven different step counts per method and measuring the Euclidean position error against the exact Kepler solution at the end, then fitting a line through the results on log-log axes.

Each of the three methods needs its own range of step counts before this scaling shows up cleanly, and that range differs by two orders of magnitude between methods. Verlet and RK4 settle into their asymptotic slope, 2 and 4, by a few hundred steps per orbit. Euler does not: at a few thousand steps per orbit its local fitted slope is still only about 0.8, and it takes tens of thousands of steps per orbit, several thousand times finer than the 100-step configuration used everywhere else in this article, before the fitted slope on a shrinking window of step counts settles at 1.000. Butcher's textbook order conditions [16] describe the limit $\Delta t \to 0$; they say nothing about how close to that limit a given problem needs to get before the limit is a useful description, and for Euler on this orbit, the honest answer is: much closer than anyone would run it in practice.

Working Notes From the Club Table

Is This Just One Lucky (or Unlucky) Orbit?

Every number so far comes from one eccentricity, 0.6, chosen for no better reason than that it looks like a real elliptical orbit and stresses all three methods enough to be interesting. A skeptical reader is entitled to ask whether the measured convergence orders in §7 are an artifact of that one choice. We do not think they should be taken on faith, so we did not take them on faith ourselves: the same one-orbit convergence fit was repeated 40 times, at eccentricities drawn independently and uniformly from $[0.10, 0.60]$, seeded once from the master seed via numpy's SeedSequence so the draw is exactly reproducible.

The Monte Carlo means agree with §7's single-orbit numbers to within two percent for every method, so the single eccentricity was not a fluke in either direction. What the larger sample buys is a standard error small enough to make the residual disagreement with theory impossible to shrug off: RK4's mean fitted order, 4.118, sits 23.1 of its own standard errors above the textbook value of 4. That is not sampling noise. We think it is a genuine, small, higher-order correction term in RK4's local truncation error that has not fully died away even at the finest step counts we used, an effect that would itself shrink with even finer resolution, the same phenomenon §7 already documents for Euler at a much larger scale. We report the 23.1-sigma gap rather than narrow the step-count range until it disappears, because doing the latter would be measuring what we wanted to see rather than what the code produced.

Buying Back Stability

Comparing methods at the same timestep is comparing them at different costs, since RK4 spends four force evaluations per step against Euler and Verlet's one. The fair comparison holds cost fixed instead: RK4 at 60 steps per orbit (240 force evaluations per orbit) against Euler and Verlet at 240 steps per orbit each (also 240 force evaluations per orbit, four times finer a timestep, for the same total spend).

Giving forward Euler four times the timestep resolution of RK4, for the same money, does not buy it a longer life at $e=0.6$ or $e=0.2$. It buys it a slightly less immediate death: the matched-timestep table (all three at 60 steps per orbit, $e=0.6$) has Euler crossing 1% at orbit 0.0002; the matched-cost table has it crossing at orbit 0.001, five times longer, measured in a unit where "five times longer than instantaneous" is still, for every practical purpose, instantaneous. This is the difference between an accuracy problem and a stability problem. A bigger budget fixes the first kind. It does not fix the second, and at this eccentricity and this order of magnitude of timestep, Euler's failure is squarely the second kind.

The Strongest Objection We Could Make

The obvious response to everything above is: nobody serious runs forward Euler on an eccentric Kepler orbit, so what has actually been demonstrated beyond a straw man? This objection has real force and deserves a real answer rather than a dismissal.

The answer is not that Euler is secretly fine. It is that "nobody serious uses it" is doing more work in that sentence than it looks like. Plenty of serious N-body and planetary-dynamics codes do not use velocity Verlet either, for reasons that have nothing to do with the failure mode documented here. The Hermite integrator, fourth order and time-symmetric like RK4 but built around continuously estimated higher derivatives rather than substeps, is a standard choice in direct-summation planet-formation codes precisely because it tolerates the close encounters and wide range of dynamical timescales that a fixed-step symplectic method handles badly [14]. It is not symplectic either, and by our own §5 logic it should show the same kind of secular drift RK4 shows, just at a rate its designers judged an acceptable trade for adaptive timestepping and better handling of close, fast passages. Large-scale cosmological N-body codes such as GADGET-2, on the other hand, default to exactly the kick-drift-kick leapfrog scheme used here, precisely for the bounded-error property §5 and §6 demonstrate [9], and modern solar-system integrators such as REBOUND ship a Wisdom-Holman symplectic mapping as their recommended default for exactly this reason [13], following the original demonstration that symplectic mapping methods hold up over solar-system-length integrations where naive methods do not [3]. The honest picture is not "symplectic good, everything else bad." It is "know which failure mode your problem can tolerate, because every fast method accepts one." A close-encounter-dominated planet-formation simulation can tolerate slow secular energy drift far better than it can tolerate a fixed fine timestep everywhere, including the 99% of the simulation where nothing close is happening; a solar-system-length stability study is the opposite case almost exactly. This article picked the second case and is honest that it did.

Where a Different Modelling Choice Would Change the Answer

Three choices made above would each have changed a headline number if made differently. The eccentricity, 0.6 by our own arbitrary pick, sets how fast the failure arrives: §9's table shows RK4 surviving 3.95 orbits at $e=0.6$ and over 1,500 at $e=0.1$ for the identical cost budget, so quoting "RK4 crosses 1% loss around orbit 4" without the eccentricity attached is quoting a number that depends on a choice we have not defended on physical grounds, only made for pedagogical contrast. The step-count range used for the convergence fit, §7's finding that Euler needs tens of thousands of steps per orbit before its slope reaches 1.000, would have reported a wrong order (0.8, not 1.0) had we stopped the fit at the same few-hundred-step range that works fine for Verlet and RK4; a study that fits all three methods over one shared step-count range, a natural and common shortcut, would understate Euler's asymptotic order without ever running a buggy line of code. And the choice to define cost in force evaluations rather than wall-clock seconds sidesteps implementation quality entirely; a poorly vectorised RK4 and a hand-tuned Euler could reverse every ranking in §9's table on wall-clock time alone, which is exactly why we charged cost the way we did and said so, rather than letting a seconds column quietly smuggle in an implementation detail dressed up as a finding about integrators.

Reproducing this

One file, one command. Python 3.12 and numpy; nothing else is imported.

python integrator-energy-drift.py > integrator-energy-drift-output.txt

Our run took 65.0 seconds on numpy 2.4.2 and Python 3.12.3, comfortably inside the ten-minute budget for a club analysis script. Master seed 20260901 governs every random draw in the Monte Carlo section (§8) and nowhere else; every other number in this article is fully deterministic and will reproduce exactly on a rerun. The interactive model runs the same three integrators, at the same default eccentricity and timestep, live in your browser, so you can watch Euler's orbit leave the frame in real time rather than read about it.

References

python integrator-energy-drift.py > integrator-energy-drift-output.txt

1. 1.

   Verlet, L. (1967). Computer "experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules. Physical Review 159, 98–103. doi:10.1103/PhysRev.159.98

2. 2.

   Swope, W. C., Andersen, H. C., Berens, P. H. & Wilson, K. R. (1982). A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: application to small water clusters. Journal of Chemical Physics 76, 637–649. doi:10.1063/1.442716

3. 3.

   Wisdom, J. & Holman, M. (1991). Symplectic maps for the n-body problem. The Astronomical Journal 102, 1528–1538. doi:10.1086/115978

4. 4.

   Yoshida, H. (1990). Construction of higher order symplectic integrators. Physics Letters A 150, 262–268. doi:10.1016/0375-9601(90)90092-3

5. 5.

   Hut, P., Makino, J. & McMillan, S. (1995). Building a better leapfrog. The Astrophysical Journal 443, L93–L96. doi:10.1086/187844

6. 6.

   Quinn, T. R., Tremaine, S. & Duncan, M. (1991). A three million year integration of the earth's orbit. The Astronomical Journal 101, 2287–2305. doi:10.1086/115850

7. 7.

   Channell, P. J. & Scovel, C. (1990). Symplectic integration of Hamiltonian systems. Nonlinearity 3, 231–259. doi:10.1088/0951-7715/3/2/001

8. 8.

   Sanz-Serna, J. M. (1992). Symplectic integrators for Hamiltonian problems: an overview. Acta Numerica 1, 243–286. doi:10.1017/S0962492900002282

9. 9.

   Springel, V. (2005). The cosmological simulation code GADGET-2. Monthly Notices of the Royal Astronomical Society 364, 1105–1134. doi:10.1111/j.1365-2966.2005.09655.x

10. 10.

    Benettin, G. & Giorgilli, A. (1994). On the Hamiltonian interpolation of near-to-the-identity symplectic mappings with application to symplectic integration algorithms. Journal of Statistical Physics 74, 1117–1143. doi:10.1007/BF02188219

11. 11.

    Hairer, E., Lubich, C. & Wanner, G. (2006). Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, 2nd edition. Springer-Verlag, Berlin. ISBN 978-3-540-30666-5.

12. 12.

    Danby, J. M. A. (1987). The solution of Kepler's equation, III. Celestial Mechanics 40, 303–312. doi:10.1007/BF01235847

13. 13.

    Rein, H. & Liu, S.-F. (2012). REBOUND: an open-source multi-purpose N-body code for collisional dynamics. Astronomy & Astrophysics 537, A128. doi:10.1051/0004-6361/201118085

14. 14.

    Kokubo, E., Yoshinaga, K. & Makino, J. (1998). On a time-symmetric Hermite integrator for planetary N-body simulation. Monthly Notices of the Royal Astronomical Society 297, 1067–1072. doi:10.1046/j.1365-8711.1998.01581.x

15. 15.

    Laskar, J. (1989). A numerical experiment on the chaotic behaviour of the Solar System. Nature 338, 237–238. doi:10.1038/338237a0

16. 16.

    Butcher, J. C. (2016). Numerical Methods for Ordinary Differential Equations, 3rd edition. Wiley, Chichester. doi:10.1002/9781119121534