Numerical methods and orbital simulation

Celestial Dynamics Results

Generated results for projectile, two-body, three-body, restricted three-body, and n-body examples.

Methods5
Systems4
Views3

Table of Contents

Published Results Site

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/.

Generated Analysis

This section is generated by Rscript analysis/generate_results.R. Do not edit the numeric tables by hand; regenerate them from the solver code.

Sun-Earth Method Summary

methodyearsstepsfinal_separation_aufinal_energy_ratiomax_abs_energy_errormax_abs_angular_momentum_drift
Euler2510001.02956e+011.45946e-018.54054e-011.09354e+00
Midpoint2510001.05514e+009.42958e-015.70418e-022.97638e-02
Heun2510001.21910e+008.49014e-011.50986e-018.44417e-02
RK42510009.99575e-011.00042e+004.21667e-042.10766e-04
Verlet2510001.00018e+009.99998e-011.38278e-042.18102e-15

Sun-Earth Convergence Summary

Errors are measured against a one-year RK4 reference run with N = 16000.

methodyearsstepsdt_daysfinal_position_error_auobserved_order
Euler12501.46100e+001.20694e+00
Euler15007.30500e-016.76416e-018.354e-01
Euler110003.65250e-013.58832e-019.146e-01
Euler120001.82625e-011.84827e-019.571e-01
Midpoint12501.46100e+002.37405e-03
Midpoint15007.30500e-015.87089e-042.016e+00
Midpoint110003.65250e-011.45910e-042.009e+00
Midpoint120001.82625e-013.63660e-052.004e+00
Heun12501.46100e+005.51744e-03
Heun15007.30500e-011.35203e-032.029e+00
Heun110003.65250e-013.34451e-042.015e+00
Heun120001.82625e-018.31604e-052.008e+00
RK412501.46100e+006.58641e-08
RK415007.30500e-013.85816e-094.094e+00
RK4110003.65250e-012.33057e-104.049e+00
RK4120001.82625e-011.43231e-114.024e+00
Verlet12501.46100e+001.32238e-03
Verlet15007.30500e-013.30656e-042.000e+00
Verlet110003.65250e-018.26677e-052.000e+00
Verlet120001.82625e-012.06671e-052.000e+00

Earth-Moon Method Summary

The same two-body methods are run for ten lunar months with N = 1000.

systemmethoddurationstepsfinal_separation_aufinal_energy_ratiomax_abs_energy_errormax_abs_angular_momentum_drift
Earth-MoonEuler10 lunar months10006.23876e-033.21768e-016.78232e-016.94754e-01
Earth-MoonMidpoint10 lunar months10002.57415e-039.98012e-011.98768e-039.71682e-04
Earth-MoonHeun10 lunar months10002.58973e-039.92271e-017.72945e-033.79312e-03
Earth-MoonRK410 lunar months10002.56867e-031.00000e+001.81136e-068.83727e-07
Earth-MoonVerlet10 lunar months10002.56904e-031.00000e+001.53235e-053.99646e-15

Special Three-Body Conservation Summary

These RK4 runs check whether known periodic configurations return close to their starting shape after one nominal period.

caseyearsstepsfinal_energy_ratiomax_abs_energy_errormax_abs_angular_momentum_driftmax_pair_separation_error_au
Figure-82.05397e+022e+041.00000e+001.93255e-141.13782e+391.70987e-08
Lagrange triangle3.33146e+022e+041.00000e+004.97494e-152.48305e-151.01387e-13
Euler collinear1.82472e+022e+041.00000e+007.95990e-154.14482e-151.22398e-15
Butterfly I2.02434e+021e+051.00003e+002.62380e-051.87407e+458.00074e-05

N-Body Conservation Summary

The n-body cases report energy and angular-momentum drift for multi-body examples beyond the dedicated two- and three-body solvers.

caseyearsstepsbodiesfinal_energy_ratiomax_abs_energy_errormax_abs_angular_momentum_drift
Sun-Earth-Mars-Jupiter1.20000e+01600041e+001.06181e-114.55554e-13
Rotating square four-body1.23993e+022000041e+005.71739e-152.81649e-15
Triangular central four-body1.62437e+022000041e+002.12602e-141.06080e-14

Runtime and Accuracy Benchmark

Runtime is measured for one-year Sun-Earth runs against the same fine RK4 reference used in the convergence table.

methodyearsstepsdt_daysruntime_secondsfinal_position_error_auerror_per_second
Euler12501.46100e+002.0e-031.20694e+006.03469e+02
Euler15007.30500e-015.0e-036.76416e-011.35283e+02
Euler110003.65250e-011.2e-023.58832e-012.99027e+01
Euler120001.82625e-013.6e-021.84827e-015.13407e+00
Midpoint12501.46100e+003.0e-032.37405e-037.91348e-01
Midpoint15007.30500e-017.0e-035.87089e-048.38699e-02
Midpoint110003.65250e-011.5e-021.45910e-049.72732e-03
Midpoint120001.82625e-014.3e-023.63660e-058.45722e-04
Heun12501.46100e+003.0e-035.51744e-031.83915e+00
Heun15007.30500e-016.0e-031.35203e-032.25338e-01
Heun110003.65250e-011.5e-023.34451e-042.22967e-02
Heun120001.82625e-014.4e-028.31604e-051.89001e-03
RK412501.46100e+006.0e-036.58641e-081.09773e-05
RK415007.30500e-011.1e-023.85816e-093.50742e-07
RK4110003.65250e-012.4e-022.33057e-109.71069e-09
RK4120001.82625e-016.4e-021.43231e-112.23799e-10
Verlet12501.46100e+003.0e-031.32238e-034.40794e-01
Verlet15007.30500e-016.0e-033.30656e-045.51093e-02
Verlet110003.65250e-011.4e-028.26677e-055.90483e-03
Verlet120001.82625e-014.2e-022.06671e-054.92075e-04

Interpretation Notes

Generated Figures

Sun-Earth energy error
Sun-Earth energy error
Sun-Earth angular momentum drift
Sun-Earth angular momentum drift
Convergence rates
Convergence rates

Comparison Examples

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.

Sun-Earth all methods
Sun-Earth all methods

The energy-error comparison uses an absolute relative error scale so the higher-accuracy methods remain visible next to Euler.

Sun-Earth all-method energy error
Sun-Earth all-method energy error

Two-Body Sun-Earth System

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 $25$ years, using a different iteration method. The simulation initializes Earth $1 AU$ along the x-axis from the Sun at the origin. We use a high number of steps, $N=1000$, to maximize integration accuracy. $25$ years is chosen because it spans multiple orbital periods and exposes long-term integration errors.

Euler Method

The Euler Method does not illustrate the full $25$ periods. The Earth also spirals out immediately, conserving very little total mechanical energy - leading to rapid orbital divergence. Therefore, the Euler method is not an accurate technique, so methods with a higher-order or greater resolution are required to trace-out celestial dynamics better. Euler fails due to accumulating first-order errors that destabilize the system rapidly, especially in energy-conserving dynamics like orbital motion.

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
Euler Method, Sun and Earth system
Euler Method, Sun and Earth system

Midpoint Method

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
Midpoint Method, Sun and Earth system
Midpoint Method, Sun and Earth system

Heun's Method

Heun’s method shows slightly less accuracy than the midpoint method in terms of final orbital radius and energy conservation. $\approx 15%$ of total energy is lost after $25$ years, compared to $\approx 6%$ for the midpoint method.

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
Heun's Method, Sun and Earth system
Heun's Method, Sun and Earth system

Runge-Kutta Method

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
Runge-Kutta Method, Sun and Earth system
Runge-Kutta Method, Sun and Earth system

Velocity Verlet Method

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.

Velocity Verlet Method, Sun and Earth system
Velocity Verlet Method, Sun and Earth system

Three-Body 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/.

Sun, Earth and Mars

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.

Sun, Earth and Mars
Sun, Earth and Mars

Earth, Moon and Spacecraft

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.

Earth, Moon and Spacecraft
Earth, Moon and Spacecraft

Binary with a Distant Third Body

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 $0.2 AU$ while the third body traces a wider path around the pair.

Binary with Distant Third Body
Binary with Distant Third Body

Figure-8 Solution

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 $1 AU$, one full figure-8 period is about $205.4$ years.

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.

Three Earths
Three Earths

Lagrange Equilateral Solution

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 $1 AU$ between three Earth-mass bodies, one full period is about $333.2$ years. RK4 returns each pairwise separation to $1 AU$ after one period, with output ``Energy conservation ratio: 1.000000``.

Lagrange Three Earths
Lagrange Three Earths

Euler Collinear Solution

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 $1 AU$, one full period is about $182.5$ years. RK4 returns the adjacent separations to $0.5 AU$ and the outer separation to $1 AU$ after one period, with output ``Energy conservation ratio: 1.000000``.

Euler Collinear Three Earths
Euler Collinear Three Earths

Butterfly I Choreography

The Butterfly I choreography uses published dimensionless equal-mass initial conditions and scales them to three Earth-mass bodies with an outer separation of $1 AU$. It needs a finer step count than the figure-8 example; with $N=100000$, RK4 returns the starting separations after one period with output ``Energy conservation ratio: 1.000026``.

Butterfly I Choreography
Butterfly I Choreography

Perturbed Special Solutions

Small perturbations were added to the figure-8 and Lagrange solutions. The figure-8 perturbation changes one initial velocity component by $2%$ and is run over three nominal periods. The Lagrange perturbation shifts one body by $0.002 AU$ and produces much larger shape changes over two periods.

Perturbed Figure-8
Perturbed Figure-8
Perturbed Lagrange Three Earths
Perturbed Lagrange Three Earths

Restricted Earth-Moon Trojan

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.

Restricted Earth-Moon Trojan
Restricted Earth-Moon Trojan

Planar Lyapunov-Like L1 Orbit

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.

Lyapunov-Like Orbit Near L1
Lyapunov-Like Orbit Near L1

Sitnikov Problem

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.

Sitnikov Three-Body Problem
Sitnikov Three-Body Problem

N-Body Systems

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.

Sun-Earth-Mars-Jupiter
Sun-Earth-Mars-Jupiter

The n-body examples also include special four-body central configurations. In these cases the shape is preserved while the bodies rotate rigidly.

Rotating square four-body solution
Rotating square four-body solution
Triangular central four-body solution
Triangular central four-body solution

Conclusions:

References

---

📘 Author: Sid Richards (SidRichardsQuantum)

LinkedIn: https://www.linkedin.com/in/sid-richards-21374b30b/

This project is licensed under the MIT License - see the LICENSE file for details.