- Published Results Site
- Generated Analysis
- Comparison Examples
- Two-Body Sun-Earth System
- Three-Body Systems
- N-Body Systems
- Conclusions
- References
This file is the canonical results narrative. GitHub Pages renders the same
docs/RESULTS.md content as the public landing page, with the artifact browser and
interactive method dashboard linked from that page.
https://sidrichardsquantum.github.io/Celestial_Dynamics_Iteration_Methods/
The same static files are generated locally under analysis/generated/ and
images/.
This section is generated by Rscript analysis/generate_results.R.
Do not edit the numeric tables by hand; regenerate them from the solver code.
| method | years | steps | final_separation_au | final_energy_ratio | max_abs_energy_error | max_abs_angular_momentum_drift |
|---|---|---|---|---|---|---|
| Euler | 25 | 1000 | 1.02956e+01 | 1.45946e-01 | 8.54054e-01 | 1.09354e+00 |
| Midpoint | 25 | 1000 | 1.05514e+00 | 9.42958e-01 | 5.70418e-02 | 2.97638e-02 |
| Heun | 25 | 1000 | 1.21910e+00 | 8.49014e-01 | 1.50986e-01 | 8.44417e-02 |
| RK4 | 25 | 1000 | 9.99575e-01 | 1.00042e+00 | 4.21667e-04 | 2.10766e-04 |
| Verlet | 25 | 1000 | 1.00018e+00 | 9.99998e-01 | 1.38278e-04 | 2.18102e-15 |
Errors are measured against a one-year RK4 reference run with N = 16000.
| method | years | steps | dt_days | final_position_error_au | observed_order |
|---|---|---|---|---|---|
| Euler | 1 | 250 | 1.46100e+00 | 1.20694e+00 | |
| Euler | 1 | 500 | 7.30500e-01 | 6.76416e-01 | 8.354e-01 |
| Euler | 1 | 1000 | 3.65250e-01 | 3.58832e-01 | 9.146e-01 |
| Euler | 1 | 2000 | 1.82625e-01 | 1.84827e-01 | 9.571e-01 |
| Midpoint | 1 | 250 | 1.46100e+00 | 2.37405e-03 | |
| Midpoint | 1 | 500 | 7.30500e-01 | 5.87089e-04 | 2.016e+00 |
| Midpoint | 1 | 1000 | 3.65250e-01 | 1.45910e-04 | 2.009e+00 |
| Midpoint | 1 | 2000 | 1.82625e-01 | 3.63660e-05 | 2.004e+00 |
| Heun | 1 | 250 | 1.46100e+00 | 5.51744e-03 | |
| Heun | 1 | 500 | 7.30500e-01 | 1.35203e-03 | 2.029e+00 |
| Heun | 1 | 1000 | 3.65250e-01 | 3.34451e-04 | 2.015e+00 |
| Heun | 1 | 2000 | 1.82625e-01 | 8.31604e-05 | 2.008e+00 |
| RK4 | 1 | 250 | 1.46100e+00 | 6.58641e-08 | |
| RK4 | 1 | 500 | 7.30500e-01 | 3.85816e-09 | 4.094e+00 |
| RK4 | 1 | 1000 | 3.65250e-01 | 2.33057e-10 | 4.049e+00 |
| RK4 | 1 | 2000 | 1.82625e-01 | 1.43231e-11 | 4.024e+00 |
| Verlet | 1 | 250 | 1.46100e+00 | 1.32238e-03 | |
| Verlet | 1 | 500 | 7.30500e-01 | 3.30656e-04 | 2.000e+00 |
| Verlet | 1 | 1000 | 3.65250e-01 | 8.26677e-05 | 2.000e+00 |
| Verlet | 1 | 2000 | 1.82625e-01 | 2.06671e-05 | 2.000e+00 |
The same two-body methods are run for ten lunar months with N = 1000.
| system | method | duration | steps | final_separation_au | final_energy_ratio | max_abs_energy_error | max_abs_angular_momentum_drift |
|---|---|---|---|---|---|---|---|
| Earth-Moon | Euler | 10 lunar months | 1000 | 6.23876e-03 | 3.21768e-01 | 6.78232e-01 | 6.94754e-01 |
| Earth-Moon | Midpoint | 10 lunar months | 1000 | 2.57415e-03 | 9.98012e-01 | 1.98768e-03 | 9.71682e-04 |
| Earth-Moon | Heun | 10 lunar months | 1000 | 2.58973e-03 | 9.92271e-01 | 7.72945e-03 | 3.79312e-03 |
| Earth-Moon | RK4 | 10 lunar months | 1000 | 2.56867e-03 | 1.00000e+00 | 1.81136e-06 | 8.83727e-07 |
| Earth-Moon | Verlet | 10 lunar months | 1000 | 2.56904e-03 | 1.00000e+00 | 1.53235e-05 | 3.99646e-15 |
These RK4 runs check whether known periodic configurations return close to their starting shape after one nominal period.
| case | years | steps | final_energy_ratio | max_abs_energy_error | max_abs_angular_momentum_drift | max_pair_separation_error_au |
|---|---|---|---|---|---|---|
| Figure-8 | 2.05397e+02 | 2e+04 | 1.00000e+00 | 1.93255e-14 | 1.13782e+39 | 1.70987e-08 |
| Lagrange triangle | 3.33146e+02 | 2e+04 | 1.00000e+00 | 4.97494e-15 | 2.48305e-15 | 1.01387e-13 |
| Euler collinear | 1.82472e+02 | 2e+04 | 1.00000e+00 | 7.95990e-15 | 4.14482e-15 | 1.22398e-15 |
| Butterfly I | 2.02434e+02 | 1e+05 | 1.00003e+00 | 2.62380e-05 | 1.87407e+45 | 8.00074e-05 |
The n-body cases report energy and angular-momentum drift for multi-body examples beyond the dedicated two- and three-body solvers.
| case | years | steps | bodies | final_energy_ratio | max_abs_energy_error | max_abs_angular_momentum_drift |
|---|---|---|---|---|---|---|
| Sun-Earth-Mars-Jupiter | 1.20000e+01 | 6000 | 4 | 1e+00 | 1.06181e-11 | 4.55554e-13 |
| Rotating square four-body | 1.23993e+02 | 20000 | 4 | 1e+00 | 5.71739e-15 | 2.81649e-15 |
| Triangular central four-body | 1.62437e+02 | 20000 | 4 | 1e+00 | 2.12602e-14 | 1.06080e-14 |
Runtime is measured for one-year Sun-Earth runs against the same fine RK4 reference used in the convergence table.
| method | years | steps | dt_days | runtime_seconds | final_position_error_au | error_per_second |
|---|---|---|---|---|---|---|
| Euler | 1 | 250 | 1.46100e+00 | 4.0e-03 | 1.20694e+00 | 3.01735e+02 |
| Euler | 1 | 500 | 7.30500e-01 | 5.0e-03 | 6.76416e-01 | 1.35283e+02 |
| Euler | 1 | 1000 | 3.65250e-01 | 1.4e-02 | 3.58832e-01 | 2.56309e+01 |
| Euler | 1 | 2000 | 1.82625e-01 | 4.6e-02 | 1.84827e-01 | 4.01797e+00 |
| Midpoint | 1 | 250 | 1.46100e+00 | 6.0e-03 | 2.37405e-03 | 3.95674e-01 |
| Midpoint | 1 | 500 | 7.30500e-01 | 1.2e-02 | 5.87089e-04 | 4.89241e-02 |
| Midpoint | 1 | 1000 | 3.65250e-01 | 1.5e-02 | 1.45910e-04 | 9.72732e-03 |
| Midpoint | 1 | 2000 | 1.82625e-01 | 5.6e-02 | 3.63660e-05 | 6.49393e-04 |
| Heun | 1 | 250 | 1.46100e+00 | 3.0e-03 | 5.51744e-03 | 1.83915e+00 |
| Heun | 1 | 500 | 7.30500e-01 | 7.0e-03 | 1.35203e-03 | 1.93147e-01 |
| Heun | 1 | 1000 | 3.65250e-01 | 1.8e-02 | 3.34451e-04 | 1.85806e-02 |
| Heun | 1 | 2000 | 1.82625e-01 | 6.0e-02 | 8.31604e-05 | 1.38601e-03 |
| RK4 | 1 | 250 | 1.46100e+00 | 6.0e-03 | 6.58641e-08 | 1.09773e-05 |
| RK4 | 1 | 500 | 7.30500e-01 | 1.6e-02 | 3.85816e-09 | 2.41135e-07 |
| RK4 | 1 | 1000 | 3.65250e-01 | 5.9e-02 | 2.33057e-10 | 3.95011e-09 |
| RK4 | 1 | 2000 | 1.82625e-01 | 8.0e-02 | 1.43231e-11 | 1.79039e-10 |
| Verlet | 1 | 250 | 1.46100e+00 | 3.0e-03 | 1.32238e-03 | 4.40794e-01 |
| Verlet | 1 | 500 | 7.30500e-01 | 7.0e-03 | 3.30656e-04 | 4.72365e-02 |
| Verlet | 1 | 1000 | 3.65250e-01 | 1.4e-02 | 8.26677e-05 | 5.90483e-03 |
| Verlet | 1 | 2000 | 1.82625e-01 | 4.9e-02 | 2.06671e-05 | 4.21779e-04 |
- Energy and angular momentum drift are the main diagnostics for whether an orbit is physically credible over long horizons.
- Euler is useful as a failure baseline: its first-order error produces visible orbital drift even when the trajectory still looks smooth.
- RK4 gives the smallest local error in these examples, while Velocity Verlet is included because symplectic methods often preserve qualitative orbital behavior over long integrations.
- Three-body periodic solutions are validation cases, not generic stability guarantees; the perturbation examples show how quickly nearby trajectories can diverge.
images/analysis/sun_earth_energy_error.pngimages/analysis/sun_earth_angular_momentum_drift.pngimages/analysis/convergence_rates.pnganalysis/generated/earth_moon_method_summary.csvanalysis/generated/three_body_special_summary.csvanalysis/generated/n_body_conservation_summary.csvanalysis/generated/runtime_benchmark.csvanalysis/generated/plot_manifest.csvanalysis/generated/artifact_baseline.csvanalysis/generated/index.htmlanalysis/generated/artifact_index.htmlanalysis/generated/method_comparison_dashboard.htmlimages/two_body/sun_earth/sun_earth_runge_kutta.htmlimages/three_body/special_solutions/three_earths.htmlimages/n_body/sun_earth_mars_jupiter.html
The comparison examples run every implemented two-body method against the same
simple setup. The Sun-Earth comparison uses one year and N = 250 steps so the
Euler method visibly drifts while Midpoint, Heun, RK4, and Velocity Verlet stay
near the expected orbit.
The method list, colors, and Sun-Earth setup are shared through
R/systems/two_body/two_body_method_registry.R.
The energy-error comparison uses an absolute relative error scale so the higher-accuracy methods remain visible next to Euler.
Our Sun and Earth files in the directory examples/two_body/sun_earth/, each plot the orbit of the Earth (blue line) around the Sun (red) over
The Euler Method does not illustrate the full
Below is the code output which includes the energy ratio:
Two-Body System Simulation Euler Method Results:
Body a mass: 1.99e+30 kg
Body b mass: 5.97e+24 kg
Total simulation time: 25.00 years
Time steps: 1000
Time step size: 9.13 days
Initial separation: 1.000 AU
Final separation: 10.296 AU
Energy conservation ratio: 0.145946
Output:
Two-Body System Simulation Midpoint Method Results:
Body a mass: 1.99e+30 kg
Body b mass: 5.97e+24 kg
Total simulation time: 25.00 years
Time steps: 1000
Time step size: 9.13 days
Initial separation: 1.000 AU
Final separation: 1.055 AU
Energy conservation ratio: 0.942958
Heun’s method shows slightly less accuracy than the midpoint method in terms of final orbital radius and energy conservation.
Output:
Two-Body System Simulation Heun's Method Results:
Body a mass: 1.99e+30 kg
Body b mass: 5.97e+24 kg
Total simulation time: 25.00 years
Time steps: 1000
Time step size: 9.13 days
Initial separation: 1.000 AU
Final separation: 1.219 AU
Energy conservation ratio: 0.849014
The most accurate technique by far, as the orbit shows no noticeable perturbations or shifts.
Output:
Two-Body System Simulation Runge-Kutta Method Results:
Body a mass: 1.99e+30 kg
Body b mass: 5.97e+24 kg
Total simulation time: 25.00 years
Time steps: 1000
Time step size: 9.13 days
Initial separation: 1.000 AU
Final separation: 1.000 AU
Energy conservation ratio: 1.000422
Velocity Verlet is a second-order symplectic method. For the Sun-Earth example it keeps the orbit close to circular and preserves energy well over the same 25-year horizon. Its main value is long-term qualitative stability: even when RK4 has a smaller local error, Verlet is often a better baseline for conservative orbital systems.
We used the RK4 method for the three-body systems for accuracy.
The example scripts are grouped under examples/three_body/ by model type:
general/, special_solutions/, perturbations/, and restricted/.
Below illustrates the orbits of both the Earth and Mars around the Sun over a year. As the planets are far away enough, and as the sun is a lot more massive, the orbits appear to be stable until errors in the Runge-Kutta method cause the system to spiral.
Another example illustrates an orbiting spacecraft between the Earth and moon. The spacecraft appears to have an unstable orbit because it is so much lighter than the celestials. Making small changes in the spacecraft's initial position and velocity drastically changes its orbit. For longer periods of time, the spacecraft escapes the Earth–Moon system due to unstable gravitational perturbations.
The hierarchical triple example places two solar-mass stars in a close binary and a Jupiter-mass body on a much wider outer orbit.
Over five years, the binary separation remains close to
A stable solution to the three-body problem exists where three identical massive bodies follow each other in a figure-8 pattern. Each body trails the previous by one-third of the orbital period, forming a periodic figure-8 path "choreography".
The simulation uses standard dimensionless initial conditions for the equal-mass figure-8 orbit and scales them into SI units for three Earth-mass bodies.
For an initial outer-body separation of
With the corrected scaling, RK4 returns the bodies to their starting separations after one period.
The energy output is Energy conservation ratio: 1.000000, so this example validates the periodic figure-8 setup rather than demonstrating immediate destabilization.
The Lagrange solution places three equal masses at the vertices of an equilateral triangle. The whole triangle rotates around its center of mass while preserving all three pairwise separations.
For a side length of Energy conservation ratio: 1.000000.
The Euler collinear solution places three equal masses along one rotating line, with one body at the center of mass and two bodies on opposite sides. This configuration is more delicate than the Lagrange triangle because the bodies remain aligned throughout the orbit.
For an outer-body separation of Energy conservation ratio: 1.000000.
The Butterfly I choreography uses published dimensionless equal-mass initial conditions and scales them to three Earth-mass bodies with an outer separation of Energy conservation ratio: 1.000026.
Small perturbations were added to the figure-8 and Lagrange solutions.
The figure-8 perturbation changes one initial velocity component by
The restricted Earth-Moon example uses normalized CR3BP coordinates. The third body starts near L4, producing a Trojan-style tadpole trajectory in the rotating frame.
The L1 example uses the CR3BP helper and starts near the Earth-Moon L1 point. It is an illustrative planar Lyapunov-like trajectory in the rotating frame, not a corrected mission-grade halo orbit.
The Sitnikov example uses two equal primaries in circular motion and a restricted third body moving perpendicular to their orbital plane. The plot shows the vertical oscillation of the third body over time.
The general n-body engine supports any number of 2D massive bodies using shared array-based position and velocity storage. The first example uses RK4 for a Sun-Earth-Mars-Jupiter system over 12 years. It is intended as a compact demonstration that the repo can now run beyond hard-coded two-body and three-body systems.
The n-body examples also include special four-body central configurations. In these cases the shape is preserved while the bodies rotate rigidly.
- The Euler method lost
$\approx 85%$ for our Sun-Earth simulation - making it the least accurate. - The Runge-Kutta (RK4) method is the most accurate, changing energy by as little as
$\approx 0.04%$ for the same sytem. - Velocity Verlet gives a symplectic comparison method with bounded long-term energy behavior for orbital systems.
- Many three-body systems are unstable, and at least one body can eventually be lost in generic simulations.
- The equal-mass figure-8, Lagrange triangle, Euler collinear, and Butterfly I cases are special periodic solutions; small perturbations can still grow over time.
- Restricted examples are useful for Lagrange points, Trojans, Lyapunov-like motion, and the Sitnikov problem without requiring a fully massive third body.
- Euler Method
- Midpoint Method
- Heun's Method
- Runge-Kutta Method
- Figure-8 Solution A special solution where three equal-mass bodies follow each other in a stable figure-8 path — discovered by Cris Moore (1993) and proven by Montgomery & Chenciner (2000).
- Lagrange and Euler Solutions Classical central-configuration solutions where the triangle or line rotates while retaining its shape.
- Periodic Three-Body Initial Conditions Cataloged equal-mass periodic orbit initial conditions, including Butterfly I.
📘 Author: Sid Richards (SidRichardsQuantum)
This project is licensed under the MIT License - see the LICENSE file for details.