The N-Body Problem in the Browser — Implementing Gravitational Simulation
Go deeper on this topic
Creative Coding
Mathematical simulations, nature visualization, and generative art in the browser
Article 6 of 12 in this series.
Next reads
Exploring the Reaction-Diffusion Parameter Space — Four Faces of the Gray-Scott Model
Varying feed/kill rates in the Gray-Scott Reaction-Diffusion model produces entirely different patterns: spots, stripes, coral, and cell division. Implementation on a 200x200 grid with exploration notes.
Sandpile Model and Self-Organized Criticality — Running the BTW Model in the Browser
Visualizing self-organized criticality in the browser through an implementation of the Bak-Tang-Wiesenfeld sandpile model. Measuring power-law avalanche distributions and the formation of fractal patterns.
Visualizing Double Pendulum Chaos — RK4 Integration and Sensitivity to Initial Conditions
Implementing a double pendulum simulation covering Lagrangian mechanics equation derivation, RK4 integration stability, NaN prevention, and visualization of sensitivity to initial conditions.
Use the topic hub as a map, then continue to the next article in the series. New here? Start at the home hub →
Gravitational Calculation
Newton's law of universal gravitation is expressed as F = G * m1 * m2 / r². This force is computed for each pair of N celestial bodies:
type CelestialBody = {
x: number
y: number
vx: number
vy: number
mass: number
radius: number
color: string
trail: { x: number; y: number }[]
}
const G = 6.674 // Scaled gravitational constant
function computeGravity(bodies: CelestialBody[]) {
for (let i = 0; i < bodies.length; i++) {
for (let j = i + 1; j < bodies.length; j++) {
const a = bodies[i]
const b = bodies[j]
const dx = b.x - a.x
const dy = b.y - a.y
const distSq = dx * dx + dy * dy
const dist = Math.sqrt(distSq)
// Softening: force diverges when distance is too small
const softening = 10
const force = (G * a.mass * b.mass) / (distSq + softening * softening)
const fx = force * dx / dist
const fy = force * dy / dist
a.vx += fx / a.mass
a.vy += fy / a.mass
b.vx -= fx / b.mass
b.vy -= fy / b.mass
}
}
}Softening is a numerical safety mechanism. When two bodies get extremely close, 1/r² diverges to infinity and velocity explodes. Adding softening² to the denominator limits the force at close range.
Choice of Integration Method
The symplectic Euler from physics.ts would work, but energy conservation matters for orbital mechanics. Velocity Verlet is used:
function integrateVerlet(bodies: CelestialBody[], dt: number) {
// Half-step velocity update (using acceleration from previous step)
for (const b of bodies) {
b.vx += b.ax * dt * 0.5
b.vy += b.ay * dt * 0.5
}
// Position update
for (const b of bodies) {
b.x += b.vx * dt
b.y += b.vy * dt
}
// Compute new acceleration
computeGravity(bodies)
// Remaining half-step velocity update
for (const b of bodies) {
b.vx += b.ax * dt * 0.5
b.vy += b.ay * dt * 0.5
}
}Velocity Verlet is symplectic (preserving phase space volume), so orbits don't drift during long-running simulations. The planets in SolarSystem tracing closed orbits is thanks to this property.
Initial Conditions for SolarSystem
The key to a solar system simulation is the initial conditions. Real planetary data is scaled and used:
const SUN: CelestialBody = {
x: 400, y: 300,
vx: 0, vy: 0,
mass: 10000,
radius: 20,
color: '#FFD700',
trail: [],
}
// Circular orbit initial velocity v = sqrt(GM/r)
function orbitalVelocity(centralMass: number, distance: number): number {
return Math.sqrt(G * centralMass / distance)
}
const EARTH: CelestialBody = {
x: 400 + 150, y: 300,
vx: 0, vy: orbitalVelocity(10000, 150),
mass: 10,
radius: 6,
color: '#4169E1',
trail: [],
}When the initial velocity computed by orbitalVelocity is applied tangentially, the body traces a stable elliptical orbit. Kepler's laws are reproduced within the code.
Orbit Rendering
The Canvas API draws each body's past positions as a trail:
function drawTrail(ctx: CanvasRenderingContext2D, trail: { x: number; y: number }[]) {
if (trail.length < 2) return
ctx.beginPath()
ctx.moveTo(trail[0].x, trail[0].y)
for (let i = 1; i < trail.length; i++) {
ctx.lineTo(trail[i].x, trail[i].y)
}
ctx.strokeStyle = 'rgba(255, 255, 255, 0.3)'
ctx.lineWidth = 1
ctx.stroke()
}
function updateTrail(body: CelestialBody, maxLength: number) {
body.trail.push({ x: body.x, y: body.y })
if (body.trail.length > maxLength) {
body.trail.shift()
}
}Limiting trail length saves memory while retaining enough to see the orbital shape.
GravityWell Interaction
GravityWell places a gravity source at the user's pointer position. Particles are visibly attracted toward the mouse:
function attractToPointer(
particles: CelestialBody[],
pointerX: number,
pointerY: number,
strength: number,
) {
for (const p of particles) {
const dx = pointerX - p.x
const dy = pointerY - p.y
const distSq = dx * dx + dy * dy
const dist = Math.sqrt(distSq)
if (dist < 1) continue
const force = strength / (distSq + 100)
p.vx += (force * dx) / dist
p.vy += (force * dy) / dist
}
}The difference between NBody and GravityWell is whether bodies attract each other. NBody involves mutual interaction between all pairs; GravityWell is unidirectional attraction toward a single fixed gravity source.
The O(n²) Limitation
The computational complexity of the N-body problem is O(n²). With 100 bodies, that's 4,950 force calculations. Scaling further would require the Barnes-Hut algorithm (using a quadtree to group distant bodies), but for a browser demo, 100 bodies is more than enough for beautiful results, so brute-force all-pairs computation was the pragmatic choice. Updating every frame with requestAnimationFrame runs smoothly at this scale.
Choosing a sufficient implementation is the balance between performance and code complexity.
Summary: Tools for N-Body Simulation
The fundamentals of celestial mechanics and vector operations are well covered by the linear algebra and physics simulation books on the tool shelf. As an exercise in translating equations to code, the N-body problem is an ideal subject.
FAQ
- Q. Why do N-body simulations become unstable so easily?
- Gravity grows very rapidly at short distances, so coarse timesteps can cause energy to explode. Stability improves with fixed timesteps, softening/clamping near-zero distances, and more stable integrators (e.g., symplectic Euler).
- Q. How far can you go with O(n²) gravity in the browser?
- Pairwise force computation scales poorly as n grows. For 60fps, keep n modest and optimize vector math first; if you need more, consider approximate methods like Barnes–Hut (O(n log n)).
- Q. Why does the choice of integrator change the orbit so much?
- Numerical error accumulates as artificial energy gain/loss, making orbits spiral in or out. Even at similar cost, more stable update schemes reduce drift and produce more believable trajectories.
- Q. How do you implement a pointer-controlled gravity source safely?
- Convert pointer position to world coordinates, compute a normalized direction vector per particle, and add acceleration toward the source. Add a minimum distance or a softening term to avoid singularities near zero distance.