All three integrators solve the same 2-body gravity ODE, starting from the identical initial state (r₀, v₀) on a Kepler ellipse with semi-major axis a and eccentricity e:
d²r/dt² = -GM · r / |r|³
Energy: E = ½|v|² - GM/|r| (conserved exactly by the true orbit)
Each method advances the state by one step h differently:
Euler: r,v += h·f(r,v) → error O(h) per step
Midpoint: k = f(r + h/2·v, ...) at half-step → error O(h²) per step
RK4: weighted average of 4 slope evaluations → error O(h⁴) per step
The convergence test re-runs each integrator from scratch at eight step sizes between 0.002 and 0.05, every run stepping the identical Kepler orbit for a quarter period — short enough that every method's error stays small relative to the orbit, which is exactly the small-h regime the O(h^p) theory describes. Because this is a 2-body problem, the true position at the exact time each run stops is not approximated — it is computed by solving Kepler's equation M = E − e·sin(E) for the eccentric anomaly with Newton's method, then evaluating the exact ellipse position. The Euclidean distance between each run's final numerical position and that exact analytic position is its global error. Plotting log(error) against log(h) and fitting a least-squares line to the eight points gives the measured convergence order directly as the line's slope — theory says it should land near 1 for Euler, 2 for RK2 and 4 for RK4, and the chart lets you check that a hundredfold change in accuracy really does come from roughly the step size squared (RK2) or to the fourth power (RK4), not just "RK4 looks smoother."
- Step size h — larger steps make truncation error explode fastest for Euler, slower for RK2, and barely visible for RK4 until h gets large.
- Eccentricity e — a more elongated ellipse swings through perihelion faster, stressing the fixed step size exactly where curvature is highest; re-run the convergence test after changing it.
- Energy drift — since the true orbit conserves energy exactly, (E(t) − E₀)/|E₀| is a direct, visible measure of each method's accumulated numerical error during the live race.
This is the same trade-off every real ODE solver (orbital mechanics, circuit simulation, molecular dynamics) has to make between step size, computational cost and long-term accuracy.