A 2D gravitational simulation of the Solar System that quantifies why integrator choice matters. Ten bodies, mutual gravity, no shortcuts — the Sun, eight planets and a satellite on an Earth-to-Mars transfer, integrated with two different numerical schemes so they can be measured against each other.
The headline result: the Beeman integrator conserves energy 4,781× better than explicit Euler over the same five simulated years, at essentially the same computational cost.
Total system energy is a conserved quantity. Any drift is pure numerical error, which makes it the cleanest available measure of integrator quality.
| Integrator | Order | Peak energy drift over 5 years |
|---|---|---|
| Explicit Euler | 1st | 1.26 % |
| Beeman | 3rd | 0.00026 % |
Euler bleeds energy monotonically — the system is visibly and steadily gaining energy, and the error compounds. Beeman oscillates around the true value without accumulating, which is what you want from an integrator you intend to run for a long time.
Nothing in the code is told how long a year is. Orbital periods emerge from the force calculation and are measured by detecting when each body crosses the positive x-axis.
| Planet | Published (yr) | Simulated (yr) | Error |
|---|---|---|---|
| Mercury | 0.2409 | 0.2411 | 0.06 % |
| Venus | 0.6152 | 0.6152 | 0.01 % |
| Earth | 0.9999 | 1.0005 | 0.06 % |
| Mars | 1.8810 | 1.8824 | 0.08 % |
All four inner planets land within 0.1 % of their published values. The four outer planets are not scored: the run covers 5.01 years and Jupiter's year is 11.86, so none of them completes an orbit.
The satellite waits for a launch window — detected by tracking the Earth–Sun–Mars phase angle until it falls between 44° and 46° — then executes a burn along Earth's velocity vector. The transfer takes 1.53 years.
The trajectory is the white dashed ellipse in the figure at the top. It reaches an aphelion of 2.06 AU, beyond Mars' 1.52 AU, so it crosses Mars' orbit rather than meeting it tangentially. See Known issues for what that means and why the closest-approach figure needs re-deriving.
git clone https://github.com/Anay2005/solar-system.git
cd solar-system
pip install -e .
solar-systemAn animation window opens showing the orbits; the orbital period table prints to
the terminal, and energy data is written to output/.
To regenerate the figures in this README:
python scripts/make_figures.pyRequires Python 3.9+, NumPy, Matplotlib and tabulate.
Force calculation. Every timestep, the gravitational force on each body is summed over every other body — an O(N²) direct sum. No Barnes-Hut approximation, no softening: at ten bodies the exact calculation is affordable and avoids introducing another source of error into an experiment that is specifically about numerical error.
Beeman integration. A predictor-corrector scheme that uses acceleration at three timesteps:
x(t+Δt) = x(t) + v(t)Δt + [4a(t) − a(t−Δt)] Δt² / 6
v(t+Δt) = v(t) + [2a(t+Δt) + 5a(t) − a(t−Δt)] Δt / 6
Position is third-order accurate. The velocity update needs a(t+Δt), so forces
are evaluated after the position update — the extra bookkeeping is what buys the
energy stability shown above.
Euler integration, for comparison, uses only the state at the start of the step:
x(t+Δt) = x(t) + v(t)Δt
v(t+Δt) = v(t) + a(t)Δt
First-order, and it shows.
Why the comparison is fair. Simulate_Euler inherits the force calculation
and energy accounting from Simulate_beeman, so the only difference between the
two runs is the integration step itself. Two separate force implementations would
have confounded the result.
Launch window detection. The Earth–Sun–Mars angle is recovered from the three mutual distances using the cosine rule, and the burn fires when it enters the 44–46° band. The burn magnitude comes from the standard Hohmann transfer expression:
Δv₁ = √[GM(2/r_earth − 1/a_transfer)] − √(GM/r_earth), a_transfer = (r_earth + r_mars)/2
solar-system/
├── src/solarsystem/
│ ├── constants.py Physical constants and published reference periods
│ ├── paths.py Data and output file locations
│ ├── data.py JSON loading; initial positions and velocities
│ ├── physics.py Hohmann transfer parameters
│ ├── plotting.py Energy-vs-time figure
│ ├── cli.py Entry point
│ └── integrators/
│ ├── beeman.py Beeman integrator, energy and period analysis
│ └── euler.py Explicit Euler, for comparison
├── data/planets.json Masses, orbital radii, timestep, iteration count
├── scripts/make_figures.py Regenerates the figures in this README
├── tests/ Smoke tests
├── docs/ Full write-up (PDF)
└── FINDINGS.md Self-review: known issues, with evidence
Configuration lives entirely in data/planets.json — masses, orbital radii,
timestep and iteration count. Adding a body means adding a JSON entry; no code
change required.
pytest # run the test suiteI audited this code after writing it and found several real problems. They are documented in full, with measurements and suggested fixes, in FINDINGS.md. The three that matter most:
-
The satellite's state arrays are off by one relative to every other body after launch, because the launch routine appends a state the integrator had already appended. The reported closest approach to Mars (167,963 km) is therefore computed between the satellite and Mars one timestep apart. Correcting the alignment gives roughly 632,012 km — which is outside Mars' sphere of influence, so this changes the conclusion, not just the number.
-
Before launch the satellite sits at the origin, on top of the Sun, and the 1/r² force produces a numerical blow-up. Its state is overwritten at launch so the post-launch trajectory is sound, but the pre-launch portion of the array is meaningless.
-
The manoeuvre is not really a Hohmann transfer. An undocumented
+ (1000, 1000)term in the launch velocity brings the applied delta-v to 5,082 m/s against the textbook Hohmann figure of 3,079 m/s — 1.65× the fuel, and an aphelion that overshoots Mars.
The energy comparison and the orbital period results are unaffected by all three and stand as reported.
Written as a computational physics project. The full write-up — derivations, method, error analysis — is in docs/Solar_system_simulation.pdf.
MIT — see LICENSE.