Skip to content
Computational Lab.

N-Body Gravity

Many bodies, one law.

Advanced / particles / gravity / integration

Every body attracts every other with the same inverse-square law, and the resulting motion is integrated step by step. Nothing is scripted: orbits, slingshots, and ejections are all consequences of the arithmetic. The integrator you choose decides whether the system keeps its energy or quietly drifts.

Tagged: orbital mechanics, integration, symplectic, chaos

Live simulation

Hold the pointer to attract bodies with a temporary mass.

Speed

Hold the pointer to attract bodies with a temporary mass.

Mathematics

Every pair of bodies attracts along their connecting line with an inverse-square law. There is no closed-form solution for N ≥ 3, so the motion exists only as the output of an integrator.

a_i = G · Σ_{j≠i} m_j (r_j − r_i) / (|r_j − r_i|² + ε²)^{3/2}

leapfrog (KDK): v += a·dt/2; x += v·dt; a(x); v += a·dt/2

energy: E = Σ ½ m v² − Σ_{i<j} G m_i m_j / r_ij

The energy drift metric is the honest test of an integrator. Forward Euler gains energy every orbit and the disk flies apart; leapfrog is symplectic, so its energy error stays bounded instead of growing; RK4 is the most accurate per step but costs four force evaluations. Switch integrators and watch the drift number, not the picture.

How it works

State is a set of parallel Float32Arrays: positions, velocities, and accelerations for up to 400 bodies plus one heavy central mass. A Plummer softening term ε prevents the singular close encounter from producing infinite forces. Initial conditions place each body on a near-circular Keplerian orbit (v = √(GM/r)) with a small dispersion.

The force loop is O(N²) with no allocations. Orbit trails come from a ring buffer per body, and previous positions enable tick-to-tick interpolation. Holding the pointer adds a temporary mass you can use to perturb or capture bodies.

Implementation

src/engine/experiments/nbody.ts · 361 lines · the module running above

TypeScript · fixed-timestep engine

1import { BaseExperiment } from "../core/base-experiment";
2import type { ParameterDef, PointerState, ThemeColors } from "../core/types";
3
4const TRAIL_LENGTH = 48;
5
6/**
7 * N-body gravity with three interchangeable integrators (forward Euler,
8 * leapfrog KDK, classical RK4) so energy behavior is measurable, not claimed.
9 * Pointer drag inserts a temporary attractor.
10 */
11export class NBodyExperiment extends BaseExperiment {
12 readonly id = "n-body";
13
14 private count = 0;
15 private x = new Float32Array(0);
16 private y = new Float32Array(0);
17 private vx = new Float32Array(0);
18 private vy = new Float32Array(0);
19 private px = new Float32Array(0);
20 private py = new Float32Array(0);
21 private ax = new Float32Array(0);
22 private ay = new Float32Array(0);
23 private mass = new Float32Array(0);
24 private trailX = new Float32Array(0);
25 private trailY = new Float32Array(0);
26 private trailHead = 0;
27 private trailTick = 0;
28 private energy0 = 1;
29 private energy = 1;
30 private time = 0;
31 private pointer = { x: 0, y: 0, down: false, inside: false };
32 // RK4 scratch: position-derivatives (p*) and accelerations (a*) per stage.
33 private k1px = new Float32Array(0);
34 private k1py = new Float32Array(0);
35 private k1ax = new Float32Array(0);
36 private k1ay = new Float32Array(0);
37 private k2px = new Float32Array(0);
38 private k2py = new Float32Array(0);
39 private k2ax = new Float32Array(0);
40 private k2ay = new Float32Array(0);
41 private k3px = new Float32Array(0);
42 private k3py = new Float32Array(0);
43 private k3ax = new Float32Array(0);
44 private k3ay = new Float32Array(0);
45 private k4ax = new Float32Array(0);
46 private k4ay = new Float32Array(0);
47 private tx = new Float32Array(0);
48 private ty = new Float32Array(0);
49
50 private readonly G = 6000;
51 private readonly centralMass = 140;
52 private readonly softening = 6;
53
54 protected params(): ParameterDef[] {
55 return [
56 { key: "bodies", label: "Orbiting bodies", min: 20, max: 400, step: 20, defaultValue: 160 },
57 {
58 key: "integrator",
59 label: "Integrator",
60 options: [
61 { value: "leapfrog", label: "Leapfrog (KDK)" },
62 { value: "rk4", label: "Runge-Kutta 4" },
63 { value: "euler", label: "Forward Euler" },
64 ],
65 defaultValue: "leapfrog",
66 },
67 {
68 key: "trails",
69 label: "Orbit trails",
70 options: [
71 { value: "on", label: "On" },
72 { value: "off", label: "Off" },
73 ],
74 defaultValue: "on",
75 },
76 ];
77 }
78
79 protected onReset(): void {
80 this.count = this.num("bodies");
81 const n = this.count + 1; // index 0 is the central mass
82 this.x = new Float32Array(n);
83 this.y = new Float32Array(n);
84 this.vx = new Float32Array(n);
85 this.vy = new Float32Array(n);
86 this.px = new Float32Array(n);
87 this.py = new Float32Array(n);
88 this.ax = new Float32Array(n);
89 this.ay = new Float32Array(n);
90 this.mass = new Float32Array(n);
91 this.trailX = new Float32Array(n * TRAIL_LENGTH);
92 this.trailY = new Float32Array(n * TRAIL_LENGTH);
93 this.trailHead = 0;
94 this.trailTick = 0;
95 this.k1px = new Float32Array(n);
96 this.k1py = new Float32Array(n);
97 this.k1ax = new Float32Array(n);
98 this.k1ay = new Float32Array(n);
99 this.k2px = new Float32Array(n);
100 this.k2py = new Float32Array(n);
101 this.k2ax = new Float32Array(n);
102 this.k2ay = new Float32Array(n);
103 this.k3px = new Float32Array(n);
104 this.k3py = new Float32Array(n);
105 this.k3ax = new Float32Array(n);
106 this.k3ay = new Float32Array(n);
107 this.k4ax = new Float32Array(n);
108 this.k4ay = new Float32Array(n);
109 this.tx = new Float32Array(n);
110 this.ty = new Float32Array(n);
111
112 this.x[0] = this.width / 2;
113 this.y[0] = this.height / 2;
114 this.mass[0] = this.centralMass;
115 for (let i = 1; i < n; i++) {
116 const radius = 40 + this.rng() * Math.min(this.width, this.height) * 0.42;
117 const angle = this.rng() * Math.PI * 2;
118 this.x[i] = this.x[0] + Math.cos(angle) * radius;
119 this.y[i] = this.y[0] + Math.sin(angle) * radius;
120 // Circular orbit speed for the central mass, plus a small dispersion.
121 const v = Math.sqrt((this.G * this.centralMass) / radius) * (0.9 + this.rng() * 0.2);
122 this.vx[i] = -Math.sin(angle) * v;
123 this.vy[i] = Math.cos(angle) * v;
124 this.mass[i] = 1;
125 this.px[i] = this.x[i];
126 this.py[i] = this.y[i];
127 }
128 this.time = 0;
129 this.energy0 = this.totalEnergy();
130 this.energy = this.energy0;
131 }
132
133 protected onParameterChange(key: string): void {
134 if (key === "bodies" || key === "integrator") this.reset();
135 }
136
137 resize(width: number, height: number): void {
138 this.width = width;
139 this.height = height;
140 this.reset();
141 }
142
143 onPointer(state: PointerState): void {
144 this.pointer = state;
145 }
146
147 /** Pairwise accelerations with Plummer softening; O(n^2) over alive bodies. */
148 private accelerations(px: Float32Array, py: Float32Array, ax: Float32Array, ay: Float32Array): void {
149 const n = this.count + 1;
150 ax.fill(0);
151 ay.fill(0);
152 const eps2 = this.softening * this.softening;
153 // Temporary attractor under pointer drag.
154 const pointerMass = this.pointer.down && this.pointer.inside ? 400 : 0;
155 for (let i = 0; i < n; i++) {
156 let axi = 0;
157 let ayi = 0;
158 for (let j = 0; j < n; j++) {
159 if (j === i) continue;
160 const dx = px[j] - px[i];
161 const dy = py[j] - py[i];
162 const r2 = dx * dx + dy * dy + eps2;
163 const inv = 1 / (r2 * Math.sqrt(r2));
164 const s = this.G * this.mass[j] * inv;
165 axi += dx * s;
166 ayi += dy * s;
167 }
168 if (pointerMass > 0) {
169 const dx = this.pointer.x - px[i];
170 const dy = this.pointer.y - py[i];
171 const r2 = dx * dx + dy * dy + 900;
172 const inv = 1 / (r2 * Math.sqrt(r2));
173 const s = this.G * pointerMass * inv;
174 axi += dx * s;
175 ayi += dy * s;
176 }
177 ax[i] = axi;
178 ay[i] = ayi;
179 }
180 }
181
182 protected onUpdate(dt: number): void {
183 const integrator = this.str("integrator");
184 if (integrator === "leapfrog") this.stepLeapfrog(dt);
185 else if (integrator === "rk4") this.stepRK4(dt);
186 else this.stepEuler(dt);
187 this.time += dt;
188 this.trailTick += 1;
189 if (this.trailTick % 2 === 0) {
190 const n = this.count + 1;
191 for (let i = 0; i < n; i++) {
192 this.trailX[i * TRAIL_LENGTH + this.trailHead] = this.x[i];
193 this.trailY[i * TRAIL_LENGTH + this.trailHead] = this.y[i];
194 }
195 this.trailHead = (this.trailHead + 1) % TRAIL_LENGTH;
196 }
197 if (this.trailTick % 20 === 0) this.energy = this.totalEnergy();
198 }
199
200 private stepEuler(dt: number): void {
201 this.accelerations(this.x, this.y, this.ax, this.ay);
202 const n = this.count + 1;
203 for (let i = 0; i < n; i++) {
204 this.px[i] = this.x[i];
205 this.py[i] = this.y[i];
206 this.vx[i] += this.ax[i] * dt;
207 this.vy[i] += this.ay[i] * dt;
208 this.x[i] += this.vx[i] * dt;
209 this.y[i] += this.vy[i] * dt;
210 }
211 }
212
213 private stepLeapfrog(dt: number): void {
214 const n = this.count + 1;
215 // Kick-drift-kick: recompute a at the midpoint for bounded energy error.
216 this.accelerations(this.x, this.y, this.ax, this.ay);
217 for (let i = 0; i < n; i++) {
218 this.vx[i] += this.ax[i] * dt * 0.5;
219 this.vy[i] += this.ay[i] * dt * 0.5;
220 this.px[i] = this.x[i];
221 this.py[i] = this.y[i];
222 this.x[i] += this.vx[i] * dt;
223 this.y[i] += this.vy[i] * dt;
224 }
225 this.accelerations(this.x, this.y, this.ax, this.ay);
226 for (let i = 0; i < n; i++) {
227 this.vx[i] += this.ax[i] * dt * 0.5;
228 this.vy[i] += this.ay[i] * dt * 0.5;
229 }
230 }
231
232 private stepRK4(dt: number): void {
233 const n = this.count + 1;
234 // Classical RK4 on the first-order system (x' = v, v' = a(x)).
235 // Stage 1: a(x)
236 this.accelerations(this.x, this.y, this.k1ax, this.k1ay);
237 for (let i = 0; i < n; i++) {
238 this.k1px[i] = this.vx[i];
239 this.k1py[i] = this.vy[i];
240 }
241 // Stage 2: a(x + dt/2 * k1p)
242 for (let i = 0; i < n; i++) {
243 this.tx[i] = this.x[i] + (this.k1px[i] * dt) / 2;
244 this.ty[i] = this.y[i] + (this.k1py[i] * dt) / 2;
245 }
246 this.accelerations(this.tx, this.ty, this.k2ax, this.k2ay);
247 for (let i = 0; i < n; i++) {
248 this.k2px[i] = this.vx[i] + (this.k1ax[i] * dt) / 2;
249 this.k2py[i] = this.vy[i] + (this.k1ay[i] * dt) / 2;
250 }
251 // Stage 3: a(x + dt/2 * k2p)
252 for (let i = 0; i < n; i++) {
253 this.tx[i] = this.x[i] + (this.k2px[i] * dt) / 2;
254 this.ty[i] = this.y[i] + (this.k2py[i] * dt) / 2;
255 }
256 this.accelerations(this.tx, this.ty, this.k3ax, this.k3ay);
257 for (let i = 0; i < n; i++) {
258 this.k3px[i] = this.vx[i] + (this.k2ax[i] * dt) / 2;
259 this.k3py[i] = this.vy[i] + (this.k2ay[i] * dt) / 2;
260 }
261 // Stage 4: a(x + dt * k3p)
262 for (let i = 0; i < n; i++) {
263 this.tx[i] = this.x[i] + this.k3px[i] * dt;
264 this.ty[i] = this.y[i] + this.k3py[i] * dt;
265 }
266 this.accelerations(this.tx, this.ty, this.k4ax, this.k4ay);
267 // Combine: k4p (position derivative) is v + dt * k3a.
268 for (let i = 0; i < n; i++) {
269 this.px[i] = this.x[i];
270 this.py[i] = this.y[i];
271 const dpx =
272 (dt / 6) *
273 (this.k1px[i] + 2 * this.k2px[i] + 2 * this.k3px[i] + (this.vx[i] + this.k3ax[i] * dt));
274 const dpy =
275 (dt / 6) *
276 (this.k1py[i] + 2 * this.k2py[i] + 2 * this.k3py[i] + (this.vy[i] + this.k3ay[i] * dt));
277 const dvx = (dt / 6) * (this.k1ax[i] + 2 * this.k2ax[i] + 2 * this.k3ax[i] + this.k4ax[i]);
278 const dvy = (dt / 6) * (this.k1ay[i] + 2 * this.k2ay[i] + 2 * this.k3ay[i] + this.k4ay[i]);
279 this.x[i] += dpx;
280 this.y[i] += dpy;
281 this.vx[i] += dvx;
282 this.vy[i] += dvy;
283 }
284 }
285
286 private totalEnergy(): number {
287 const n = this.count + 1;
288 let kinetic = 0;
289 let potential = 0;
290 const eps2 = this.softening * this.softening;
291 for (let i = 0; i < n; i++) {
292 kinetic += 0.5 * this.mass[i] * (this.vx[i] ** 2 + this.vy[i] ** 2);
293 for (let j = i + 1; j < n; j++) {
294 const dx = this.x[j] - this.x[i];
295 const dy = this.y[j] - this.y[i];
296 const r = Math.sqrt(dx * dx + dy * dy + eps2);
297 potential -= (this.G * this.mass[i] * this.mass[j]) / r;
298 }
299 }
300 return kinetic + potential;
301 }
302
303 render(ctx: CanvasRenderingContext2D, theme: ThemeColors, alpha: number): void {
304 ctx.fillStyle = theme.bg;
305 ctx.fillRect(0, 0, this.width, this.height);
306 const n = this.count + 1;
307 const showTrails = this.str("trails") === "on";
308 if (showTrails) {
309 ctx.lineWidth = 1;
310 for (let i = 1; i < n; i++) {
311 ctx.strokeStyle = theme.viz[1];
312 ctx.globalAlpha = 0.35;
313 ctx.beginPath();
314 for (let t = 0; t < TRAIL_LENGTH; t++) {
315 const idx = (this.trailHead + t) % TRAIL_LENGTH;
316 const tx = this.trailX[i * TRAIL_LENGTH + idx];
317 const ty = this.trailY[i * TRAIL_LENGTH + idx];
318 if (t === 0) ctx.moveTo(tx, ty);
319 else ctx.lineTo(tx, ty);
320 }
321 ctx.stroke();
322 }
323 ctx.globalAlpha = 1;
324 }
325 for (let i = 0; i < n; i++) {
326 let ix = this.px[i] + (this.x[i] - this.px[i]) * alpha;
327 let iy = this.py[i] + (this.y[i] - this.py[i]) * alpha;
328 if (Math.abs(ix - this.x[i]) > this.width / 2) ix = this.x[i];
329 if (Math.abs(iy - this.y[i]) > this.height / 2) iy = this.y[i];
330 if (i === 0) {
331 ctx.fillStyle = theme.accent;
332 ctx.beginPath();
333 ctx.arc(ix, iy, 4, 0, Math.PI * 2);
334 ctx.fill();
335 } else {
336 ctx.fillStyle = theme.fg;
337 ctx.fillRect(ix - 1, iy - 1, 2, 2);
338 }
339 }
340 }
341
342 getMetrics() {
343 const drift = this.energy0 !== 0 ? ((this.energy - this.energy0) / Math.abs(this.energy0)) * 100 : 0;
344 return {
345 bodies: this.count + 1,
346 integrator: this.str("integrator") === "leapfrog" ? "Leapfrog" : this.str("integrator") === "rk4" ? "RK4" : "Euler",
347 time: `${this.time.toFixed(1)} s`,
348 energyDrift: `${drift >= 0 ? "+" : ""}${drift.toFixed(3)} %`,
349 };
350 }
351
352 describe(): string {
353 const drift = this.energy0 !== 0 ? ((this.energy - this.energy0) / Math.abs(this.energy0)) * 100 : 0;
354 return `${this.count} bodies orbiting a central mass; total energy has drifted ${drift.toFixed(2)}% since launch.`;
355 }
356
357 entities(): number {
358 return this.count + 1;
359 }
360}
361

Performance

Time per tick
O(N²) force pairs per step; 4× for RK4
Space
Typed arrays plus a 48-sample trail per body
Notes
At 160 bodies a leapfrog step is a few hundred thousand pair interactions, comfortably inside one 60 Hz tick.

Live numbers (FPS, frame time, steps per second, microseconds per step, entity count) are in the performance overlay on the simulation above: the pulse icon in its top-right corner.