Runge–Kutta for Space Mechanics

Why this matters

Numerical error can look like real orbital behaviour

In an ideal two-body orbit, the spacecraft is only acted on by central gravity. There is no drag, no thrust, and no perturbation. That means the specific mechanical energy should remain constant. If a simulation shows the orbit spiraling outward or decaying inward, the first question should be: is this physics, or is this numerical error?

This page uses Runge–Kutta methods as a space-mechanics diagnostic tool. The goal is not to memorize solver formulas. The goal is to learn how propagation accuracy, timestep size, and conservation checks affect trajectory prediction.

ideal two-body orbit → energy should be conserved
numerical propagation → energy must be checked

Core message: Higher-order integration can improve orbit prediction, but long-term accuracy depends on timestep, conservation checks, and whether the numerical method preserves orbital physics.

Problem roadmap

Orbit propagation as a numerical experiment

A · Two-body orbit propagation
Define the ideal gravitational model and compute initial orbit quantities.

B · Orbit paths
Compare Euler, RK4, RK7, and a reference-like solution visually.

C · Energy conservation
Use specific mechanical energy as a physics diagnostic.

D · Trajectory prediction
Measure final position error after several revolutions.

E · Step-size sensitivity
Sweep timestep values and see how error grows.

Final link
Treat timestep selection as an engineering decision, not a cosmetic setting.

Part A — Two-body orbit propagation

The model: central gravity only

The two-body model assumes a spacecraft moving under Earth’s central gravitational acceleration. In vector form, the acceleration points toward Earth’s centre.

\[\ddot{\mathbf r} = -\frac{\mu}{\lVert \mathbf r \rVert^3}\mathbf r\]

For a clean educational setup, the spacecraft starts at perigee. The initial velocity is computed from the chosen altitude and eccentricity using the vis-viva equation. When eccentricity is zero, the orbit is circular. When eccentricity is above zero, the starting point is perigee of an elliptical orbit.

Initial perigee altitude [km]

Eccentricity e

Number of orbits

Timestep Δt [s]

Perigee radius—

Semi-major axis—

Initial speed—

Orbital period—

What to check: In this ideal model, energy should remain nearly constant. If it does not, the issue is likely numerical rather than physical.

Part B — Euler vs RK4 vs RK7 orbit paths

Seeing orbit distortion directly

Orbit plots make numerical error visually obvious. Euler can spiral outward or distort because it does not preserve the orbital geometry well. RK4 and RK7 stay closer to the reference trajectory, especially when the timestep is reasonable.

Perigee altitude [km]

Eccentricity e

Timestep Δt [s]

Number of orbits

Euler final position error—

RK4 final position error—

RK7 final position error—

Euler radius drift—

Visual insight: A distorted orbit path is not automatically a new physical effect. It may be the numerical method failing to preserve the orbit.

Part C — Long-term energy conservation

Energy drift as a physics diagnostic

Specific mechanical energy is a powerful diagnostic because the ideal two-body problem should conserve it. This plot shows how much the numerical energy changes over time.

\[\epsilon = \frac{v^2}{2} - \frac{\mu}{r}\]

Perigee altitude [km]

Eccentricity e

Timestep Δt [s]

Number of orbits

Euler energy drift—

RK4 energy drift—

RK7 energy drift—

Best conservation—

Core insight: Energy drift is a diagnostic. It tells you whether your numerical solution is respecting the physics.

Part D — Trajectory prediction error

Small per-step errors compound

A propagator can look acceptable after one orbit but become poor after many revolutions. This section compares final-state prediction error against a reference-like RK7 solution with a very small timestep.

Perigee altitude [km]

Eccentricity e

Timestep Δt [s]

Maximum orbit count

Euler final error—

RK4 final error—

RK7 final error—

Error growth—

Insight: Orbit propagation errors compound. A method that looks acceptable for one orbit may be unacceptable after many revolutions.

Part E — Step-size sensitivity map

Timestep is an engineering choice

This is the strongest diagnostic section. Instead of testing one timestep, sweep a range of timesteps and compare final position error and energy drift. This shows how quickly a numerical method becomes unreliable as the timestep increases.

Perigee altitude [km]

Eccentricity e

Minimum Δt [s]

Maximum Δt [s]

Number of tested step sizes

Number of orbits

Worst Euler error—

Worst RK4 error—

Worst RK7 error—

Recommended action—

Insight: Step size is an engineering choice, not just a numerical setting.

Python / MATLAB snippets

Code patterns for orbit propagation checks

These snippets show the core checks used in this page: two-body acceleration, Runge–Kutta stepping, and specific mechanical energy monitoring.

% Two-body acceleration in 2D
mu = 398600.4418; % km^3/s^2

r = x(1:2);
v = x(3:4);
rnorm = norm(r);

a = -mu * r / rnorm^3;
xdot = [v; a];

% Fixed-step RK4 update
k1 = f(t, x);
k2 = f(t + dt/2, x + dt*k1/2);
k3 = f(t + dt/2, x + dt*k2/2);
k4 = f(t + dt,   x + dt*k3);

x_next = x + dt*(k1 + 2*k2 + 2*k3 + k4)/6;

% Specific mechanical energy diagnostic
rnorm = norm(r);
vnorm = norm(v);
energy = 0.5*vnorm^2 - mu/rnorm;

% For an ideal two-body orbit, energy should remain nearly constant.
% If energy drifts strongly, check timestep and integration method.

Engineering takeaway

Orbit propagation must be checked against invariants

1 · Orbit path
A visual orbit plot can reveal spiraling, distortion, or phase error.

2 · Energy conservation
Specific mechanical energy is a diagnostic for ideal two-body propagation.

3 · Prediction error
Small per-step errors can become large final-state errors over many orbits.

4 · Step-size sensitivity
A timestep that works for one orbit may fail for long propagation.

5 · RK4 and RK7
Higher-order methods help, but only when timestep choice remains reasonable.

6 · Engineering judgement
The solver is part of the simulation assumption set.

Final insight: For orbit propagation, accuracy is not only about choosing RK4 or RK7. It is about checking whether the numerical trajectory preserves the physical invariants of the problem.

← Back to Problem Bank