N-body gravity

1Newton's law of gravitation

Every pair of bodies attracts along the line between them with a force proportional to both masses and falling off with the square of the distance. Add up the pulls from everyone else and divide by your own mass to get your acceleration.

A tiny softening length ε stops the force blowing up when two bodies pass through each other. With … bodies there are … pairs to add up on every force evaluation, which is why direct summation costs n².

F = G m₁ m₂ / r²G = 1 here; masses and distances in page unitsfirst planet: r = …, G M / r² = …

aᵢ = Σⱼ G mⱼ (rⱼ − rᵢ) / (|rⱼ − rᵢ|² + ε²)^3/2sum over every other body j|a| = … (sun alone gives …)

2Euler: the obvious step

Move the planet along its current velocity, then change the velocity by the current acceleration. It is first-order: halve Δt and the error per orbit only halves. Worse, on an orbit it always overshoots outwards, so every step adds a little energy and the planet spirals away.

xₙ₊₁ = xₙ + vₙ Δt
vₙ₊₁ = vₙ + a(xₙ) Δtboth updates use the old statefirst planet, x: …drift so far |ΔE / E₀| = …

3Velocity Verlet: half a kick, a drift, half a kick

Give the velocity half a step of acceleration, move with that velocity, recompute the forces at the new place and add the other half kick. It costs the same single force evaluation per step as Euler but is second-order and time-reversible.

It is also symplectic: it keeps the area of phase space exactly, so it follows the true orbit of a slightly different "shadow" system. Its energy wobbles but never wanders off, however long you run it. This is leapfrog, the workhorse of astronomy and molecular dynamics.

v½ = vₙ + ½ a(xₙ) Δt
xₙ₊₁ = xₙ + v½ Δt
vₙ₊₁ = v½ + ½ a(xₙ₊₁) Δtkick, drift, kickfirst planet: v½ = …drift so far |ΔE / E₀| = …

4RK4: four looks before you leap

Classic Runge-Kutta samples the slope four times across the step (start, two midpoints, end) and takes a weighted average. It is fourth-order, so halving Δt cuts the error sixteen times, and it is beautifully accurate over a few orbits.

But it is not symplectic: its tiny error always leans the same way, and over thousands of orbits the energy creeps steadily while Verlet's stays put. It also needs four force evaluations per step, so at the same cost Verlet could take four smaller steps.

k₁ = f(yₙ), k₂ = f(yₙ + ½Δt k₁)
k₃ = f(yₙ + ½Δt k₂), k₄ = f(yₙ + Δt k₃)
yₙ₊₁ = yₙ + Δt (k₁ + 2k₂ + 2k₃ + k₄) / 6y = (x, v), f(y) = (v, a(x))drift so far |ΔE / E₀| = …force evaluations: …

5Energy: the honest referee

A real gravitating system never gains or loses energy, so the total of kinetic and potential energy is a free check on any integrator: whatever it drifts by is pure numerical error. The drift chart plots that error on a log scale, one line per method, all three started from exactly the same state. Angular momentum L is conserved too.

A bound orbit has negative total energy. Euler's error is always positive on an orbit, so it pushes the energy up towards zero, and past it the planet escapes.

E = Σᵢ ½ mᵢ vᵢ² − Σᵢ<ⱼ G mᵢ mⱼ / rᵢⱼkinetic plus potential (Verlet run shown)K = …, U = …E = …, E₀ = …

L = Σᵢ mᵢ (xᵢ vᵧᵢ − yᵢ vₓᵢ)angular momentumL = …

6Barnes-Hut: a tree of boxes

For thousands of stars, n² pairs is too slow. Barnes-Hut puts the stars in a quadtree: split the square into four, and keep splitting any box that holds more than one star. Each box remembers its total mass and centre of mass.

To find the force on a star, walk down from the top. If a box of width s is far enough away, at distance d, that s / d < θ, treat everything inside it as one heavy star at its centre of mass; otherwise open it and look at its four children. θ = 0 opens everything (exact, n²); θ near 1 is fast and a little rough, about n log n. The tree is rebuilt every step in a Web Worker so the page never stalls. Binary search halves a list the same way the tree halves space, and Boids can find its neighbours with a quadtree too, next to its flat grid, for a head-to-head count.

s / d < θopen the box unless it is small for its distanceθ = …; one star's biggest accepted box: s = …, d = …, s / d = …

cost ≈ n log n, not n(n − 1)body-box interactions per step… instead of … (…)tree: … boxes, … levels deep

7Three bodies

Two bodies make an ellipse; three have no general formula and are usually chaotic, sooner or later flinging one out. A few special starts repeat for ever, like the figure-eight found by Moore in 1993 and proved by Chenciner and Montgomery in 2000, where three equal stars chase each other round one loop. Nudge it and watch how long it survives.

The challenge drops three stars at rest: drag from each to give it a velocity, then launch. They stay together only if the total energy is negative; a tight pair with a distant third (a hierarchical triple) is the safest bet. Score is the number of figure-eight periods before a star escapes.

How fragile is an orbit? A shadow copy rides along a hair's breadth (10⁻⁸) away; after every frame the page measures how far it has strayed and pulls it back, and the average log stretch is the Lyapunov exponent λ, the same estimate Chaos lab uses. "λ in Chaos lab" opens the current orbit there, with dozens of copies parting side by side.

x₁ = −x₂ = (0.9700, −0.2431), x₃ = 0
v₁ = v₂ = −½ v₃, v₃ = (−0.9324, −0.8647)the figure-eight start, G = m = 1, period T ≈ 6.326

E < 0 to stay bounda star escapes when ½ μ v² − G m M / r > 0 relative to the other twoE = …, time = … periods, closest pass …

λ ≈ (1 / t) Σ ln(dₖ / d₀)shadow renormalised to d₀ = 10⁻⁸ each frameλ = …