Deltat = -- YOSHIDA-4
PE = --
H = --
Overview v2.0
GRAVITAS is a fully client-side, real-time N-body gravitational simulator embedded as a single self-contained HTML file. It runs entirely in the browser with no server, no dependencies, and no network requests after initial font loading.
Unlike simplified orbit visualisers that solve only two-body Kepler equations or use patched-conic approximations, GRAVITAS performs genuine numerical integration of the full N-body gravitational equations of motion, computing all N(N-1)/2 pairwise interactions at every integration stage.
Key Capabilities
- Four symplectic integrators -- Yoshida-4, Stormer-Verlet leapfrog, RKN-4, and Forest-Ruth-4, all switchable at runtime
- Adaptive timestepping -- Richardson extrapolation dynamically scales Deltat to maintain a user-defined tolerance
- Barycentric frame -- centre-of-mass position and velocity are zeroed at each preset load, preventing spurious drift
- Live conservation monitoring -- relative energy error DeltaE/E0 and angular momentum error DeltaL/L0 computed every frame
- Softened gravity -- regularised force prevents numerical blowup at close approaches
- 8 built-in presets -- from IAU-element Solar System to the Chenciner-Montgomery figure-8 choreography
- Interactive injection -- click the viewport to add a test particle with approximately circular orbital velocity
Units and Conventions
These Gaussian gravitational units are standard in solar system dynamics. They cause G.MSun = 4pi^2, so a body at 1 AU orbiting one solar mass has exactly a 1-year period -- no unit conversions needed in the code.
Physics Engine
Equations of Motion
Each body i of mass mi at position ri experiences the gravitational acceleration:
This is Newton's inverse-square law with a Plummer softening parameter eps that regularises the force at short range, preventing the divergence that would otherwise cause numerical blowup during close encounters.
Softening (Plummer Kernel)
The softened potential for a pair (i, j) is:
As eps->0 this recovers the exact Newtonian potential. The default eps = 0.01 AU is appropriate for planetary-scale simulations; for compact or stellar systems it should be reduced. The softening is adjustable via the slider.
Hamiltonian (Total Energy)
For a symplectic integrator on an autonomous Hamiltonian system, H should be nearly conserved -- it oscillates around a "shadow Hamiltonian" that is exponentially close to the true H for small enough Deltat. This near-conservation is what makes symplectic integrators superior for long-term orbital simulation.
Angular Momentum
In an isolated system with no external torques, L is exactly conserved. Deviations DeltaL/L0 measure rotational symmetry violation, which can arise from timestep truncation error or close encounters.
Barycentric Frame
After loading any preset, the code computes the centre-of-mass position and velocity:
and subtracts both from all bodies. This ensures no net linear drift of the system, matching the convention used in solar system ephemerides (ICRF barycentric frame).
Force Complexity
All N(N-1)/2 pairwise interactions are computed each force evaluation. This is O(N^2) per step -- exact but expensive. For N <= 15 on mobile hardware it runs at thousands of steps per second. For N > ~30, performance degrades; a Barnes-Hut tree (O(N log N)) would be needed but is outside this widget's scope.
Adaptive Timestep -- Richardson Extrapolation
Each candidate step of size h is validated by comparing:
- One step of size h (low-accuracy reference)
- Two steps of size h/2 (high-accuracy result, accepted)
The maximum relative position error across all bodies estimates the local truncation error. The next timestep is scaled:
where p is the integrator order (2 for leapfrog, 4 for Yoshida/RKN/Ruth). The safety factor 0.9 prevents over-ambitious steps. The result accepted is always the two-half-step solution, which is one order more accurate than the full-step solution (Richardson's principle).
Integrators
All four integrators are symplectic (or near-symplectic) -- they preserve the phase-space volume of the Hamiltonian flow, which is the mathematical foundation of long-term stability in orbital mechanics.
Why Symplectic Matters
Standard ODE solvers (RK4, Adams-Bashforth) are not symplectic. They introduce a secular energy drift -- a systematic one-way error growth -- that causes orbits to spiral in or out over long timescales. Symplectic integrators instead produce a bounded oscillation of the energy around a nearby "shadow Hamiltonian," keeping orbits stable for billions of timesteps.
1 . Stormer-Verlet Leapfrog (2nd order)
The simplest symplectic integrator. Uses a Kick-Drift-Kick (KDK) splitting:
Requires 1 force evaluation per step (the two half-kicks share one evaluation). Time-reversible and symplectic. Energy error oscillates bounded as O(h^2). Best for quick runs where accuracy is secondary.
2 . Yoshida 4th-Order Symplectic
Yoshida (1990) constructed higher-order symplectic integrators by composing leapfrog steps with specific coefficients. The 4th-order method uses 3 force evaluations with coefficients derived from the real root of 2w^3 - 1 = 0:
Energy error is O(h^4). For planetary orbits with Deltat ~ 0.005 yr this achieves DeltaE/E0 ~ 10^-^8 per step. The default integrator -- best balance of accuracy and speed.
3 . Runge-Kutta-Nystrom 4 (RKN-4)
RKN methods exploit the special structure of second-order systems x_ddot = f(x) -- they do not require storing half-step velocities. Three force evaluations:
4th-order accurate. Not strictly symplectic (it is an explicit Runge-Kutta method), but the near-symplectic structure of NKN methods gives it very good energy behaviour in practice. Excellent for high-accuracy short-run benchmarks.
4 . Forest-Ruth / Ruth-4 (4th order symplectic)
Independently derived by Ruth (1983) and Forest & Ruth (1990). Uses the same composition parameter theta = 1/(2-cbrt2) as Yoshida but in a 4-stage drift-kick-drift structure. Strictly symplectic, 4 force evaluations per step. Slightly more expensive than Yoshida-4 but with different error coefficients -- useful for cross-checking results.
Integrator Comparison
| Method | Order | Force evals | Symplectic | Best for |
|---|---|---|---|---|
| Leapfrog | 2nd | 1 | Yes | Fast runs, exploration |
| Yoshida-4 | 4th | 3 | Yes | General use (default) |
| RKN-4 | 4th | 3 | Near | Short accuracy checks |
| Ruth-4 | 4th | 4 | Yes | Cross-validation |
Controls Reference
Presets
| Preset | Description |
|---|---|
| Solar System (9) | Sun + 8 planets using approximate IAU J2000 orbital elements. Masses from DE430 ephemeris. Placed at perihelion with vis-viva velocity. |
| Figure-8 Choreo | The famous Chenciner & Montgomery (2000) choreographic solution: 3 equal masses chasing each other on a figure-8 curve. Exact published initial conditions. |
| Pythagorean 3-body | Three bodies with masses 3, 4, 5 at the vertices of a right triangle, released from rest. Classic chaotic problem -- highly sensitive to initial conditions. |
| Binary Star + Planets | Two unequal stars in mutual orbit with 4 test-mass planets in circumbinary orbits at various radii. |
| Lagrange Points | Sun-Jupiter system with test particles placed at L4 and L5 Trojan points, one slightly perturbed to illustrate long-term stability. |
| Chaotic 7-body | Seven bodies of varying mass arranged in a near-circular ring with tangential velocities -- rapidly becomes chaotic. |
| Mini Galaxy (12) | Massive central body + 11 lighter stars in a rough disk. Illustrates differential precession and orbital resonances. |
| Random (8) | Randomly generated 8-body system. Different every reset. |
Sliders
- Time scale -- multiplier on wall-clock time to simulation time. 1x = 1 yr/s real-time (approximately). Range 0.1x to 1000x.
- Softening eps -- Plummer softening in AU. Larger values smooth close encounters; smaller values are more physically realistic but risk numerical instability.
- Trail length -- number of position samples stored per body. 0 disables trails entirely.
- Zoom -- viewport scale multiplier. Also controllable via mouse wheel or pinch gesture.
- Tolerance 10^-n -- adaptive step error target. n=7 (10^-^7) is the default. Higher n = tighter tolerance = smaller steps = slower but more accurate.
- Base Deltat -- initial candidate timestep in years. The adaptive controller will modify this. Typical values 0.001-0.01 yr.
Buttons
- Pause / Play -- halts or resumes the integration loop. Canvas rendering continues while paused.
- Reset -- reloads the current preset from scratch. All trails, time, and step count are cleared.
- Trails -- toggles orbital trail rendering. Clearing trails frees memory.
- Vectors -- overlays velocity arrows on each body, scaled to the current zoom level.
- CoM -- draws a crosshair marker at the current centre-of-mass. Should remain near the viewport centre if the barycentric correction was applied correctly.
Canvas Interaction
- Click -- injects a low-mass test particle (m = 10^-^5 MSun) at that position with approximately circular orbital velocity around the total system mass.
- Mouse wheel -- zoom in/out.
- Pinch (touch) -- two-finger pinch to zoom on mobile.
Telemetry Colours
| Colour | Meaning for DeltaE/E0 |
|---|---|
| Green | Excellent -- < 10^-^6 (symplectic regime) |
| Blue | Good -- 10^-^6 to 10^-^4 |
| Amber | Caution -- 10^-^4 to 10^-^2 |
| Red | Poor -- > 10^-^2 (reduce Deltat or increase tolerance) |
Provenance & Governance
Authorship and Generation
This simulator was developed by Astrophyzix.
The physics algorithms implemented (Yoshida coefficients, RKN-4 formulation, Richardson extrapolation scheme, Plummer softening) are standard textbook methods whose mathematical content predates this tool by decades. The implementations are original code produced by Astrophyzix
Intended Use
GRAVITAS is intended for:
- Educational demonstration of N-body gravitational dynamics in browser environments
- Qualitative exploration of orbital mechanics, chaos, and integrator behaviour
It is not intended for:
- Production ephemeris computation or spacecraft trajectory planning
- Publication of quantitative results without independent verification
- Replacing validated astrodynamics software (REBOUND, MERCURY6, JPL Horizons)
Accuracy and Limitations
The Solar System preset uses approximate J2000 orbital elements for mass and semi-major axis but places each planet at perihelion with vis-viva velocity. Actual positions at any calendar date are not reproduced. Long-term (> 100 yr) integration of the solar system is chaotic (Laskar 1989) and any simulation diverges from reality on that timescale regardless of numerical method.
The Figure-8 choreography uses published exact initial conditions (Chenciner & Montgomery 2000) in normalised units. With Yoshida-4 at tolerance 10^-^8, the solution remains recognisable for ~100 orbital periods before chaotic divergence.
Numerical Validation
The following self-consistency checks are performed at runtime:
- Energy conservation: DeltaE/E0 displayed live. Yoshida-4 at default settings typically achieves 10^-^8 to 10^-^6 for smooth orbits.
- Angular momentum conservation: DeltaL/L0 displayed live. Should track energy error in magnitude.
- Barycentric check: CoM overlay should remain stationary near canvas centre after preset load.
External validation was not performed against a reference integrator (e.g., REBOUND's IAS15) as part of this build. Users who require validated results should cross-check against established sources.
Licence and Redistribution
This file is provided as-is for educational and personal use by Astrophyzix Digital Observatory . No warranty is given regarding numerical accuracy or fitness for any purpose. The module code and documentation must not be copied or redistributed without written permission. For non-commercial educational use on the Astrophyzix platform. Copyright Registered 2026
Change Log
References
The following publications are the primary sources for the algorithms implemented in GRAVITAS.