Runge-Kutta methods

The formula

The problem is always the same: you have dy/dt = f(t, y) with a known starting state y(t₀) = y₀, and no closed-form solution. The classical fourth-order method — RK4, the one people mean when they say "Runge-Kutta" — advances one step of size h like this:

k₁ = f(tₙ, yₙ)
k₂ = f(tₙ + h/2, yₙ + (h/2)·k₁)
k₃ = f(tₙ + h/2, yₙ + (h/2)·k₂)
k₄ = f(tₙ + h,   yₙ + h·k₃)

y_(n+1) = yₙ + (h/6)·(k₁ + 2k₂ + 2k₃ + k₄)

Reading it line by line: k₁ is the slope where you stand. k₂ is the slope at the midpoint of a trial step taken with k₁. k₃ is the midpoint slope again, retried with the better estimate k₂. k₄ is the slope at the far end of the step. The update averages them with 1-2-2-1 weights — the same weighting Simpson's rule uses to integrate a parabola exactly.

The payoff is the error scaling. Each step is wrong by O(h⁵); the accumulated error over a fixed interval is O(h⁴) ^[1]. Halve the step and the answer gets sixteen times better.

The production solvers — Fehlberg's RKF45, Dormand-Prince, MATLAB's ode45 — add one idea: compute a fourth- and a fifth-order answer from the same slope evaluations, and use their disagreement as a free error meter ^[7]^[8]:

err   ≈ |y⁽⁵⁾ − y⁽⁴⁾|
h_new = h · (tol / err)^(1/5)

First line: the gap between the two answers estimates the error of the step you just took. Second line: grow the step where the dynamics are gentle, shrink it where they aren't. That is the entire secret of adaptive step-size control.

The Dormand-Prince pair inside ode45 uses seven stages, but the last slope of an accepted step is the first slope of the next one — the "first same as last" trick — so it pays for six evaluations per step in steady running ^[8].

Where you meet it

Every ode45 call in a trajectory or 6-DOF simulation. Propagating a missile, a reentry body, or a spacecraft state vector means integrating equations of motion with no analytical solution. ode45 is the Dormand-Prince 5(4) pair under the hood ^[8]^[9], and when two sims disagree at a review board, the solver tolerances are one of the first places the argument goes.
Fixed-step RK4 in hardware-in-the-loop benches. Simulink's ode4 solver is classical RK4. Real-time HIL rigs can't tolerate an adaptive stepper deciding to take 400 substeps mid-frame, so the plant model runs fixed-step RK4 at the frame rate — and somebody has to verify that step size is small enough by checking against an adaptive run offline.
Thermal and battery transient models on the test stand. Predicting how long a test article takes to soak to temperature, or how a battery model responds to a load profile, is an initial-value problem — the same integrator, smaller stakes.
Monte Carlo dispersion runs. A thousand-case dispersion study is a thousand ode45 calls, and solver tolerance becomes a cost knob: loosen it and the run finishes overnight, but check first that the miss-distance statistics don't move when you tighten it back.

How it works

RK4 buys four orders of accuracy with four evaluations of f per step because the staged samples are arranged to cancel Taylor-series error terms through h⁴. That's an unusually good trade, and it isn't a coincidence that RK4 became the default: past order four, the price goes up — an explicit method of order five needs at least six stages, not five ^[1]. Order four is where accuracy per function evaluation peaks for simple methods, which is why a 1901 formula still ships in modern flight software.

The numbers make the case. Take dy/dt = y, y(0) = 1, exact answer e ≈ 2.718282 at t = 1. One RK4 step of h = 0.1 gives 1.10517083 against the exact 1.10517092 — error under 1×10⁻⁷ in a single step. Integrating to t = 1 at h = 0.1 costs 40 slope evaluations and lands within 2.1×10⁻⁶ of e. Forward Euler given the same 40 evaluations misses by 0.033 — about 16,000 times worse for identical cost. Halving the RK4 step cuts the error by a factor of 15.3, right on the order-four prediction of 16.

Now the fine print, which is where programs lose time.

The tolerance is not an accuracy guarantee. An adaptive solver controls the local error of each step against tol. Global error accumulates step over step and can be far larger, especially on long integrations of unstable or chaotic dynamics. The honest check is to rerun at a tolerance ten times tighter and see what moves. If the answer to a review-board question changes at RelTol = 1e-8, you didn't have an answer at 1e-6.

Explicit RK methods choke on stiff systems. The stability region of any explicit Runge-Kutta method is bounded ^[1]; for RK4 on a decaying mode dy/dt = λy, the step must satisfy h·|λ| ≲ 2.785 on the negative real axis or the numerical solution blows up while the true solution decays quietly to zero. A system with one fast time constant and one slow one — a stiff thermal network, a circuit model, chemistry — forces the step to track the fastest mode even after it has died out. The symptom is unmistakable: ode45 crawls, taking absurdly small steps on a solution that looks smooth. That is not a tuning problem. Switch to an implicit solver (ode15s and kin) and get on with the day.

Discontinuities break the error estimate. Thruster on/off, staging, contact, a relay closing — a right-hand side that jumps mid-step violates the smoothness the whole Taylor argument rests on. The stepper will still produce numbers; they'll just be less accurate than the error estimate claims, after a burst of rejected steps. The correct move is event detection: integrate exactly to the switching instant, change the model, restart.

Long propagations drift. RK4 is very slightly dissipative. Integrate an undamped oscillator (ÿ = −y) at h = 0.1 for 1,000 time units — about 159 full periods — and the system's energy falls by 0.014%. Harmless in a ten-second flight sim; not harmless when propagating an orbit for years, where the fake decay masquerades as drag. Orbit-lifetime work uses symplectic integrators or high-order methods at tight tolerance for exactly this reason.

The classic mistake, seen on real benches: a fixed-step RK4 loop with the step size chosen by folklore and no error check anywhere. RK4 with too coarse a step fails smoothly — the trajectory looks plausible and is simply wrong. Always run the step-halving check once. It costs a minute and it's the whole difference between a simulation and a cartoon.

History

Carl Runge was a Göttingen physicist-mathematician who spent years measuring atomic spectral lines, work that left him with a working scientist's need for numerical answers rather than existence proofs ^[4]. In 1895 he published "Über die numerische Auflösung von Differentialgleichungen" in Mathematische Annalen — the first method to evaluate the slope at trial points inside the step instead of only at its ends ^[2]^[1]. In 1904 Felix Klein brought him to the first chair of applied mathematics at Göttingen, a discipline Runge effectively imported into Germany ^[4].

Wilhelm Kutta generalized the idea in his 1901 paper "Beitrag zur näherungsweisen Integration totaler Differentialgleichungen," which laid out the general multi-stage framework and contains the classical fourth-order scheme engineers now call RK4 ^[3]^[10]. Aerospace readers already know Kutta from the other half of his career: the Kutta-Joukowski theorem tying airfoil lift to circulation is his too ^[3]^[5]. The same man supplied both the lift model and the integrator that propagates it.

The methods stayed hand-computation tools until digital computers made step-size selection the practical bottleneck — and the fix came from Huntsville. Erwin Fehlberg, a mathematician at NASA's Marshall Space Flight Center, published NASA TR R-315 in 1969, introducing embedded pairs: a fourth- and fifth-order formula sharing six slope evaluations, with the difference driving automatic step control ^[6]^[7]. RKF45 became the workhorse of trajectory codes. In 1980, J.R. Dormand and P.J. Prince published a refined embedded 5(4) pair that propagates the fifth-order solution and squeezes out more accuracy per evaluation ^[8]. When Shampine and Reichelt built the MATLAB ODE suite in the 1990s, they made the Dormand-Prince pair the engine of ode45 ^[9]^[11] — which is how a formula with an 1895 birthday and a Huntsville middle chapter ended up as the default button several hundred thousand engineers press every day.

Related tools

/tools/projectile-range — the no-drag closed form; add drag and you're integrating numerically, with this as the sanity check
/tools/drag-force — the right-hand side f(t, y) of most trajectory sims
/tools/terminal-velocity — the equilibrium your falling-body integration should converge to
/tools/orbital-velocity-period — a circular orbit is the standard truth case for checking a propagator
/tools/tsiolkovsky-delta-v — the closed-form bookend on the burns your simulation integrates through
/tools/rc-charge-time — a first-order ODE with an exact answer, the cheapest possible integrator test case
/tools/vibration-natural-freq — the highest natural frequency in your model sets the step size a fixed-step RK4 can get away with

Sources

Written by HE in our own words from the cited sources — engineering judgment included, your stamp still required. All entries →