N-Body Gravity SandboxBeta

About N-Body Gravity Sandbox

This sandbox sums Newtonian gravity over every pair of bodies on screen, with no tree or grid approximation (direct summation, so the cost grows with the square of the body count). Positions and velocities are advanced with velocity Verlet in its kick-drift-kick form: give each velocity half a step of acceleration, move every position a full step, recompute forces at the new positions, then apply the second half-kick.

The Energy drift readout compares the current total energy (kinetic plus gravitational potential) with its value when the system was set up. In a system where only gravity acts that number should never change, so whatever you see there is integration error. The baseline resets whenever bodies merge or you add one.

Units are game units: a Sun-class star has mass 4000 and the default gravitational constant is 100. Orbit shapes and conservation laws behave as they do in reality, but you cannot convert the on-screen seconds or masses into astronomical values.

Why Euler's method lets planets escape

The simplest integrator, forward Euler, updates position with the current velocity and velocity with the current acceleration, both at once. On an elliptical orbit it adds a little energy every lap. Its errors all point the same way, so they pile up instead of cancelling.

The numbers make the point. Take GM = 1, start at distance 1 with a tangential speed of 0.8 (eccentricity 0.36, period about 3.96) and step with dt = 0.01. Forward Euler gains 22.9% in energy after one orbit, 57% after ten and 89% after a hundred. The orbit swells lap by lap and drifts toward escape. Kick-drift-kick under identical conditions never strays more than 0.021% at any instant, and that ceiling does not grow over a hundred orbits.

A smaller time step does not close that gap. Verlet-type schemes are symplectic: they exactly conserve a slightly modified energy that stays close to the true one, so the error oscillates within each orbit instead of accumulating (Hairer, Lubich and Wanner, 2003). Loup Verlet used the method for molecular dynamics in 1967 for the same reason. His runs lasted hundreds of thousands of steps, and energy could not be allowed to leak. Set this tool's Integrator to Euler (for comparison) and you can watch the Energy drift keep climbing and the same orbit unwind outward.

What pushes the drift number up

Verlet still loses accuracy when steps are coarse. The tool splits each 1/60 s frame into substeps, and it uses more of them (up to 16) when there are fewer bodies. Raising Time speed makes each substep cover more simulated time. Run the Solar system preset at speed 1 and then at 4, and compare the drift readout.

The worst errors come from close encounters. Because force goes as 1/r², acceleration can change enormously between two steps when bodies skim past each other. Softening replaces r with √(r² + ε²) in the force law, which caps that spike. Inside ε you are no longer simulating Newtonian gravity. If you set softening to zero and fling a planet past a star at point-blank range, you may see it leave at an absurd speed. That slingshot comes from numerical error, not physics.

When a body flies far off (beyond coordinate 50,000) it is deleted. As with merges, shattering, adding and erasing, the energy baseline is reset whenever the body count changes, so the drift readout carries on measured against the bodies that remain.

The three-body problem and the figure-eight orbit

Two bodies follow Kepler orbits with an exact solution. Add a third and there is no general closed-form answer. For most starting conditions the motion is chaotic, so a tiny change in the initial positions leads to completely different paths later on. Throw three stars of the same class at similar speeds and you will often see them tangle for a while, then fling one out while the other two settle into a tighter binary.

There are still some remarkable periodic solutions. In one of them, three equal masses chase each other around a single figure-eight curve. Cristopher Moore found it numerically in 1993, and Alain Chenciner and Richard Montgomery proved it exists with a variational argument in 2000. The widely quoted initial conditions (computed by Carles Simó) in units where G = m = 1 are: two bodies at (±0.97000436, ∓0.24308753), the middle body moving at (−0.93240737, −0.86473146), the outer two each moving at minus half that velocity. The period is about 6.3259.

Integrating those conditions with dt = 0.001 shows the same contrast as before. With kick-drift-kick the middle body is back at the origin after ten periods and the energy error sits around 10⁻¹¹%. With forward Euler the energy is off by 3% after one period and by 20% after ten, and the figure eight falls apart. The Figure-eight (three-body) preset scales these initial conditions to the tool's units: each star has mass M = 4000, the length unit L is 0.42 times half the shorter side of the view (larger on small screens, so the stars never touch even at closest approach), velocities are multiplied by √(G·M/L), and the period becomes 6.3259 × L ÷ √(G·M/L). Switch the integrator to Euler and reload the preset to watch the figure eight fall apart.

The collision model is an approximation

In Realistic mode the ratio of impact energy to binding energy decides between merging, a hit-and-run graze, partial disruption and total destruction. Mass, momentum and the centre of mass are conserved. The number and size of fragments and the mass thresholds for supernovae and kilonovae are game choices, though. Research codes for planet formation and giant impacts use methods such as smoothed-particle hydrodynamics, which this is not.

Collisions are detected by assuming both bodies move in a straight line during the substep and checking whether their closest approach along that path came within the sum of their radii (a segment closest-point test). So even a small, very fast body cannot jump clean over another in one substep. When a collision is caught after the bodies have already passed each other, they are moved back to their closest-approach arrangement, keeping the centre of mass fixed, so the impact direction is right.

Worked examples

How to use it

1. Load a preset (Solar system, Binary star, Protoplanetary disk and others) or start from an empty sky.
2. Choose a body type, then drag and release to throw it. A quick tap places it at circular-orbit speed around the heaviest body.
3. Change the gravitational constant, time speed, softening and integrator (Velocity Verlet or Euler), and watch how the orbits and the Energy drift readout respond.
4. Switch to 3D to see systems whose orbital planes are tilted.

Key terms

Symplectic integrator

An integrator that preserves phase-space volume. Its energy error oscillates instead of accumulating. Velocity Verlet and leapfrog are the standard examples.

Softening (ε)

Using √(r² + ε²) instead of r in the force law so that close passes do not produce near-infinite forces.

Choreography

A periodic solution in which several bodies follow the same closed curve, evenly spaced in time. The figure eight is the famous first example.

Frequently asked questions

Notes

In 2D mode every body lies in the z = 0 plane and the vertical force is exactly zero. Switching from 3D to 2D flattens bodies onto that plane, which changes their orbits.