N-Body Gravity Simulator

Every pair of masses pulls on each other. The pull grows with both masses and falls with the square of the distance between them. That, plus one measured constant, is all the physics a gravity simulation contains. Everything else is bookkeeping.

The bookkeeping is where simulations go wrong. A picture of planets circling a star looks convincing whether the arithmetic behind it is careful or careless, because the eye forgives a great deal. Orbits that slowly swell, a moon that quietly gains speed, a cluster that heats itself up from nothing: all of these look like physics. They are errors, and the only way to tell is to check the quantities that gravity is not allowed to change.

One law, one constant

The constant is G. The value NIST publishes from the CODATA adjustment is 6.67430 x 10^-11 in units of cubic meters per kilogram per second squared, with an uncertainty in the last two digits. What an orbit feels, though, is G times a mass. NASA’s fact sheets for the Sun, the Earth and the Moon list that product on its own line beside the mass, and the simulator multiplies the two published figures together to get it.

Before the law there were the rules it explains. Kepler found that each planet moves on an ellipse with the Sun at one focus, and that the line from Sun to planet sweeps equal areas in equal times. In 1619 he published the third: a precise relation between a planet’s distance from the Sun and the time it takes to go around. NASA’s account notes that Kepler knew nothing of gravity. Newton supplied the force behind the pattern, and his version of the third law lets astronomers weigh any two bodies from their separation and their period.

For two bodies that is the end of the story. The ellipse is the exact answer, forever, and its period follows from the size of the orbit and the summed mass alone. A simulation of two bodies is therefore a test with a known answer. Run a circular orbit for one predicted period and the bodies should be back where they started. The size of the miss is a direct measurement of the method.

Three bodies and the end of formulas

Add a third body and the formulas stop. Šuvakov and Dmitrašinović put the reason in one line: three bodies have 18 degrees of freedom but only 10 integrals of motion, as Bruns showed late in the nineteenth century. Those ten are the energy, three components of momentum, three of angular momentum, and three that describe the steady drift of the center of mass. So the count is 10 of 18, which is not enough to pin the motion down, and the problem cannot be solved in the sense that the two-body problem can.

Henri Poincaré found out how badly it cannot. King Oscar II of Sweden and Norway set up a prize competition in 1887, and Poincaré won it with a memoir on the three-body problem. As the memoir went to press an editor found an error. The MacTutor biography records that the memoir gave the first mathematical description of chaotic motion, and that the error is now regarded as marking the birth of chaos theory.

Chaos does not mean that no orderly solutions exist. It means they are precious, which is the word Poincaré himself used for periodic orbits: the only opening into a place thought to be inaccessible. Lagrange had already found one family in 1772: three bodies at the corners of an equilateral triangle that turns rigidly. A small body at such a corner sits at what is now called the L4 or L5 point of the two large ones.

Those triangular points are stable only when one large mass heavily outweighs the other. NASA gives the condition as a mass ratio above 24.96, and the exact figure, 25 plus three times the square root of 69, all halved, has its own entry in the encyclopedia of integer sequences. The Sun and Jupiter pass easily. That is why asteroids called Trojans gather ahead of and behind Jupiter, and why NASA’s Lucy is described as the first mission to explore them. The three points on the line through the two masses are unstable whatever the ratio.

The strangest orderly solution is recent. In 1993 Moore found three equal masses chasing each other around a single figure-eight. In 2000 Alain Chenciner and Richard Montgomery showed by a variational argument that such an orbit exists. Their paper prints starting positions and velocities computed by Carles Simó to eight decimal places, along with the period, 6.32591398 in units where the masses and G are one. Simó’s computations also indicated that the orbit is stable. A numerical search published in 2013 reported 15 further sets of starting conditions, 13 of them distinct orbits.

What a time step does to an orbit

With no formula, a computer does the only thing left. It freezes the positions, adds up the pull on every body from every other, moves everything forward by a small step of time, and repeats. The work in each pass is one calculation per pair, and N bodies make N (N - 1) / 2 pairs, so doubling the bodies quadruples the cost. Two hundred bodies already mean 19,900 pairs per pass.

How the step is taken matters more than its size. The leapfrog method, also called Störmer-Verlet, nudges the velocities by half a step, moves the positions a whole step, recomputes the forces, and nudges the velocities again. It is only second order: halve the step and the error falls by four. The classical fourth-order Runge-Kutta method samples the forces four times per step, and halving its step cuts the error by sixteen.

On that count Runge-Kutta should always win. It does not, because leapfrog has a property that accuracy alone does not buy. It is symplectic, which means each step is, very nearly, an exact move of some slightly different mechanical system. Hairer’s lecture notes prove the consequence. A symplectic method of a given order keeps its energy error at the size set by the step, over extremely long times, provided the computed solution stays in a bounded region. Leapfrog is also time reversible: run it backward and it retraces its path.

The difference shows in a ledger. Take two bodies on an ellipse of eccentricity 0.5 and give both methods 200 steps per orbit for a thousand orbits. Leapfrog’s energy error peaks at 0.27 percent in the first hundred orbits and at the same 0.27 percent in the last hundred. Runge-Kutta begins nearly thirty times better, near 0.009 percent, and ends ten times worse than it began, still falling. Its orbit is slowly shrinking.

There is a way to have both. Run leapfrog three times in a row with step lengths in a particular ratio, the middle one negative, and the leading errors cancel. Hairer calls it the triple jump and credits it to several authors, Yoshida among them. The result is fourth order, still symplectic, and costs three force calculations where Runge-Kutta needs four.

A worked example: when the ledger says no

The promise of a bounded error has a condition, and a simulator should show the case where the condition fails. In the tool, choose the setup called Chaotic trio: three masses of 400 with a time step of 0.5 and a softening length of 5. Press Compare with the span at 20,000 steps.

All three methods fail. Leapfrog’s largest energy error is about 1.7, meaning 170 percent of the starting energy. The fourth-order symplectic method reaches 3 and Runge-Kutta nearly 15. Nothing on the screen looks wrong: three balls swing around each other. But early in the run two of them pass very close, the fixed step is far too coarse for that passage, and the energy jumps. Symplectic or not, the run after that jump describes a different system.

Now type 0.05 into the time step and compare again. Leapfrog’s worst error falls from 1.7 to 0.025, and the two fourth-order methods land below 0.0002. The step now resolves the close pass, and the guarantee is back. No integrator is trustworthy in general. A step is trustworthy for the closest approach the run actually contains.

Softening is the other way out, and it has a price. It replaces each point mass with a small fuzzy ball, so the force stops growing once two bodies are closer than the softening length. Dehnen and Read describe both sides: close encounters are suppressed, which makes the equations far easier to integrate, and gravity is systematically too weak at short range. Any orbit tighter than the softening length is wrong by construction.

What a sandbox cannot tell you

The model is Newton’s and nothing more. There is no relativity. Bodies are points, so there are no tides, no spin, and no collisions except an optional merge that conserves mass and momentum and throws the lost energy away.

Nor is a solar system built from a fact sheet an ephemeris. The tool starts each planet at its closest distance to the Sun on a flat ellipse made from NASA’s published perihelion and aphelion. That gets periods about right: within 1.3 percent of the sheet’s own period column, with Neptune the worst. It says nothing about where any planet is on a given date. For that there is JPL’s Horizons service, which exists to produce highly accurate ephemerides.

What the sandbox can do is make the bookkeeping visible. Momentum and the center of mass should hold to the last digit the machine carries, and with every method here they do. Angular momentum holds that well for the two symplectic methods and visibly less well for Runge-Kutta. Energy is the one to watch. A flat line means the step suits the motion. A staircase means a close pass was skipped over. A steady slope means the method itself is leaking.