Playground · research instrument
In 2000, Chenciner and Montgomery proved three equal masses can chase each other forever around a figure-eight. It exists in real gravity — this page integrates it live with a symplectic scheme accurate enough to hold the dance for hundreds of orbits, then lets you nudge it and watch how long the choreography survives.
Independent research instrument — not claimed as MakerPortal shipped product code. Methods, equations, assumptions, and limitations are disclosed so you can inspect what the page does and does not establish.
Running. Nudge a body and a faint ghost keeps dancing the unperturbed choreography for comparison.
Newtonian gravity between three bodies has no general closed-form solution — this is the original chaos problem. Ordinary integrators (like RK4) leak energy steadily, so an orbit that should repeat forever spirals in or escapes as pure numerical artifact. A symplectic integrator instead conserves a slightly-perturbed energy exactly, so the error stays bounded for hundreds of periods instead of accumulating.
The figure-eight is special: it is a proven periodic solution, and mildly stable — tiny perturbations wobble around it rather than destroying it instantly. The Lagrange triangle is its opposite: an exact solution that is linearly unstable for equal masses. The integration here is clean enough that the triangle holds until you seed the instability yourself — nudge it, then watch the same |δv| that barely dents the figure-eight tear the triangle apart. That asymmetry is real physics, not integration error.
Integration scheme
Yoshida 4th-order composition of velocity Verlet, dt = 2·10⁻³
Three Verlet sub-steps with weights (w₁, w₀, w₁), w₀ = −∛2·w₁. Validated: figure-eight returns to its start after one period (T ≈ 6.3259) to within 10⁻⁴, energy drift ~10⁻¹¹ over 50 periods.
Collision safety
r² → r² + ε², ε² = 10⁻⁶
Plummer softening caps the force during close encounters so the integration never blows up. At choreography distances (~1) it shifts positions by less than 10⁻⁴; the min-separation readout tells you when it is actually active.
Conserved quantities
E = Σ ½mv² − Σ mᵢmⱼ/rᵢⱼ, P = Σ mv, L = Σ m(r × v)
All three are monitored live and shown without smoothing. Initial conditions are shifted to the center-of-momentum frame, so |P| sits at machine precision until a nudge adds momentum.
Honest limitations
Fixed-step integration loses accuracy during very close encounters (the softening is what keeps it finite); a production ephemeris would switch to adaptive or regularized steps. Everything is planar (2D), and the ghost comparison shares the same integrator — it measures sensitivity to the nudge, not absolute truth.
Five preset systems, one symplectic integrator, and a ghost twin that never gets nudged. Here is what each piece of the simulation actually does under the hood.
Yoshida 4th-order symplectic composition. Three Verlet sub-steps per full step with weights [w₁, w₀, w₁] where w₀ = −∛2 · w₁, w₁ = 1/(2 − ∛2). Each sub-step is a standard velocity Verlet: half-kick, drift, half-kick. The composition cancels error terms up to O(dt⁴) while preserving the symplectic structure. dt = 2 × 10⁻³ in dimensionless time units.
Force accumulation. double-nested loop over body pairs, accumulating ax[i], ay[i] via Newton's law: aᵢ = Σⱼ mⱼ (rⱼ − rᵢ) / (|rⱼ − rᵢ|² + ε²)^(3/2). Plummer softening with ε² = 10⁻⁶ caps forces at close approach. The acceleration arrays are zeroed fresh each step via Float64Array.fill(0) — no floating state between steps.
Center-of-momentum framing. Initial conditions are shifted so Σ mᵢv⃗ᵢ = 0 and Σ mᵢr⃗ᵢ = 0. This keeps the view stable around the origin — without it, the whole system would drift slowly across the canvas. The |P| metric readout confirms that momentum stays at machine precision (~10⁻¹⁵) until you nudge a body, which adds the nudge size to |P|.
Trails with age fade. Positions are recorded every 4th step, maximum 700 points per body. Each trail line is drawn with increasing opacity from oldest (0.05) to newest (0.55). This gives a natural motion-blur effect — the denser the arc, the slower the body was moving through that region.
Nudge + ghost. Nudge button stores a deep copy of the current system as the ghost, then adds a random-direction velocity impulse of magnitude |δv| to one randomly selected body. The ghost continues integrating with the same Yoshida step but never gets nudged. Ghost divergence is the RMS separation between the two systems — starts near zero and grows after the nudge.
Frame-rate-independent stepping. The rAF loop accumulates elapsed real time in an accumulator. Each frame processes floor(accumulator / DT) steps, capping at 120 to prevent spiral of death. Time scale slider multiplies accumulation rate — 4× runs physics 4× faster without changing DT, preserving integration accuracy.
Verlet sub-step within Yoshida composition
One Verlet sub-step costs 2 force evaluations. Three sub-steps per Yoshida iteration costs 6 per full step — about twice the cost of RK4 (4 evaluations) but with energy conservation that makes it the right choice for any integration beyond a few orbits.
Each body feels the softened Newtonian pull of the others (G = 1, softening ε² = 10⁻⁶). The energy below is the exact potential of that softened force, so the conserved-quantity readouts test the same law the integrator steps.
Acceleration
Energy
Angular momentum
Both are integrated when the page is built, with the Yoshida scheme the simulation runs (dt = 2·10⁻³, G = 1), and checked in the site's test suite against the published figure-eight values and an independent RK4 integration. Lengths, times and energies are in G = 1 units.
A · One period of the figure-eight
Three unit masses from the preset, run for one period, T = 6.32591398, in 3163 steps.
The energy moves by at most 1.1e-11 (relative) inside the period and ends 1.8e-14 from E₀; the bodies come back to within 1.3e-5 of their starting positions.
B · A 10⁻³ nudge, twenty periods later
The same start, run for 20 periods (t = 126.5), with |δv| = 0.001 added to body 2 along +x, against the unperturbed twin.
The impulse adds momentum |δv| to a system of total mass 3, so the whole figure-eight drifts at |δv|/3: 0.042 of the 0.099 RMS separation is that drift. Repeating the run with a 10⁻² nudge gives 0.96.
Orbital mechanics stack · 7 picks
$39.99ApparatusConservation of momentum/energy demo — the same conservation your symplectic integrator preserves while Newtonian gravity chaotically exchanges momentum among three equal masses.
$220.49BookCanonical derivation of the double-pendulum equations of motion via Euler-Lagrange — exactly the T-V Lagrangian this simulator integrates symplectically.
$86.00BookDirect N-body with Hermite integration — compares to this page's symplectic scheme showing why 4th-order composition keeps figure-eight orbit stable vs Euler blow-up.
$69.47BookChapter 6 explicitly integrates the double pendulum Lagrangian — Lyapunov exponent, phase-space Poincaré section, and the sensitivity that this simulator paints as geometry.
$99.95BookDerives n-body equations, figure-eight choreography existence proof, and symplectic energy conservation — exactly the trajectory this page verifies stays bounded for 500 periods.
BookSGP4 and TLE propagation — this globe computes ISS position client-side from real TLE elements using the same Kepler + J2 math Montenbruck derives.
$399.00Model600+ precision metal parts, motor-driven planets on brass gear trains — a buildable mechanical computer tracing the same Newtonian choreography this page integrates with symplectic Yoshida 4th order.
Prices shown were retrieved from the Amazon Product Advertising API on 8 October 2026 and are indicative only — the price and availability on Amazon at the time of purchase apply.
More gear across every app: the full Gear list →
Fixed-step integrator cannot handle true close encounters
A real three-body close encounter (two bodies passing within 0.001 distance units at high speed) changes velocity by orders of magnitude within a single timestep. The Plummer softening prevents infinity, but accuracy degrades dramatically — the min-separation readout lets you monitor when this is happening. Production ephemeris codes (Mercury, REBOUND) switch to adaptive or regularized steps during close encounters. This demo keeps a fixed dt for clarity and lets the softening handle the edge case.
Yoshida composition weights are exact algebraic numbers
The weights w₀ = −∛2/(2−∛2) and w₁ = 1/(2−∛2) are not chosen for numerical convenience — they are the unique solution to the order conditions for a 3-stage symplectic composition. If you tweak them by even 0.1%, the symplectic property is lost and the energy drift jumps from 10⁻¹¹ to 10⁻⁴. Always use these exact constants. The script computes them as Math.cbrt(2) to machine precision.
The exact integration loop powering this simulation. Port to any language — the core is three Verlet sub-steps and a pair-loop force calculation.
JavaScript — Yoshida step with force accumulation
const DT = 2e-3, EPS2 = 1e-6;
const CBRT2 = Math.cbrt(2);
const W1 = 1 / (2 - CBRT2);
const W0 = -CBRT2 / (2 - CBRT2);
const WEIGHTS = [W1, W0, W1];
function accumulate(bodies, ax, ay) {
ax.fill(0); ay.fill(0);
for (let i = 0; i < bodies.length; i++) {
for (let j = i + 1; j < bodies.length; j++) {
const dx = bodies[j].x - bodies[i].x;
const dy = bodies[j].y - bodies[i].y;
const r2 = dx*dx + dy*dy + EPS2;
const inv = 1 / (r2 * Math.sqrt(r2));
ax[i] += bodies[j].m * dx * inv;
ay[i] += bodies[j].m * dy * inv;
ax[j] -= bodies[i].m * dx * inv;
ay[j] -= bodies[i].m * dy * inv;
}
}
}
function yoshidaStep(bodies) {
for (const w of WEIGHTS) {
accumulate(bodies, scratchAX, scratchAY);
for (let i = 0; i < bodies.length; i++) {
bodies[i].vx += 0.5 * w * DT * scratchAX[i];
bodies[i].vy += 0.5 * w * DT * scratchAY[i];
bodies[i].x += w * DT * bodies[i].vx;
bodies[i].y += w * DT * bodies[i].vy;
}
accumulate(bodies, scratchAX, scratchAY);
for (let i = 0; i < bodies.length; i++) {
bodies[i].vx += 0.5 * w * DT * scratchAX[i];
bodies[i].vy += 0.5 * w * DT * scratchAY[i];
}
}
}Figure-eight initial conditions (Chenciner-Montgomery)
// Period T ≈ 6.32591398, equal masses m = 1
const figureEight = [
{ m: 1, x: 0.97000436, y: -0.24308753, vx: 0.4662036850, vy: 0.4323657300 },
{ m: 1, x: -0.97000436, y: 0.24308753, vx: 0.4662036850, vy: 0.4323657300 },
{ m: 1, x: 0, y: 0, vx: -0.93240737, vy: -0.86473146 },
];
// Shift to center-of-momentum frame before integrating:
// vx[i] -= sum(m*vx)/sum(m), vy[i] -= sum(m*vy)/sum(m)The Chenciner-Montgomery figure-eight is a proven periodic solution — if you place the three bodies at exactly the correct initial conditions, they trace the same loop forever. It is mildly stable: a small nudge does not tear it apart at once. Measured with this page's integrator, a 10⁻³ velocity nudge to one body leaves it 0.099 (RMS position) from the unperturbed twin after 20 periods, and a 10⁻² nudge leaves it 0.96. For small nudges the separation grows roughly in proportion to |δv| × time rather than exponentially. The nudge slider lets you compare the sizes yourself.
The equal-mass Lagrange triangle is an exact solution — three bodies at the vertices of an equilateral triangle orbiting the center at constant angular velocity — but it is linearly unstable. Small perturbations grow exponentially with a Lyapunov timescale of roughly one orbital period. The figure-eight is in a narrow stable island of phase space; the Lagrange triangle sits on an unstable ridge. This is not integration error — it is real celestial mechanics. The "ghost" comparison proves it.
A symplectic integrator exactly conserves a slightly-perturbed energy instead of an approximation of the true energy. Ordinary Runge-Kutta (RK4) leaks energy monotonically — orbits spiral in or out as a pure numerical artifact. The Yoshida 4th-order composition used here has energy drift ~10⁻¹¹ per period over 50 periods for the figure-eight. That means the choreography stays intact for hundreds of orbits before numerical error accumulates, whereas RK4 would visibly distort it within 5-10 orbits.
Hard-sphere collisions (bouncing bodies off each other at zero distance) create a singularity: the force goes to infinity as r → 0, and any fixed-step integrator will diverge. Plummer softening replaces r⁻¹ with (r² + ε²)⁻¹/², capping the force at a finite maximum near ε. Here ε² = 10⁻⁶ means the cap only activates at separations < 0.001 — about 4 orders of magnitude below the typical choreography scale. The integrator stays stable during close encounters while the physics remains accurate everywhere else.
Energy drift = (current_E − initial_E) / |initial_E| — should stay near zero for the symplectic integrator, shown in scientific notation. Momentum |P| tracks the total linear momentum after center-of-mass correction — at machine precision (~10⁻¹⁵) until you nudge a body, which adds |δv|. Angular momentum drift |L − L₀| measures the fidelity of rotational conservation. Ghost divergence = RMS separation between the nudged system and the unperturbed ghost — zero until you nudge, then growing. The nudge also adds momentum |δv| to the system (total mass 3), so after a nudge |P| equals the nudge size, not machine precision, and the whole system drifts at |δv|/3; the two worked examples below separate that drift from the change in shape.
With G = 1, equal masses of 1 and this page's initial conditions (softening ε² = 10⁻⁶), the total energy is E = -1.28714 and the total angular momentum is L = 0, as is the linear momentum after the shift to the centre-of-momentum frame. The period is T ≈ 6.32591398. Over one period the integrator moves the energy by at most 1.1e-11 (relative) and returns the three bodies to within 1.3e-5 (RMS position) of where they started.
Shareable still
This 16:9 frame is rendered from the real browser instrument above. It is the page's canonical preview for image search, link unfurls, and posts that need to show what the tool actually does.
Download 1280 × 720 JPEG