Three-Body Problem

About this piece

The three-body problem asks how three masses move when each pulls on the other two by gravity. Two bodies are solved exactly: each goes round an ellipse for ever. Add a third and there is no general formula. The answer has to be computed step by step, and for most starting positions it is chaotic: a tiny change in the start grows until the whole story is different. This plate runs two famous starting positions through that computation and draws what comes out. Units are the ones the papers use, with the gravitational constant G = 1, so a “time unit” is not seconds or years; it is the problem’s own clock.

The loop opens on a still of Act 1, one lap already traced. Act 1 plays one more lap of the figure eight at 0.5 time units per second, about 13 seconds. Act 2 then runs the Pythagorean problem from t = 0 to t = 70 at 2 time units per second, 35 seconds, and the loop starts again. Each act runs at one steady speed throughout. Close approaches are not slowed down: when two bodies swing past each other in Act 2 they really do move that fast, and you see it as a flick.

Act 1: the figure eight

In 2000 Alain Chenciner and Richard Montgomery proved, in the Annals of Mathematics, that three equal masses can chase each other round one fixed eight-shaped curve, with zero total angular momentum. In their words, apart from motions along a straight line, the only other known motion along a fixed curve was Lagrange’s: three bodies at the corners of an equilateral triangle that turns rigidly. The paper’s Figure 1 gives initial conditions computed by Carles Simó, and this plate copies them digit for digit:

x1 = −x2 = ( 0.97000436, −0.24308753)     x3 = (0, 0)
v3 = −2 v1 = −2 v2 = (−0.93240737, −0.86473146)
m1 = m2 = m3 = 1        period T = 6.32591398

Body 3 starts at the crossing point with the other two at the far ends of the lobes, one on each side. Each body follows the one ahead of it a third of a period behind. That is why each tail here is exactly one third of a period long: the three bright tails join end to end and cover the whole eight between them. The paper also reports that Simó’s numerical computations indicate the orbit is stable, so a slightly disturbed start should still trace something close to an eight.

The test suite holds the plate’s own integrator to this. After one published period every body is back within 10⁻⁴ of where it started. The measured gap is 3.0×10⁻⁸, about as close as eight printed digits allow. Energy drifts by less than 10⁻¹² of itself over the lap, and angular momentum stays at zero to 2×10⁻¹⁴.

Act 2: Burrau’s Pythagorean problem

Put masses 3, 4 and 5 at rest on the corners of a right triangle with sides 3, 4 and 5, each mass opposite the side of its own length, and let go. According to a project page on the UC Observatories (Lick) site, the mathematician Meissel proposed this in 1893 and expected the motion to be periodic. Carl Burrau integrated it numerically in 1913 with a variable step size, but over the stretch he could compute he found no sign that it was. It was settled in 1967 by Victor Szebehely and C. Frederick Peters, with high-speed numerical integration and a two-body regularisation that handles near-collisions. The motion is neither periodic nor quasi-periodic. After a long, complicated dance the two heavy bodies end up as a bound pair and the lightest leaves on a hyperbolic path.

The positions used here are the ones tabulated by Roa, Urrutxua and Peláez (MNRAS, 2016), which put the centre of mass at the origin: mass 3 at (1, 3), mass 4 at (−2, −1), mass 5 at (1, −1). A 2020 paper on N-body integrators by Antoñana, Chartier, Makazaga and Murua uses the same numbers as a test problem and credits Szebehely and Peters with the result. Roa et al. describe the outcome: the bodies first fall towards the origin, and after a number of close approaches the mass-3 body “is ejected along a trajectory in the first quadrant”, while masses 4 and 5 “form a binary that escapes in the opposite direction. The escape occurs at approximately t_esc ∼ 60.”

That is what the plate shows, and what its test checks against its own step function. On screen the lightest body is drawn in coral and the pair left behind in the two cobalts. By t = 60 the coral body is unbound from the pair: its energy relative to them is positive, so it is never coming back. It crosses 5 units from the centre of mass at t ≈ 61.3 and is 8.2 units out and still moving away at t = 63, heading up and to the right. It leaves the frame shortly after. The 4–5 pair drifts off towards the lower left, turning about each other as they go.

The closest approach in the whole act comes at t ≈ 15.83, when masses 4 and 5 pass 0.00041 units apart, about a ten-thousandth of the triangle’s size. That moment is why the integrator has to change its step size.

How it is computed

Each body’s acceleration is the sum of mⁱ(rⁱ − r)/|rⁱ − r|³ over the other two. The equations are stepped forward by the Dormand–Prince 5(4) method. Each step computes a fifth-order and a fourth-order answer, and the gap between them estimates the error. The step is rejected and shortened whenever that estimate is above 10⁻¹² of the size of the numbers involved, and lengthened when there is room. In Act 2 the accepted steps run from 1.0×10⁻⁸ time units at the tightest near-collision up to 0.019 on the quiet stretches: 49,408 steps to reach t = 70. The worst energy error along the way is 1.3×10⁻⁸ of the total. A figure-eight lap takes 1,522 steps.

A fixed step cannot do this. The test suite runs the same Act 2 with a fixed step of 0.01 as a planted defect, and the energy ends up wrong by more than a hundred times its own size. That is a physically impossible result, and the check that catches it goes red, as it should. Tolerance matters even with an adaptive step. At 10⁻¹⁰ this code throws out the wrong body, mass 4. At 10⁻¹¹, 10⁻¹² and 10⁻¹³ it throws out mass 3, as the literature says. The plate uses 10⁻¹².

The frames you see are not the steps. The integrator’s steps are chosen by the error control alone, never by the screen. What is drawn at time t is a cubic curve through the positions and velocities at both ends of the step that contains t. So a 30-frames-per-second phone and a 144 Hz monitor show the identical state at the same moment, bit for bit, which the test suite also checks. The faint line under each body is its whole path so far, sampled every 1/64 of a time unit. The bright part is its recent tail. Each disc is sized by the cube root of its mass.

No colour is named in the drawing code. The background and the four ink slots come from the shared palette, so the gallery tile, this page and an embedded copy follow the same light and dark settings. Under prefers-reduced-motion no animation runs. You get one composed still instead: Act 1 after exactly one period, the whole eight traced and the three bodies back at their starting points.

The two plates it sits nearest

• Double Pendulum is also sensitive dependence you can watch, but the chaos is mechanical: rods and a fixed pivot, integrated with a fixed step of 1/480 of a second. There are no close approaches for a fixed step to fail on. Here nothing is fixed in place, the force grows without limit as two bodies meet, and the step has to adapt.
• Strange Attractor draws a chaotic orbit too, but of an iterated map: a formula applied over and over, with no forces, no masses and no time in between. The three-body plate is a real physical system, and it shows the two faces of that system side by side: one perfectly periodic orbit, and one that ends by throwing a body out.

Honest limits

Point masses in a plane. The bodies have no size, cannot collide and never leave the plane, and gravity is Newton’s with no relativity. The discs are drawn at a size you can see, not to scale.

The ending depends on accuracy. Act 2 is chaotic, and the tolerance experiment above shows its outcome is only as good as the integration. Agreement with a 1967 result and at three tolerances is strong evidence, not proof.

The frame is fixed. The view is set to hold the dance before the escape, so the escaping body and then the pair leave the picture before t = 70. The figure eight is shown for two laps only, and every loop restarts from the published numbers, so error never builds up from one loop to the next.

Nothing here is adjustable. You cannot set your own masses or positions. This gallery ships finished compositions rather than a simulator.

Draw your own

The figure eight is short enough to fly yourself. The snippet below is Act 1 of this plate in 30 lines of HTML and plain JavaScript with no library, and with the plate’s own model: Newton’s gravity with G = 1, three unit masses, Simó’s eight-digit starting positions and velocities from Figure 1 of Chenciner and Montgomery’s paper (arXiv math/0011268), and the same Dormand–Prince 5(4) integrator with the same tolerances and the same step-size rule. There is no fixed time step to copy, because the plate has none: its step is whatever the error control allows. The snippet makes the same choices in the same order, so it takes the same 1,522 steps per lap, and the test suite checks that they land on the plate’s numbers to the last bit. It draws Act 1 the way the plate does: a path point every 1/64 of a time unit, a cubic across the step for the moments in between, time running at 0.5 units per second, a faint line for the lap and a bright tail one third of a period long. What it leaves out is everything that exists for the gallery rather than for the physics: Act 2 (Burrau’s problem), the fades, the restart after two laps (this one keeps flying and keeps only the last lap’s line), and sharp drawing on high-density screens (it draws in CSS pixels).

• State. Y0 is twelve numbers: the x and y of the three bodies, then their velocities. The digits are the paper’s: x1 = −x2 = (0.97000436, −0.24308753), x3 = (0, 0), v3 = (−0.93240737, −0.86473146) and period T = 6.32591398. The paper gives the other two velocities as v1 = v2 = −v3/2, which is (0.466203685, 0.43236573) exactly, so the total momentum starts at zero.
• Equations and integrator. deriv() adds up, for each body, G mⁱ(rⁱ − r)/|rⁱ − r|³ from the other two. step() is one Dormand–Prince step: six new calls to deriv(), a fifth-order answer, and the gap to the fourth-order one as the error estimate. A try is rejected and shrunk while that gap is above RTOL = 10⁻¹² of the numbers (plus ATOL = 10⁻¹⁴); an accepted step sets the next one to 0.9 × error^−1/5 times its own size, kept between 0.2 and 5 times. Measured on this code: after one period every body is back within 4.1×10⁻⁸ of where it started, about as close as eight printed digits allow; energy stays within 7.9×10⁻¹³ of itself over the lap; and the total momentum never rises above 6.3×10⁻¹⁵.
• Step and speed. A lap is 1,522 accepted steps. The first trial step is 0.001 time units and the longest step is 0.0057. A whole lap takes about 10 ms in Node on the machine that built this page. On screen a lap lasts T / 0.5 = 12.7 seconds, so a 60 fps screen takes about two steps per frame. at() steps until it has passed the time asked for, then reads positions off the cubic through both ends of that step, so the frame rate never changes the steps.
• Colour. The ground and the three inks are this site’s palette, light or dark to match your device. Body 1 is coral, body 2 the deeper cobalt and body 3 cobalt, in the plate’s order. The opacities for line, tail and disc are the plate’s too: 0.24, 0.95 and 1 on dark; 0.3, 1 and 1 on light. Reduced-motion visitors get one still frame, like the plate’s: exactly one period computed and drawn once, the whole eight traced and the three bodies back where they began.

A few edits show how special these numbers are. Each result was measured on this code. Set G = 2 and the eight never closes: after one period the bodies are 0.98 units from their starting points, they tangle, and by t = 60 one of them is 13 units from the middle. At G = 0.5 the three spread apart and are 29 units out by t = 60. Make one body 1% heavier with M = [1, 1, 1.01] and the lap misses its start by 0.16, but the bodies keep to a wobbling eight: none goes further than 1.32 units from the middle in 60 time units, against 1.08 for the true orbit. That fits the paper’s report that Simó’s computations indicate the orbit is “stable”. Loosen RTOL to 1e-6 and the gap after one lap grows from 4.1×10⁻⁸ to 2.3×10⁻⁶.

<canvas id="threebody" style="position:fixed; inset:0"></canvas>
<script>
const dark = matchMedia('(prefers-color-scheme: dark)').matches, still = matchMedia('(prefers-reduced-motion: reduce)').matches;
const [ground, inks, LOOK] = dark ? ['11,13,18', ['255,120,84', '120,158,255', '150,180,255'], [0.24, 0.95, 1]]
                                  : ['231,226,213', ['190,68,28', '52,88,214', '40,72,205'], [0.3, 1, 1]];  // bodies 1-3: coral, deep cobalt, cobalt; LOOK: path, tail, disc opacity
const G = 1, M = [1, 1, 1], T = 6.32591398, RATE = 0.5, SAMPLE = 1 / 64, RTOL = 1e-12, ATOL = 1e-14;  // T and Y0: Chenciner & Montgomery (2000), Figure 1, digits by Carles Simó
const Y0 = [0.97000436, -0.24308753, -0.97000436, 0.24308753, 0, 0, 0.466203685, 0.43236573, 0.466203685, 0.43236573, -0.93240737, -0.86473146];  // x1 y1 x2 y2 x3 y3, then vx1 vy1 ... vy3
function deriv(y) { const o = [...y.slice(6), 0, 0, 0, 0, 0, 0];  // Newton: body i is pulled by G m_j (r_j - r_i) / |r_j - r_i|^3 towards each other body j
  for (let i = 0; i < 3; i++) for (let j = i + 1; j < 3; j++) { const dx = y[2 * j] - y[2 * i], dy = y[2 * j + 1] - y[2 * i + 1], r2 = dx * dx + dy * dy, f = G / (r2 * Math.sqrt(r2));
    o[6 + 2 * i] += M[j] * dx * f; o[7 + 2 * i] += M[j] * dy * f; o[6 + 2 * j] -= M[i] * dx * f; o[7 + 2 * j] -= M[i] * dy * f; } return o; }
const A = [[], [1 / 5], [3 / 40, 9 / 40], [44 / 45, -56 / 15, 32 / 9], [19372 / 6561, -25360 / 2187, 64448 / 6561, -212 / 729], [9017 / 3168, -355 / 33, 46732 / 5247, 49 / 176, -5103 / 18656], [35 / 384, 0, 500 / 1113, 125 / 192, -2187 / 6784, 11 / 84]];
const B5 = [...A[6], 0], B4 = [5179 / 57600, 0, 7571 / 16695, 393 / 640, -92097 / 339200, 187 / 2100, 1 / 40];  // Dormand-Prince 5(4) weights, 5th and 4th order
let t0 = 0, y0 = Y0, f0 = deriv(Y0), t1 = 0, y1 = Y0, f1 = f0, h = 1e-3;  // the accepted step runs from (t0, y0) to (t1, y1); f = deriv(y)
function step() { for (;;) { const K = [f1], ys = [y1];  // one accepted step: retried smaller while the 5th-4th order gap is above RTOL of the numbers
    for (let s = 1; s < 7; s++) { ys[s] = y1.map((v, i) => A[s].reduce((a, c, j) => a + h * c * K[j][i], v)); K[s] = deriv(ys[s]); }
    const err = Math.max(...ys[6].map((v, i) => Math.abs(h * K.reduce((a, k, s) => a + (B5[s] - B4[s]) * k[i], 0)) / (ATOL + RTOL * Math.max(Math.abs(y1[i]), Math.abs(v)))));
    if (err <= 1) { [t0, y0, f0, t1, y1, f1] = [t1, y1, f1, t1 + h, ys[6], K[6]]; h *= Math.min(5, Math.max(0.2, err > 0 ? 0.9 * err ** -0.2 : 5)); return; } h *= Math.max(0.2, 0.9 * err ** -0.2); } }
function at(t) { while (t1 < t) step(); const d = t1 - t0, u = d > 0 ? Math.max(0, (t - t0) / d) : 1, u2 = u * u, u3 = u2 * u;  // positions at t: a cubic across the step holding t
  return y1.slice(0, 6).map((_, i) => (2 * u3 - 3 * u2 + 1) * y0[i] + (u3 - 2 * u2 + u) * d * f0[i] + (-2 * u3 + 3 * u2) * y1[i] + (u3 - u2) * d * f1[i]); }
const cv = document.getElementById('threebody'), ctx = cv.getContext('2d'), path = [], LAP = Math.round(T / SAMPLE), TAIL = Math.round(T / 3 / SAMPLE); let n = 0, t = 0, last = performance.now();
function advance(to) { for (; n * SAMPLE <= to + 1e-12; n++) path.push(at(n * SAMPLE)) > LAP && path.shift(); }  // a point every 1/64 time unit, the last lap kept
function draw(now) { const W = cv.width = innerWidth, H = cv.height = innerHeight, S = Math.min(W * 0.46 / 1.12, H), pts = [...path, now], P = (p, i) => [W / 2 + p[2 * i] * S, H / 2 - p[2 * i + 1] * S];
  const ink = (i, a) => 'rgba(' + inks[i] + ',' + a + ')'; ctx.fillStyle = 'rgb(' + ground + ')'; ctx.fillRect(0, 0, W, H); ctx.lineCap = ctx.lineJoin = 'round';
  for (let i = 0; i < 3; i++) { ctx.lineWidth = 1; ctx.strokeStyle = ink(i, LOOK[0]); ctx.beginPath(); pts.forEach((p) => ctx.lineTo(...P(p, i))); ctx.stroke(); }  // the whole lap, faint
  for (let i = 0; i < 3; i++) for (let j = Math.max(1, pts.length - TAIL); j < pts.length; j++) { const q = 1 - (pts.length - 1 - j) / TAIL;  // tails: a third of a period each
    ctx.lineWidth = 0.8 + 2.2 * q; ctx.strokeStyle = ink(i, LOOK[1] * q * q); ctx.beginPath(); ctx.moveTo(...P(pts[j - 1], i)); ctx.lineTo(...P(pts[j], i)); ctx.stroke(); }
  for (let i = 0; i < 3; i++) for (const [R, a] of [[2.4, 0.14], [1, 1]]) { ctx.fillStyle = ink(i, LOOK[2] * a); ctx.beginPath(); ctx.arc(...P(now, i), R * Math.max(2.4, 0.009 * Math.min(W, H)), 0, 2 * Math.PI); ctx.fill(); } }
const frame = (ms) => { t += Math.min(0.1, Math.max(0, ms - last) / 1000) * RATE; last = ms; advance(t); draw(at(t)); requestAnimationFrame(frame); };
if (still) { advance(T); draw(at(T)); } else requestAnimationFrame(frame);  // reduced motion: one still, one full period traced
</script>

Curious how the loop and canvas fit together? Read how it works →