What is the N-body problem?
The N-body problem asks how N point masses move under their mutual gravitational attraction. Every body pulls on every other according to Newton’s law of universal gravitation, F = G·m₁m₂ / r², so the acceleration of any body is the vector sum of the pull from all the others. For N = 2 the problem has a closed-form solution (Kepler’s orbits). For N ≥ 3 there is no general closed-form solution — Poincaré proved the motion can be chaotic — so the equations of motion are integrated numerically, which is exactly what this simulator does.
Why it matters in celestial mechanics
Almost every real gravitational system is an N-body system: the Solar System, globular clusters, galaxies, and the rings and moons of the giant planets. Numerical N-body integration is used to compute planetary ephemerides, to test the long-term stability of orbits, to model galaxy mergers, and to design spacecraft trajectories that exploit gravity assists. The three-body case underlies Lagrange points, horseshoe and tadpole orbits, and the sensitive dependence on initial conditions that makes long-term prediction fundamentally limited.
How to use the simulator
Click and drag on the canvas to launch a body: the drag direction and length set the initial velocity vector v, and the mass slider sets m. The presets load a simplified solar system, a binary star with a circumbinary planet, and an unstable three-body configuration. The readout reports the body count, the current G, and the kinetic, potential and total energy (Eₖ, Eₚ, Eₜₒₜ); the drift figure shows how far Eₜₒₜ has moved from its initial value — a direct diagnostic of integration accuracy. Scroll to zoom, switch to navigate mode to pan, and toggle trails to trace the orbital paths.
Note: G is rescaled to canvas units (pixels, seconds, arbitrary mass) so orbits stay visible on screen; it is not the astronomical value 6.674×10⁻¹¹ N·m²·kg⁻². The dynamics are physically faithful: pairwise Newtonian forces, fixed-step Velocity Verlet integration, and momentum-conserving inelastic collisions.
Frequently asked questions
What numerical integration method does this simulator use?
Velocity Verlet — a second-order symplectic integrator. It advances position using the current acceleration, recomputes the forces, then updates velocity with the average of the old and new accelerations. Because it is symplectic it conserves a nearby 'shadow' Hamiltonian, so total energy oscillates within a bounded band instead of drifting away as it would under naive forward Euler or even a standard Runge-Kutta scheme.
Why is a two-body orbit always an ellipse?
For an inverse-square central force, solving the equation of motion in polar coordinates yields a conic section with the central mass at one focus. Bound orbits (negative total energy) are ellipses — Kepler's first law — with the circle as the zero-eccentricity special case; parabolas and hyperbolas correspond to zero and positive total energy. Closed bound orbits are a special property of the 1/r^2 force (shared only with the harmonic potential), by Bertrand's theorem.
Why does the total energy drift slightly?
Any fixed-step integrator makes a small truncation error each step. Velocity Verlet keeps that error bounded and oscillatory rather than cumulative, so the drift readout usually stays within a fraction of a percent. Drift grows during very close encounters, where accelerations spike and the fixed step under-resolves the motion; the Plummer softening term limits those spikes so the step stays adequate.
What is the softening parameter?
It is a small length epsilon added inside the force law, so the effective force is G*m1*m2 / (r^2 + epsilon^2)^(3/2) — Plummer softening. It caps the otherwise singular force when two bodies get very close, avoiding the numerical blow-ups a fixed time step cannot resolve. The trade-off is that gravity is slightly underestimated at separations comparable to epsilon.