Course contents · Foundations · Selection guide
These explanations describe the local implementations. Read Euler, RK4, and Leapfrog first; return to the other methods as their experiments introduce them. IAS15 and external BS are separately called reference implementations.
Contents: Overview · Euler · RK4 · Leapfrog · RK45 · Hermite4 · BS · GL4 · WH2 · WH4 · External references
| Method | Main idea | Time control here | Key implementation boundary |
|---|---|---|---|
| Euler | One forward derivative estimate | Fixed shared step | First-order teaching baseline |
| RK4 | Four weighted derivative stages | Fixed shared step | Explicit, fourth-order; not symplectic |
| Leapfrog | Kick, drift, kick | Fixed shared step | Second-order; position-dependent forces; no endpoint-force cache |
| RK45 | Embedded Dormand–Prince 5(4) pair | Adaptive shared step | Local error controller, not a global-error guarantee |
| Hermite4 | Acceleration-and-jerk predictor/corrector | Fixed shared step | One PEC correction; no block or individual steps |
| BS | Modified midpoint plus extrapolation | Adaptive shared step | Variable extrapolation work; local implementation distinct from external BS |
| GL4 | Two implicit Gauss collocation stages | Fixed shared step | Dense finite-difference Newton; small educational systems |
| WH2 | Jacobi Kepler/interaction splitting | Fixed shared step | Positive-mass unsoftened planetary path, with a distinct restricted path |
| WH4 | Symmetric composition of WH2 | Fixed shared step | Includes a negative substep; not WHFast |
Let y collect positions and velocities and let
Read: implementation · 1A
For gravity, update position with the old velocity and velocity with the old acceleration. This is the tangent-line approximation at the current state. With smooth solutions over a fixed interval, its global error is first order. It is useful because each operation is visible, not because it is a competitive precision orbit solver.
The local implementation uses one force-field evaluation per step. A short Euler run illustrates accumulated error; do not compare its cost with a longer task performed by another method.
Read: implementation · 1B
Each stage estimates a derivative at an intermediate state. Position and velocity must advance together in those stage states. The weights resemble a quadrature pattern, but RK4 is not obtained by applying Simpson's rule to known acceleration samples: the intermediate states themselves are estimates coupled to the differential equation.
The local implementation evaluates four force fields per step. Fourth-order convergence is a resolved smooth-problem property, not a promise of exact energy or long-term phase.
Read: implementation · 1D
For position-dependent acceleration, the kick–drift–kick form is
Each kick changes velocity while holding position fixed; the drift changes position with the intermediate velocity. For the separable autonomous Hamiltonian setting and fixed h, this is a symmetric symplectic second-order method. Symplectic means preserving the canonical phase-space form; it does not mean preserving the original energy exactly. See Hairer's geometric-integration notes for the mathematical framework.
This teaching implementation makes two force calls every step and does not cache the final acceleration for the next step. The stored counts describe that implementation. An optimized cached implementation can change cost without changing the nominal order.
Read: implementation · driver · 1C
The Dormand–Prince pair forms fifth- and fourth-order approximations from related stages. Their difference estimates local error. The driver scales each component using absolute and relative tolerances, accepts or rejects the trial, and adjusts the next step. The published tableau is attributed in Sources.
The important distinction is between trial and accepted state. Calling an adaptive trial function repeatedly and accepting every return would bypass the error controller. The experiments use the integration driver for adaptive methods.
One shared timestep advances all particles. Componentwise error scaling can change under rotations or added zero coordinates, which is explored in 1F. Requested tolerance, accepted dt, and achieved global error are three different measurements. Rejected work still counts.
Read: implementation · 2E
Jerk is the derivative of acceleration,
Evaluate predicted acceleration
This local fourth-order PEC implementation makes two acceleration-plus-jerk evaluations, with one correction. It has no derivative cache, iteration to time symmetry, Ahmad–Cohen neighbor scheme, or individual/block timestep scheme. Do not transfer those production-Hermite advantages to this fixed-shared-step experiment.
Read: implementation · 1C
The local Bulirsch–Stoer-family method divides a proposed macro-step into modified-midpoint substeps. It repeats with increasing substep counts, then extrapolates the approximations toward zero substep size. For a smooth problem the symmetric midpoint construction has a useful even-power error expansion, allowing extrapolation to cancel successive terms.
The solver compares successive extrapolation estimates and adjusts the macro-step through the driver. One accepted step can therefore require many force evaluations and extrapolation columns. Compare those costs and rejected trials, not merely the number of accepted output points.
This is an independently written educational polynomial-extrapolation implementation. The external REBOUND BS reference is a different implementation; the name alone does not make their costs or behavior identical.
Read: implementation · 3B
Two-stage Gauss–Legendre collocation solves coupled stage states:
with
Unlike explicit RK4, each stage depends on the other unknown stage. The code uses Newton iteration, a finite-difference Jacobian, a dense linear solve, and damping. The exact collocation method is fourth-order, symmetric, and symplectic; finite nonlinear residuals and roundoff perturb those properties. The standard theory is discussed in Hairer's notes.
In 3D there are 12N stage unknowns. Thus two nominal stages do not mean two cheap force evaluations. Track Jacobians, Newton iterations, and linear solves. GL4's large-N budget exclusions in 3B are planned omissions, not proof of mathematical inapplicability.
Read: implementation · Jacobi transforms · Kepler solver · 2A
A planetary system can be decomposed into dominant Kepler motion A plus interactions B. WH2 composes a half interaction kick, a full Kepler drift, and a half kick. The local positive-mass path uses Jacobi coordinates: each body's coordinate is relative to the center of mass of all preceding bodies. Ordering the hierarchy therefore matters.
The published Wisdom–Holman idea is Hamiltonian splitting; the local code is a teaching implementation, not a rebranding of WHFast. In the two-body limit B vanishes, so the Kepler solve dominates the remaining error. With nonzero perturbations, 2B measures splitting convergence.
A close encounter can make the interaction part large. The current unsoftened positive-mass planetary path does not apply to the softened equal-mass cluster in 3D/3E. A separate restricted implementation handles prescribed moving primaries and time-dependent kicks in 2R.
Read: implementation · 2B
Compose the symmetric second-order map
The coefficients satisfy
WH4 performs more subflows and Kepler work than WH2. Higher order can improve accuracy when the split is appropriate and sufficiently resolved, but does not make every encounter safe or every long integration cheapest. The coefficient source is Yoshida (1990), linked in Sources.
The course calls REBOUND IAS15 and external BS through references.py. IAS15 is a high-order adaptive production method described in the original paper; its presence here supplies an external comparison, not automatic ground truth.
Part 1 can test numerical references against an analytic solution. Later lessons compare reference settings and an independent implementation, then limit pointwise interpretation to their checked agreement window. The adapter also records the force-count convention and unavailable counters. Production references are labeled separately from the nine course methods; implementation language is not used to create an algorithmic wall-time ranking.
Order describes how error changes under refinement. It does not count the work in a step. RK4 evaluates explicit stages, GL4 solves a coupled nonlinear system, Hermite4 evaluates jerk, and WH4 composes Kepler and interaction flows. Their assumptions and error constants differ too. The useful comparison is the measured work needed to reach the same observable target.
Apply these ideas in the selection guide
Lesson 4A adds a third-order Wisdom corrector to the existing Jacobi WH2 kernel through an explicit BAB/ABA conjugacy. The dedicated processed.py loop transforms the initial state once and inverse-transforms output copies. The original nine-method registry and step functions are unchanged.
The kernel remains second order. A finite processor can suppress a leading perturbation-linear term while leaving perturbation-quadratic and higher-step terms. Initial and output work is counted separately. This implementation supports fixed steps, positive ordered masses, point gravity, and the existing Jacobi split. Step changes and close-encounter validity are not established by the pilot.
REBOUND WHFast-C3, MERCURIUS and TRACE are production comparison methods, not additions to the nine local registry entries. WHFast uses the declared Jacobi coordinates; hybrid adapters retain their declared democratic-heliocentric architecture and switching settings. IAS15 and external BS construct checked numerical references. Their role is distinct from a candidate even when the same library provides both. Unavailable external work counters remain unavailable. See 4C, 4E and sources.