Planetary motion: comparing numerical integrators
How much does an integration method change the orbit we predict? Starting from Newtonian gravity, we implement six methods, compare energy and trajectory error, and animate a three-body system. The examples use dimensionless units with $G=1$. They are idealised dynamics, not Solar System ephemerides.
Prerequisites: vectors, Newton's laws, basic differential equations and NumPy.
Install projects/requirements.txt, then run all cells in order.
The original filename is retained so existing links continue to work.
Equations and conserved quantities
For positions $r_i$, momenta $p_i=m_i v_i$ and masses $m_i$, $$\dot r_i=p_i/m_i,\qquad \dot p_i=Gm_i\sum_{j\ne i}m_j\frac{r_j-r_i}{(|r_j-r_i|^2+\epsilon^2)^{3/2}}.$$ We use a fixed Plummer softening length $\epsilon=0.01$ to avoid a singular force at close approach. The consistent potential is $$V=-\sum_{i<j}\frac{Gm_i m_j}{\sqrt{|r_i-r_j|^2+\epsilon^2}}.$$ This changes the physical model at small separation; it is not a collision model. Kinetic energy is $T=\sum_i |p_i|^2/(2m_i)$, and the Hamiltonian is $H=T+V$. The Lagrangian $L=T-V$ gives the same equations through Euler–Lagrange equations. Total linear momentum $\sum p_i$ and angular momentum $\sum r_i\times p_i$ are conserved in this isolated model. We compare numerical errors against these laws.
From Newton's force law to the equations of motion
For two point masses separated by a distance $r$, the magnitude of the Newtonian force is $F=Gm_1m_2/r^2$. The force is attractive, so its vector form points from each body towards the other. Adding the contributions from all other bodies gives $$m_i\ddot{\mathbf r}_i= -Gm_i\sum_{j\ne i}m_j\frac{\mathbf r_i-\mathbf r_j} {|\mathbf r_i-\mathbf r_j|^3}.$$ With $\mathbf p_i=m_i\dot{\mathbf r}_i$, each second-order position equation becomes two first-order equations: one for position and one for momentum. Those are the quantities advanced by the numerical methods below.
The implementation replaces the singular denominator with $(|\mathbf r_i-\mathbf r_j|^2+\epsilon^2)^{3/2}$. The following derivations use that same softened potential, so the force calculation and energy diagnostic refer to the same physical model. Setting $\epsilon=0$ recovers the point-mass expressions away from a collision.
Lagrangian formulation
The kinetic and potential energies are $$T=\sum_i\frac12m_i|\dot{\mathbf r}_i|^2,\qquad V=-\sum_{i<j}\frac{Gm_im_j} {\sqrt{|\mathbf r_i-\mathbf r_j|^2+\epsilon^2}}.$$ Each pair occurs once in $V$. The Lagrangian is $L=T-V$. For a Cartesian component $q_i$ of body $i$, the Euler–Lagrange equation is $$\frac{d}{dt}\left(\frac{\partial L}{\partial\dot q_i}\right) -\frac{\partial L}{\partial q_i}=0.$$ Because the kinetic energy is quadratic in velocity, $$\frac{\partial L}{\partial\dot q_i}=m_i\dot q_i,\qquad \frac{d}{dt}\left(\frac{\partial L}{\partial\dot q_i}\right)=m_i\ddot q_i.$$ Differentiating the potential gives $$\frac{\partial L}{\partial q_i}=-\frac{\partial V}{\partial q_i} =-\sum_{j\ne i}\frac{Gm_im_j(q_i-q_j)} {(|\mathbf r_i-\mathbf r_j|^2+\epsilon^2)^{3/2}}.$$ Substituting these derivatives reproduces Newton's equation for that component. Repeating for $x$, $y$ and $z$ supplies the full equations for every body.
Hamiltonian formulation and the numerical state
The canonical momentum is $\mathbf p_i=\partial L/\partial\dot{\mathbf r}_i =m_i\dot{\mathbf r}_i$. The Legendre transform gives $$H=\sum_i\mathbf p_i\cdot\dot{\mathbf r}_i-L=T+V =\sum_i\frac{|\mathbf p_i|^2}{2m_i} -\sum_{i<j}\frac{Gm_im_j}{\sqrt{|\mathbf r_i-\mathbf r_j|^2+\epsilon^2}}.$$ Hamilton's equations then separate the two derivatives needed by the code: $$\dot{\mathbf r}_i=\frac{\partial H}{\partial\mathbf p_i} =\frac{\mathbf p_i}{m_i},$$ $$\dot{\mathbf p}_i=-\frac{\partial H}{\partial\mathbf r_i} =-\sum_{j\ne i}\frac{Gm_im_j(\mathbf r_i-\mathbf r_j)} {(|\mathbf r_i-\mathbf r_j|^2+\epsilon^2)^{3/2}}.$$ The kinetic term moves positions while holding momenta fixed; the potential term changes momenta while holding positions fixed. This split motivates the drift and kick steps in symplectic Euler and Verlet. Symplectic structure does not imply exact conservation of the original Hamiltonian at a finite step size.
The isolated continuous model conserves energy because $H$ has no explicit time dependence. Its pairwise central forces also conserve total linear and angular momentum. Tracking these quantities gives checks on the numerical solution, while a trajectory comparison tests errors that conservation alone can miss.
Units and assumptions
A length scale $L$, mass scale $M$ and time scale $\sqrt{L^3/(GM)}$ set $G=1$. All plot axes below use these dimensionless units. We omit relativity, radiation, finite-body collisions and external forces. Pairwise force evaluation costs $O(n^2)$, which is adequate for the two- and three-body examples here.
Set up the initial conditions
Show or hide code
import numpy as np
import numpy.typing as npt
import matplotlib.pyplot as plt
import scipy.integrate as integrate
from enum import IntEnum
class Integrator(IntEnum):
EULER=0
SYMPLECTIC_EULER=1
VERLET=2
RK2=3
RK4=4
IVP=5
class SolarSystem:
"""
Class holding the solar system state
:param name: Name of the solar system
:param n_bodies: Number of bodies in the solar system
:param G: Gravitational constant
:param M: Mass of the bodies
:param r: Position of the bodies
:param p: Momentum of the bodies
:param R: Radius of the bodies
:param rho: Density of the bodies -- calculated based on other properties assumed constant
"""
def __init__(self, name, n_bodies, G = 1., M = None, r = None, p = None, R = None, rho = None):
self.name = name
self.n_bodies = n_bodies
self.G = G
self.softening = 0.01
self.M = M if M is not None else np.ones( shape=n_bodies )
self.r = r if r is not None else np.zeros( shape=(n_bodies, 3))
self.p = p if p is not None else np.zeros( shape=(n_bodies, 3))
self.R = R if R is not None else np.ones( shape=n_bodies) * 1e-3
self.rho = 3*self.M/(4*np.pi*self.R**3)
# ───────── derived quantities ─────────
def energy(self): return calc_energy(self)
def kinetic(self ): return calc_kinetic_energy(self)
def potential(self): return calc_potential_energy(self)
def centre_of_mass(self): return calc_centre_of_mass(self)
def angular_momentum(self): return calc_angular_momentum(self)
def step(self, dt: float, scheme: Integrator = Integrator.VERLET) -> None:
"""
Advance the system by one time-step *dt* using the chosen Integrator.
All methods use the same softened force law.
"""
# Fixed-step branches
if scheme == Integrator.EULER:
self.r, self.p = calc_step_euler(self, dt)
elif scheme == Integrator.SYMPLECTIC_EULER:
self.r, self.p = calc_step_euler_symplectic(self, dt)
elif scheme == Integrator.VERLET:
self.r, self.p = calc_step_Stormer_Verlet(self, dt)
elif scheme == Integrator.RK2:
self.r, self.p = calc_step_rk2(self, dt)
elif scheme == Integrator.RK4:
self.r, self.p = calc_step_rk4(self, dt)
# Adaptive solver
elif scheme == Integrator.IVP:
#
# Everything inside the `with objmode()` block runs in normal
# CPython space, so SciPy is allowed.
#
self.r, self.p = calc_step_ivp(self, dt)
else:
raise ValueError(f"Unknown integrator id {scheme}")
Conserved quantities
Show or hide code
def calc_energy(ss: SolarSystem) -> float:
"""
Calculates the energy of the system
:param ss: SolarSystem object
:return: energy of the system
"""
return calc_kinetic_energy(ss) + calc_potential_energy(ss)
def calc_potential_energy(ss: SolarSystem) -> float:
"""
Calculates the potential energy of the system
:param ss: SolarSystem object
:return: potential energy of the system
"""
U = 0.
for i in range(ss.n_bodies):
for j in range(i+1, ss.n_bodies):
distance = np.sqrt(np.sum((ss.r[j]-ss.r[i])**2) + ss.softening**2)
U -= ss.G*ss.M[i]*ss.M[j]/distance
return U
def calc_grav_potential(ss: SolarSystem, pos: np.array) -> float:
"""
Calculates the potential energy of the system
:param ss: SolarSystem object
:param pos: position of the body to calculate potential energy for
:return: potential energy of the body at the given position
"""
distance = np.sqrt(np.sum((ss.r-pos)**2, axis=1) + ss.softening**2)
return -ss.G*np.sum(ss.M/distance)
def calc_kinetic_energy(ss: SolarSystem) -> float:
"""
Calculates the kinetic energy of the system
:param ss: SolarSystem object
:return: kinetic energy of the system
"""
return 0.5 * np.sum( np.sum(ss.p**2, axis=1)/ss.M )
def calc_centre_of_mass(ss: SolarSystem) -> np.array:
"""
Calculates the centre of mass of the system
:param ss: SolarSystem object
"""
return np.sum(ss.r*ss.M[:,None], axis=0)/np.sum(ss.M)
def calc_angular_momentum(ss: SolarSystem) -> np.array:
"""
Calculates the angular momentum of the system
:param ss: SolarSystem object
"""
return np.sum(np.cross(ss.r, ss.p), axis=0)
Show or hide code
# Simple test system
two_body = SolarSystem("Two body", 2)
two_body.r[0] = np.array( [-0.5, 0, 0] )
two_body.r[1] = np.array( [ 0.5, 0, 0] )
two_body.p[0] = np.array( [ 0,-0.5, 0] )
two_body.p[1] = np.array( [ 0,+0.5, 0] )
two_body.M[0] = 1
two_body.M[1] = 1
print( 'Energy ' , two_body.energy() )
print( 'T ' , two_body.kinetic() )
print( 'V ' , two_body.potential() )
print( 'COM ' , two_body.centre_of_mass() )
print( 'PA ' , two_body.angular_momentum() )
Show or hide code
# The system we are going to look at
three_body = SolarSystem("Three body", 3)
three_body.r[0] = np.array( [ -0.97000436, 0.24308753, 0])
three_body.r[1] = np.array( [ +0.97000436, -0.24308753, 0])
three_body.r[2] = np.array( [ 0, 0, 0])
three_body.p[0] = np.array( [0.4662036850, 0.4323657300, 0 ] )
three_body.p[1] = np.array( [0.4662036850, 0.4323657300, 0 ] )
three_body.p[2] = np.array( [ -0.93240737, -0.86473146, 0 ] )
Update momentum and position
Show or hide code
# ───────── packing helpers (vectorised) ─────────
def pack_state(r: np.ndarray, p: np.ndarray) -> np.ndarray:
"""Flatten and concatenate positions + momenta into one 1-D array."""
return np.hstack((r.ravel(), p.ravel()))
def reshape_state(arr: np.ndarray) -> np.ndarray:
"""Reshape a 2-D array into a 1-D array with positions and momenta."""
return arr.ravel()
def unpack_state(v: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
"""Inverse of pack_state → (r, p)."""
n = v.size // 6
return v[:3*n].reshape((n, 3)), v[3*n:6*n].reshape((n, 3))
Show or hide code
def calc_dpdt ( ss: SolarSystem ) -> np.array:
"""
Calculates the derivative of the position vector with respect to time
:param ss: Solar system
:return: derivative of position vector with respect to time
"""
state = pack_state(ss.r, ss.p)
return calc_dpdt_v( state, ss )
def calc_dpdt_v(state: np.ndarray, ss: SolarSystem) -> np.ndarray:
"""
Calculates force (the momentum derivative)
:param state: state vector
:param ss: Solar system
:return: derivative of position vector with respect to time
"""
r, p = unpack_state(state)
n = ss.n_bodies
dpdt = np.zeros_like(r)
# broadcast pair-wise differences
for i in range(n):
dr = r - r[i] # shape (n, 3)
d2 = np.sum(dr * dr, axis=1) # squared distance
mask = np.arange(n) != i # skip self
inv_r3 = ss.G * ss.M[mask] / np.power(d2[mask] + ss.softening**2, 1.5)
dpdt[i] = ss.M[i] * (dr[mask].T * inv_r3).T.sum(axis=0)
return dpdt
def calc_drdt ( ss: SolarSystem ) -> np.array:
"""
Calculates the derivative of the position vector with respect to time
:param ss: Solar system
:return: derivative of position vector with respect to time
"""
state = pack_state(ss.r, ss.p)
return calc_drdt_v( state, ss )
def calc_drdt_v ( state: np.array, ss:SolarSystem ) -> np.array:
r, p = unpack_state(state)
return p / ss.M[:,None]
def calc_derivatives( t, y0 : np.array, ss : SolarSystem ) -> np.array :
"""
Calculate the derivatives of the state vector.
""" # calculate the derivatives of the state vector
dpdt = calc_dpdt_v( y0, ss )
drdt = calc_drdt_v( y0, ss )
return pack_state(drdt, dpdt)
Show or hide code
drdt = calc_drdt(two_body)
dpdt = calc_dpdt(two_body)
print(drdt)
print(dpdt)
# Unequal masses expose accidental confusion between force and acceleration.
unequal = SolarSystem('Unequal masses', 2)
unequal.M[:] = [2., 3.]
unequal.r[:] = [[-0.5, 0., 0.], [0.5, 0., 0.]]
forces = calc_dpdt(unequal)
np.testing.assert_allclose(forces.sum(axis=0), np.zeros(3), atol=1e-12)
expected = 6./(1.+unequal.softening**2)**1.5
np.testing.assert_allclose(forces[0], [expected, 0., 0.])
unequal.p[:] = [[0., -2., 0.], [0., 2., 0.]]
np.testing.assert_allclose(unequal.angular_momentum(), [0., 0., 2.])
Explicit Euler
$$r_{n+1}=r_n+h\,p_n/m,\qquad p_{n+1}=p_n+hF(r_n).$$ Both derivatives are evaluated at the old state. The method is first order and usually drifts in energy on orbital problems.
Show or hide code
def calc_step_euler( ss: SolarSystem, step_size: float) -> tuple[npt.NDArray,npt.NDArray]:
"""
Calculates a step in time using the Euler method
:param ss: Solar system
:param step_size: time step size
"""
y0=pack_state(ss.r, ss.p)
dydt = calc_derivatives( 0, y0 , ss )
return unpack_state( y0 + (dydt * step_size) )
Symplectic Euler (drift then kick)
$$r_{n+1}=r_n+h\,p_n/m,\qquad p_{n+1}=p_n+hF(r_{n+1}).$$ Updating momentum with the new position gives a first-order symplectic method. Symplectic does not mean exact energy conservation.
Show or hide code
def calc_step_euler_symplectic( ss: SolarSystem, step_size: float) -> tuple[npt.NDArray,npt.NDArray]:
"""
Calculates a step in time using the first-order symplectic methods Euler-A
:param ss: Solar system
:param step_size: time step size
"""
y0=pack_state(ss.r, ss.p)
drdt = calc_drdt_v( y0, ss )
size = drdt.size
dydt = np.zeros( shape = size * 2)
dydt[0:size] = reshape_state( drdt )
y05 = y0 + dydt * step_size
dpdt = calc_dpdt_v( y05, ss )
dydt[0:size] = 0
dydt[size:size*2] = reshape_state( dpdt )
return unpack_state( y05 + (dydt * step_size) )
Störmer–Verlet (kick, drift, kick)
$$p_{n+1/2}=p_n+\tfrac h2F(r_n),$$ $$r_{n+1}=r_n+h\,p_{n+1/2}/m,$$ $$p_{n+1}=p_{n+1/2}+\tfrac h2F(r_{n+1}).$$ This is second order and time reversible for a fixed step. Its energy error usually oscillates rather than drifting steadily on a well-resolved bound orbit.
Show or hide code
def calc_step_Stormer_Verlet( ss: SolarSystem, step_size: float) -> tuple[npt.NDArray,npt.NDArray] :
"""
A second-order accurate symplectic method is the Störmer-Verlet method.
Algorithm:
1. p_{1/2} = p_0 + 0.5 * dt * F(r_0)
2. r_1 = r_0 + dt * p_{1/2} / m
3. p_1 = p_{1/2} + 0.5 * dt * F(r_1)
"""
# Step 1: Half momentum update
force_0 = calc_dpdt_v(pack_state(ss.r, ss.p), ss)
p_half = ss.p + 0.5 * step_size * force_0
# Step 2: Full position update
r_new = ss.r + step_size * p_half / ss.M[:, None]
# Step 3: Second half momentum update
force_1 = calc_dpdt_v(pack_state(r_new, p_half), ss)
p_new = p_half + 0.5 * step_size * force_1
return r_new, p_new
Runge–Kutta and adaptive integration
Heun's RK2 averages the initial derivative with an Euler trial derivative. Classical RK4 combines four stage evaluations. These methods can give small short-term trajectory errors but do not preserve symplectic structure. The SciPy comparison uses adaptive RK45 with explicit tight tolerances.
Show or hide code
def calc_step_rk2( ss: SolarSystem, step_size: float) -> tuple[npt.NDArray,npt.NDArray]:
"""
Calculates a step in time using RungeKutta 2
:param ss: Solar system
:param step_size: time step size
"""
y0=pack_state(ss.r, ss.p)
k1 = calc_derivatives( 0, y0 , ss ) * step_size
k2 = calc_derivatives( 0, y0+k1 , ss ) * step_size
return unpack_state( y0+(k1+k2)*0.5 )
def calc_step_rk4( ss: SolarSystem, step_size: float) -> tuple[npt.NDArray,npt.NDArray]:
"""
Calculates a step in time using RungeKutta 4
:param ss: Solar system
:param step_size: time step size
"""
y0=pack_state(ss.r, ss.p)
k1 = calc_derivatives( 0, y0 , ss )
k2 = calc_derivatives( 0, y0+k1*step_size*0.5 , ss )
k3 = calc_derivatives( 0, y0+k2*step_size*0.5 , ss )
k4 = calc_derivatives( 0, y0+k3*step_size , ss )
return unpack_state( y0 + (k1 + 2*k2 + 2*k3 + k4)*step_size/6 )
Show or hide code
def calc_step_ivp(ss, dt):
return unpack_state(calc_step_ivp_v(pack_state(ss.r, ss.p), dt, ss))
def calc_step_ivp_v(y0, dt, ss):
sol = integrate.solve_ivp(calc_derivatives, [0., dt], y0, args=(ss,),
method='RK45', rtol=1e-10, atol=1e-12)
if not sol.success:
raise RuntimeError(sol.message)
return sol.y[:, -1]
Show or hide code
# Compare one small step without changing the original system.
for method in Integrator:
# Each function returns new position and momentum arrays.
functions = [calc_step_euler, calc_step_euler_symplectic, calc_step_Stormer_Verlet,
calc_step_rk2, calc_step_rk4, calc_step_ivp]
r_next, p_next = functions[method.value](two_body, 0.001)
print(f'{method.name:18s} |delta r|={np.linalg.norm(r_next-two_body.r):.6e}; '
f'|delta p|={np.linalg.norm(p_next-two_body.p):.6e}')
Show or hide code
def run_iters_and_store(ss_0: SolarSystem,
method, # Can be IntegrationMethod enum or string
num_iterations: int,
step_size: float) -> dict:
"""
Integrate a *copy* of ss_0 with the chosen scheme.
Args:
method: Can be either IntegrationMethod enum or string
"""
ss = SolarSystem(ss_0.name, ss_0.n_bodies, ss_0.G,
ss_0.M.copy(), ss_0.r.copy(), ss_0.p.copy(),
ss_0.R.copy(), ss_0.rho.copy())
results_r = np.zeros((num_iterations, ss.n_bodies, 3))
results_p = np.zeros_like(results_r)
results_E = np.zeros((num_iterations, 1))
results_T = np.zeros_like(results_E)
results_V = np.zeros_like(results_E)
results_COM = np.zeros((num_iterations, 3))
results_AngMom = np.zeros_like(results_COM)
# initial state
results_r[0] = ss.r
results_p[0] = ss.p
results_E[0] = ss.energy()
results_T[0] = ss.kinetic()
results_V[0] = ss.potential()
results_COM[0] = ss.centre_of_mass()
results_AngMom[0] = ss.angular_momentum()
if num_iterations < 2 or step_size <= 0:
raise ValueError('Use at least two samples and a positive time step.')
grid = np.arange(num_iterations)*step_size
if method == Integrator.IVP:
sol = integrate.solve_ivp(calc_derivatives, [0., grid[-1]], pack_state(ss.r, ss.p),
t_eval=grid, args=(ss,), rtol=1e-10, atol=1e-12)
if not sol.success:
raise RuntimeError(sol.message)
for i in range(1, num_iterations):
if method == Integrator.IVP:
ss.r, ss.p = unpack_state(sol.y[:, i])
else:
ss.step(step_size, method)
results_r[i] = ss.r
results_p[i] = ss.p
results_E[i] = ss.energy()
results_T[i] = ss.kinetic()
results_V[i] = ss.potential()
results_COM[i] = ss.centre_of_mass()
results_AngMom[i] = ss.angular_momentum()
return dict(time=grid, r=results_r, p=results_p, E=results_E, T=results_T,
V=results_V, COM=results_COM, AngMom=results_AngMom)
The fixed-step methods below all start from independent copies of the same system. The adaptive reference integrates the full interval once and samples the same time grid.
Show or hide code
x_vals = np.linspace(-1.5, 1.5, 400)
plt.figure(figsize=(7, 3))
plt.plot(x_vals, [calc_grav_potential(two_body, np.array([xx, 0., 0.])) for xx in x_vals])
plt.xlabel('x [dimensionless]'); plt.ylabel('Potential per unit mass')
plt.title('Softened two-body potential along the x axis'); plt.tight_layout(); plt.show()
Compare the methods on a bound two-body orbit
The initial orbit is eccentric. We use the same time grid for every method and record energy, centre of mass and angular momentum. No method is selected by runtime alone.
Show or hide code
import time
num_iterations, step_size = 2001, 0.01
runs = {}
for method in Integrator:
start = time.perf_counter()
runs[method] = run_iters_and_store(two_body, method, num_iterations, step_size)
drift = np.max(np.abs(runs[method]['E']-runs[method]['E'][0]))
print(f'{method.name:18s} max |delta E|={drift:.3e}')
results = run_iters_and_store(three_body, Integrator.VERLET, 2001, 0.005)
Show or hide code
def plot_results(result, num_iterations=None):
sl = slice(None, num_iterations)
t = result['time'][sl]
fig, ax = plt.subplots(2, 3, figsize=(12, 6))
r, p = result['r'][sl], result['p'][sl]
ax[0, 0].plot(r[:, 0, 0], r[:, 0, 1]); ax[0, 0].plot(r[:, 1, 0], r[:, 1, 1])
ax[0, 0].set(xlabel='x', ylabel='y', title='Orbital paths', aspect='equal')
ax[0, 1].plot(t, p[:, 0, 1]); ax[0, 1].set(title='Body 0 momentum', ylabel='p_y')
ax[0, 2].plot(t, result['E'][sl]); ax[0, 2].set(title='Total energy')
ax[1, 0].plot(t, (result['E'][sl]-result['E'][0])); ax[1, 0].set(title='Energy error')
ax[1, 1].plot(t, np.linalg.norm(result['COM'][sl]-result['COM'][0], axis=1)); ax[1, 1].set(title='Centre-of-mass displacement')
ax[1, 2].plot(t, np.linalg.norm(result['AngMom'][sl]-result['AngMom'][0], axis=1)); ax[1, 2].set(title='Angular-momentum error')
for axis in ax.flat:
if axis is not ax[0, 0]: axis.set_xlabel('Time')
fig.tight_layout(); plt.show()
Show or hide code
fig, axes = plt.subplots(1, 2, figsize=(11, 4))
reference = runs[Integrator.IVP]
for method, result in runs.items():
if method == Integrator.IVP:
continue
axes[0].plot(result['time'], (result['E']-result['E'][0]).ravel(), label=method.name)
error = np.linalg.norm(result['r']-reference['r'], axis=(1, 2))
axes[1].semilogy(result['time'][1:], np.maximum(error[1:], 1e-16), label=method.name)
axes[0].set(xlabel='Time', ylabel='Energy error', title='Energy conservation')
axes[1].set(xlabel='Time', ylabel='Position error vs RK45', title='Trajectory agreement')
for ax in axes: ax.legend(fontsize=8)
plt.tight_layout(); plt.show()
# Check RK4 convergence over a shorter fixed interval; this measures trajectory error.
coarse = run_iters_and_store(two_body, Integrator.RK4, 101, 0.02)
fine = run_iters_and_store(two_body, Integrator.RK4, 201, 0.01)
ref = run_iters_and_store(two_body, Integrator.IVP, 201, 0.01)
coarse_error = np.linalg.norm(coarse['r'][-1]-ref['r'][-1])
fine_error = np.linalg.norm(fine['r'][-1]-ref['r'][-1])
assert fine_error < coarse_error/5
print('RK4 final-position error, h / h/2:', coarse_error, fine_error)
Show or hide code
from IPython.display import HTML, display
import matplotlib.animation as animation
def plot_xy(result, n_to_plot=1000):
fig, ax = plt.subplots(figsize=(5, 5))
for i in range(result['r'].shape[1]):
ax.plot(result['r'][:n_to_plot, i, 0], result['r'][:n_to_plot, i, 1], label=f'Body {i}')
ax.set(xlabel='x', ylabel='y', aspect='equal', title='Three-body paths')
ax.legend(); plt.tight_layout()
return fig, ax
plot_xy(results, len(results['time'])); plt.show()
fig, ax = plt.subplots(figsize=(5, 5))
positions = results['r']
limit = np.max(np.abs(positions[:, :, :2]))*1.1
ax.set(xlim=(-limit, limit), ylim=(-limit, limit), xlabel='x', ylabel='y', aspect='equal')
lines = [ax.plot([], [], label=f'Body {i}')[0] for i in range(3)]
ax.legend()
def update(frame):
stop = frame*30+1
for i, line in enumerate(lines):
line.set_data(positions[max(0, stop-250):stop, i, 0], positions[max(0, stop-250):stop, i, 1])
return lines
ani = animation.FuncAnimation(fig, update, frames=60, interval=50)
plt.close(fig)
display(HTML(ani.to_jshtml()))
Reading the comparison
Energy error and trajectory error answer different questions. A symplectic method can have bounded energy error while accumulating phase error. A higher-order method can match the short-term trajectory closely yet have long-term drift. Adaptive RK45 is a numerical reference, not an exact solution. Halving the step helps expose integration error, but it cannot remove the modelling error from softening or the omission of other forces.
A three-dimensional example
The preceding initial conditions also permit a non-planar three-body orbit. Small numerical errors may grow rapidly in interacting systems, so conservation diagnostics should accompany an attractive trajectory plot.
Show or hide code
# The system we are going to look at
chaotic = SolarSystem("chaotic", 3)
chaotic.r[0] = np.array( [ -1/np.sqrt(2), 1/np.sqrt(2), 0])
chaotic.r[1] = np.array( [ 0, -1/np.sqrt(2), 1/np.sqrt(2)])
chaotic.r[2] = np.array( [ 1/np.sqrt(2), 0, -1/np.sqrt(2)])
chaotic.p[0] = np.array( [ +0.00, -0.1,0.1])
chaotic.p[1] = np.array( [ 0.1, 0.0, -0.1])
chaotic.p[2] = np.array( [ -0.1, 0.1, 0])
num_iterations = 5000 # You can adjust this as needed
results_chaotic = run_iters_and_store(chaotic, Integrator.IVP, num_iterations, 0.001)
Show or hide code
plot_results(results_chaotic,2500)
plot_xy( results_chaotic, n_to_plot = 2500)
Show or hide code
fig = plt.figure(figsize=(6, 5)); ax = fig.add_subplot(111, projection='3d')
for i in range(3):
ax.plot(*results_chaotic['r'][:, i, :].T, label=f'Body {i}')
ax.set(xlabel='x', ylabel='y', zlabel='z', title='Three-dimensional trajectories')
ax.legend(); plt.tight_layout(); plt.show()
assert np.all(np.isfinite(results_chaotic['r']))
assert np.max(np.linalg.norm(results_chaotic['COM']-results_chaotic['COM'][0], axis=1)) < 1e-8
Exercise: check convergence
Repeat the two-body comparison with half the fixed step size. Which methods show first-, second- and fourth-order reductions in short-term trajectory error? Keep the final time fixed, and compare with the same reference force law.
Show or hide code
# Answer scaffold for Verlet: about second-order convergence for this smooth orbit.
verlet_coarse = run_iters_and_store(two_body, Integrator.VERLET, 101, 0.02)
verlet_fine = run_iters_and_store(two_body, Integrator.VERLET, 201, 0.01)
errors = [np.linalg.norm(r['r'][-1]-ref['r'][-1]) for r in [verlet_coarse, verlet_fine]]
print('Verlet error ratio (expected near 4):', errors[0]/errors[1])
assert errors[1] < errors[0]/2
Interactive experiment: compare integrators
Switch between the fixed-step integrators and vary the step size. Use the energy diagnostic alongside the orbit, and pause or scrub through the animation to compare methods at the same simulated time.
The controls use a separate, seeded synthetic example and run in your browser
in the exported HTML article. Some notebook viewers do not execute embedded
JavaScript; run scripts/check_notebooks.py --html and open the article preview
to use the controls.
Show or hide code
from pathlib import Path
import sys
interactive_project = next(p / "projects" for p in [Path.cwd(), *Path.cwd().parents]
if (p / "projects" / "src").is_dir())
sys.path.insert(0, str(interactive_project / "src"))
from notebook_interactive import display_interactive
display_interactive("orbit")
Conclusion
A plausible-looking orbit is only the first check of a numerical simulation. The two-body example lets us compare trajectories, conservation diagnostics and step-size convergence against the same softened gravitational model.
For the run shown here ($h = 0.01$, $0 \leq t \leq 20$), explicit Euler has a maximum absolute energy error of about $0.54$, compared with $0.028$ for symplectic Euler and $0.0013$ for velocity Verlet. RK4 reduces this to roughly $2.1 \times 10^{-6}$. These are results for this orbit, step size and time interval; the energy plot alone does not establish which method will give the best trajectory over a much longer run.
The convergence check provides a second piece of evidence. At the same final time, halving the step reduces Verlet's position error by about a factor of four, consistent with second-order convergence. RK4's error falls by about a factor of eighteen, close to the factor of sixteen expected in its fourth-order regime. Both comparisons use a tightly toleranced adaptive RK45 solution as a numerical reference, which still has its own integration error.
Choose the diagnostics to match the question. For short-term positions, compare trajectory errors and repeat the step-size check. For long-term orbital behaviour, extend the run and inspect energy, angular momentum and phase error; symplectic methods such as Verlet are useful candidates for that comparison. For close encounters, check that the integration resolves the faster motion and test sensitivity to the softening length. The three-body plots illustrate more complicated motion, but do not by themselves demonstrate chaos.
Finally, reducing the integration error only solves the chosen model more accurately. These examples use dimensionless units and softened point-mass gravity, without collisions, relativity or additional bodies. A realistic planetary prediction would also need physical units, initial conditions tied to an epoch, and a force model appropriate to the required accuracy.
References
SciPy solve_ivp documents adaptive integration, tolerances and sampled output. Hairer, Lubich and Wanner, Geometric Numerical Integration is a reference for symplectic integration. For real Solar System initial conditions, use JPL Horizons with explicit units and epochs; these examples do not fetch or reproduce an ephemeris.