nbody-sim — Gravitational N-Body Simulator

A gravitational N-body simulator written from scratch in Python + NumPy, whose core is a pluggable numerical integrator validated by conservation laws. The thesis of the project: which integrator you use matters more than per-step accuracy — explicit Euler makes orbits spiral out, symplectic integrators conserve energy forever.

---

The signature figure

The same Sun–Earth orbit (eccentric, e ≈ 0.19), the same step dt, ~400 orbits — only the integrator changes. The vertical axis is the relative energy drift $|\Delta E / E_0|$ on a logarithmic scale:

• Explicit Euler (red) grows in a secular, monotonic way up to ~91% — every step injects spurious energy into the orbit, which spirals outward. This is what happens in a naive simulator. It is not an implementation bug: it is a structural property of the method.
• Semi-implicit Euler and Leapfrog (orange/blue) are symplectic: the energy error stays bounded and oscillating around the true value, with no accumulated drift — forever. Leapfrog (2nd order) has a much narrower error band than semi-implicit Euler (1st order).
• RK4 (green) is the most accurate in the short term (4th order), but it is not symplectic: a slow secular drift creeps in and, over very long integrations, it eventually loses to the symplectic schemes.

This is the central argument of the project, made visible. The rest of the code exists to back this figure up with rigor.

The visualization (Pygame, "space-grade")

Inner system (Sun → Mars)	Full system (out to Neptune)

---

What it is / why it's interesting

Most "weekend" gravitational simulators use explicit Euler because it is the obvious method: update position with velocity, update velocity with acceleration. The problem is that the orbits spiral — energy grows without bound and the simulation lies about the physics after a few dozen orbits.

This project treats the time integrator as a first-class object of study. The four schemes (explicit Euler, semi-implicit Euler, Leapfrog/Velocity-Verlet, RK4) live behind a common interface and can be swapped in one line — or live, with a keystroke, inside the visualization. This makes it possible to compare symplectic vs. non-symplectic rather than just assert the difference, and to use the conservation laws (energy, linear momentum, angular momentum, center of mass) as an objective correctness criterion.

The result is validated at two milestones: the analytic Sun–Earth system (Milestone 1) and the real Solar System, with NASA JPL Horizons ephemerides at the J2000 epoch (Milestone 2).

---

Theoretical background

1. Gravitation and the ODE system

The gravitational force on body $i$ due to the others (Newton's law) produces the acceleration

$$ \mathbf{a}_i ;=; \frac{\mathbf{F}_i}{m_i} ;=; G \sum_{j \neq i} m_j , \frac{\mathbf{r}_j - \mathbf{r}_i}{\lVert \mathbf{r}_j - \mathbf{r}_i \rVert^{3}}. $$

Each body is described by position and velocity. Writing the state as $\mathbf{y} = (\mathbf{x}, \mathbf{v})$, the problem becomes a first-order system of ODEs:

$$ \dot{\mathbf{x}} = \mathbf{v}, \qquad \dot{\mathbf{v}} = \mathbf{a}(\mathbf{x}). $$

This is exactly the form the integrators receive: a function accel_fn(pos) -> acc that evaluates $\mathbf{a}(\mathbf{x})$, and nothing more. The integrator never knows what the force law is.

2. AU / year / M☉ nondimensionalization and why G = 4π²

Working in SI carries the measurement uncertainty of $G$ and enormous numbers. We adopt the astronomical system: length in astronomical units (AU), time in years, mass in solar masses (M☉). In this system $G$ stops being a measured constant and becomes an exact number.

By Kepler's third law, a circular orbit of radius $a$ around a mass $M$ has period

$$ T^2 = \frac{4\pi^2}{G,M}, a^3. $$

By definition of our units, the Earth orbits with $a = 1,\text{AU}$, $T = 1,\text{yr}$ and the Sun has $M = 1,\text{M}_\odot$. Substituting:

$$ 1 = \frac{4\pi^2}{G \cdot 1} \cdot 1 ;;\Longrightarrow;; \boxed{G = 4\pi^2}. $$

Practical consequence: the circular orbital velocity at 1 AU is exactly $v_\text{circ} = \sqrt{G M / r} = 2\pi$ AU/yr, and the Earth's orbit closes at $t = 1$ yr by construction. Every constant lives in src/units.py — no magic number leaks into the rest of the code.

3. Softening ε — taming the 1/r² singularity

The Newtonian force $\propto 1/r^2$ diverges when two bodies approach each other ($r \to 0$): a close encounter produces an arbitrarily large acceleration that "blows up" the integrator. The standard solution is Plummer softening: replace the singular denominator with a regularized one,

$$ \mathbf{a}_i ;=; G \sum_{j \neq i} m_j , \frac{\mathbf{r}_j - \mathbf{r}_i}{\bigl(\lVert \mathbf{r}_j - \mathbf{r}_i \rVert^{2} + \varepsilon^{2}\bigr)^{3/2}}. $$

With $\varepsilon > 0$ the force is bounded from above for separations smaller than $\varepsilon$; with $\varepsilon = 0$ exact Newtonian gravity is recovered. The potential used in the diagnostics uses the same $\sqrt{r^2 + \varepsilon^2}$, so that force and energy remain mutually consistent (see gradient check below). The Solar System milestones run with $\varepsilon = 0$ (exact Newtonian); softening exists for dense scenarios.

4. Symplectic vs. non-symplectic integrators

The flow of the N-body Hamiltonian system preserves the symplectic structure of phase space (the Liouville volume, the 2-form $\sum dp \wedge dq$). An integrator is symplectic when its one-step map also preserves this structure. The deep consequence: a symplectic integrator solves exactly a neighboring Hamiltonian (the shadow Hamiltonian) that differs from the true one by $O(\Delta t^p)$. That is why the simulation energy oscillates around the true value with bounded amplitude forever, instead of drifting.

A non-symplectic integrator has no such protection: the energy error accumulates in a secular (monotonic) way, even if each individual step is very accurate.

Scheme	Global order	Symplectic?	Energy behavior
Explicit Euler	$O(\Delta t)$	no	secular drift, spirals outward
Semi-implicit Euler	$O(\Delta t)$	yes	bounded (wide band)
Leapfrog / Velocity-Verlet (default)	$O(\Delta t^2)$	yes	bounded and oscillating (narrow band)
RK4	$O(\Delta t^4)$	no	tiny in the short term, slow secular drift

The difference between explicit and semi-implicit Euler is a single line: the semi-implicit one uses the already updated velocity to move the position. This trivial swap is what makes the scheme symplectic — a perfect textbook example that the structure of the method, not just its order, governs long-term behavior.

Leapfrog in kick-drift-kick form (the default integrator):

$$ \begin{aligned} \mathbf{v}_{n+1/2} &= \mathbf{v}_n + \tfrac{1}{2},\mathbf{a}(\mathbf{x}_n),\Delta t && \text{(half kick)}\\ \mathbf{x}_{n+1} &= \mathbf{x}_n + \mathbf{v}_{n+1/2},\Delta t && \text{(full drift)}\\ \mathbf{v}_{n+1} &= \mathbf{v}_{n+1/2} + \tfrac{1}{2},\mathbf{a}(\mathbf{x}_{n+1}),\Delta t && \text{(half kick)} \end{aligned} $$

The scheme is time-symmetric (reversible) and symplectic; the local error is $O(\Delta t^3)$, the global one $O(\Delta t^2)$. Why does RK4, with a much smaller per-step error, end up losing? Because its drift, although slow, is secular: it grows without bound with the number of periods. Leapfrog's bounded energy beats any monotonic drift over long horizons — exactly what the signature figure shows.

5. Conservation as a correctness test

In an isolated system, linear momentum $\mathbf{P}$, center of mass, and angular momentum $\mathbf{L}$ are conserved by symmetry (translational/rotational invariance of the force law), independently of the discretization error — so any sane integrator must preserve them to machine precision. Energy $E = T + U$ is the finest test, because it separates symplectic from non-symplectic methods. These four invariants are measured at every step in src/diagnostics.py and used as validation asserts in the experiments and in the test suite.

As proof that force and potential are mutually consistent, there is a gradient check: verifying that $\mathbf{a}i = -\tfrac{1}{m_i}\nabla{\mathbf{r}_i} U$ by central finite differences. Maximum measured error < 1e-7 in the Newtonian case ($\varepsilon = 0$) and < 1e-6 with softening ($\varepsilon &gt; 0$) — confirming that force and potential use the same $\sqrt{r^2 + \varepsilon^2}$ in the regularized regime too.

---

Architecture

Structure-of-Arrays (SoA) layout, ready for 3D ($z = 0$ in the plane for planar problems): pos[N,3] (AU), vel[N,3] (AU/yr), mass[N] (M☉). SoA keeps the arrays contiguous and lets all the physics be expressed via NumPy broadcasting, with no Python loops over particle pairs.

Module	Responsibility
src/units.py	AU/year/M☉ unit system, $G = 4\pi^2$ exact, planetary masses (JPL DE440), SI conversions. Single source of truth — no magic numbers in the rest of the code.
src/forces.py	Fully vectorized $O(N^2)$ gravitational accelerations (difference tensor [N,N,3] by broadcasting) with Plummer softening.
src/integrators.py	Pluggable integrator behind a common interface (Integrator.step). Explicit Euler, semi-implicit Euler, Leapfrog (default), RK4. get_integrator(name) selects by name.
src/diagnostics.py	Conserved quantities: energy $T+U$, linear momentum $\mathbf{P}$, center of mass, angular momentum $\mathbf{L}$. These are the correctness criteria.
src/simulation.py	Simulation engine (simulate) + State and History. Decoupled from any rendering — it only knows about state, integrator and force law.
src/initial_conditions.py	Scenario builders (e.g. sun_earth), with momentum balancing and recentering on the CM.
src/data_loader.py	Real Solar System ephemerides via JPL Horizons (astroquery), converted to internal units, .npz cache + built-in offline fallback.
src/render.py	2D Pygame visualization. 100% decoupled from the physics — every drawn position comes from the existing engine, it never reimplements gravity.

The central design is the pluggable integrator. The interface is deliberately minimal:

Integrator.step(pos, vel, mass, dt, accel_fn) -> (pos_new, vel_new)

where accel_fn(pos) -> acc[N,3] is a callback. The integrator does not know the force law, which lets you swap the scheme without touching the physics and drives the comparative experiment.

---

How to run

Installation

pip install -r requirements.txt

NumPy and Matplotlib are required. astroquery (Horizons), pygame (visualization), pillow (GIF export), hypothesis (property-based tests) are optional — the code degrades gracefully when they are absent (offline fallback, clear messages, skips).

Run everything from the repository root. The scripts also accept direct execution (python experiments/file.py), but python -m experiments.file is the recommended path.

Milestone 1 — Sun–Earth (Leapfrog, 10 years)

python -m experiments.marco1_two_body

Integrates the circular Sun–Earth orbit with Leapfrog, prints the validation numbers (measured period, Leapfrog vs. Euler energy drift, conservation of P/CM/L) and saves experiments/energia_marco1.png.

Milestone 2 — Real Solar System (JPL Horizons, 100 years)

python -m experiments.marco2_solar_system

Loads Sun + 8 planets at the J2000 epoch via Horizons, converts to AU/yr/M☉ and integrates 100 years with Leapfrog. Notable points:

• Cache: the 1st run queries Horizons and writes data/cache/solar_system_J2000.npz; subsequent runs load from the cache and do not touch the network (use force_refresh=True to re-query).
• Offline fallback: without astroquery/network, it uses approximate embedded J2000 vectors (flagged on stdout; not ephemeris-grade) — the simulation runs anywhere.
• Saves experiments/marco2_solar_system.png (orbits + energy drift) and the trajectory history data/cache/marco2_trajectory.npz for the renderer.

Signature experiment — energy drift of the 4 integrators

python -m experiments.energy_drift_comparison

Integrates the same eccentric Sun–Earth orbit for ~400 orbits with the four schemes and plots $|\Delta E / E_0|(t)$ on log-y. Generates the figure at the top, experiments/energy_drift_comparison.png.

Visualization (Pygame)

python -m experiments.visualize                  # live integration (default)
python -m experiments.visualize --mode replay    # replay the cached trajectory
python -m experiments.visualize --view full      # full system (out to Neptune)
python -m experiments.visualize --integrator rk4 # start with another integrator

In live mode, one engine step per (sub)frame starting from the data_loader state; switching the integrator with I lets you watch symplectic vs. non-symplectic stability diverge on the live energy readout in the HUD. Replay mode plays back data/cache/marco2_trajectory.npz (falls back to live if the cache does not exist).

Controls:

Key	Action	Key	Action
SPACE	pause / resume	T	toggle trails
+ / -	speed up / slow down time	L	toggle labels
I	cycle integrator (euler→semi→leapfrog→rk4)	C	clear trails
V	toggle inner / full view	R	reset simulation
scroll	zoom	drag	pan
H	controls overlay	ESC / Q	quit

Headless capture (no display — used to generate the README images):

python -m experiments.visualize --headless --frames 360 --substeps 40 \
    --view full --screenshot assets/screenshot_full.png --gif assets/orbits.gif

Tests

python -m pytest tests/ -q

35 tests. They cover: the core gradient check (force = $-\nabla U$/m, with and without softening), conservation parametrized across the 4 integrators, the two-body problem, golden/regression by hash, property-based tests with Hypothesis, and the offline loader. (test_smoke.py also runs as a script: python tests/test_smoke.py.)

---

Validation / Results

All the numbers below are measured by the scripts in this repository.

Core gradient check

$\mathbf{a} = -\tfrac{1}{m}\nabla U$ by central finite differences — maximum error < 1e-7. Proves that forces.py and diagnostics.py are mutually consistent.

Milestone 1 — Sun–Earth, Leapfrog, 10 years

Metric	Measured value
Orbital period	1.000021 yr (error 0.002%)
Max energy drift $	\Delta E/E_0
Linear momentum, CM, angular momentum	conserved to machine precision

Milestone 2 — Real Solar System (Horizons, 100 years)

Planet	Measured period	Real period	Error
Mercury	0.241	0.2408	0.05%
Venus	0.615	0.6152	0.002%
Earth	1.000	1.0000	0.005%
Mars	1.880	1.8808	0.04%

Conservation: max $|\Delta E/E|$ 3.5e-7 (bounded, no secular drift); angular momentum and barycenter to machine precision.

Signature experiment — energy drift (~400 orbits)

| Integrator | Order | Symplectic | final $|\Delta E/E_0|$ | |---|:---:|:---:|---:| | Explicit Euler | 1 | no | 0.91 (spirals outward) | | Semi-implicit Euler | 1 | yes | 2.5e-3 (bounded) | | Leapfrog / Verlet | 2 | yes | 7.5e-6 (bounded) | | RK4 | 4 | no | 2.9e-7 (best in the short term; slow secular drift) |

Test suite

35 passed — parametrized conservation (4 integrators), two-body, golden by hash, property-based (Hypothesis), offline loader.

---

Roadmap / next steps

The core is validated; the next milestones are not yet implemented and define the direction of the project:

• Milestone 4 — Interactive sandbox: gravity assist scenarios, Lagrange points, live editing of initial conditions.
• Milestone 5 — Barnes-Hut $O(N \log N)$: octree to scale from dozens to thousands of bodies, beating the $O(N^2)$ bottleneck of the direct sum.
• Milestone 6 — Post-Newtonian relativistic correction: 1st-order PN term to reproduce Mercury's perihelion precession (43"/century) — the classic test that pure Newtonian gravity does not capture.

Known limitations and technical debt

• The force sum is $O(N^2)$ and materializes the [N,N,3] tensor: great for small N (the current milestones), but infeasible for thousands of bodies until Milestone 5.
• Leapfrog does two force evaluations per step; reusing the second as the first of the next step was left out for clarity (pending micro-optimization).
• Fixed time step; there is no adaptive error control nor per-body individual timestep (relevant for systems with very disparate time scales).
• The data_loader offline fallback is approximate (not ephemeris-grade): good for running on any machine, but ephemeris-grade numbers require the Horizons query.

---

Credits / license

• Ephemeris data: NASA JPL Horizons (via astroquery), J2000 epoch.
• Planetary masses: JPL planetary GM (DE440) / GM☉ — physical constants, not derived from the Horizons state query (which returns kinematics only).
• Licensed under the MIT License.

Portfolio project — Eduardo, Computer Engineering.