🛰️Ordinary Differential Equations: Initial Value Problems

Runge–Kutta Methods and Adaptive Step Size

Midpoint and RK4, Butcher tableaux, embedded pairs, step-size control and solve_ivp

1980: the method behind ode45

By the 1970s Runge–Kutta methods were standard, but they had an awkward weakness: a solver would take whatever step you gave it and return an answer, with no idea whether that answer was good. Choosing the step was guesswork. Erwin Fehlberg at NASA's Marshall Space Flight Center had shown in 1969 how to compute two Runge–Kutta solutions of different orders from almost the same function evaluations; their difference estimates the error, and the step can be adjusted to keep it at a tolerance. In 1980 John Dormand and Peter Prince, at Teesside Polytechnic in England, published a better-tuned pair: seven stages giving a fifth-order solution and a fourth-order one, with the extra stage reused as the first stage of the next step. In 1997 Lawrence Shampine and Mark Reichelt built MATLAB's ODE suite around it as ode45. SciPy's solve_ivp uses it by default as RK45. Most trajectories, circuits and reaction networks computed today with a general-purpose tool pass through the Dormand–Prince coefficients.

Euler uses the slope at one point. A Runge–Kutta method samples the slope at several points inside the step — stages — and combines them so that the Taylor series of the result matches the true solution to higher order. The simplest improvement is the midpoint method: take half an Euler step, evaluate the slope there, and use that slope for the full step. Two evaluations, second order.

Derivation

Why the midpoint method is second order

The midpoint step is

y_{n+1} = y_n + h\,f\big(t_n + \tfrac h2,\ y_n + \tfrac h2f\big)

, with

f = f(t_n, y_n)

. Expand the inner evaluation to first order in its small shifts (Taylor in two variables):

f\big(t_n + \tfrac h2,\ y_n + \tfrac h2f\big) = f + \tfrac h2\big(f_t + f_yf\big) + O(h^2).

Along a solution, the chain rule gives

y'' = \frac{d}{dt}f(t, y) = f_t + f_y\,y' = f_t + f_yf

. So the bracket is

y''

, and

y_{n+1} = y_n + h\,y' + \tfrac{h^2}{2}y'' + O(h^3),

which is the Taylor series of the true solution through

h^2

. The local error is

O(h^3)

and, as in lesson 9-1, the global error one order lower:

O(h^2)

. Euler's step stops at

hy'

; sampling the slope halfway along captures the curvature term too. RK4 does the same bookkeeping up to

h^4

, which needs four stages.

k_1 = f(t_n, y_n),\ k_2 = f\big(t_n + \tfrac h2, y_n + \tfrac h2k_1\big),\ k_3 = f\big(t_n + \tfrac h2, y_n + \tfrac h2k_2\big),\ k_4 = f(t_n + h, y_n + hk_3),\quad y_{n+1} = y_n + \tfrac h6(k_1 + 2k_2 + 2k_3 + k_4)

The classical fourth-order Runge–Kutta method (Kutta, 1901). Four evaluations per step; global error O(h⁴). For y′ = y it reproduces the Taylor series of eʰ through h⁴/24 — and the weights 1, 2, 2, 1 are Simpson's rule.

Derivation

RK4 on

y' = \lambda y

, and where Simpson hides

Apply one RK4 step to

y' = \lambda y

and write

z = h\lambda

. Each stage is

\lambda

times its input:

k_1 = \lambda y_n, \quad k_2 = \lambda y_n\big(1 + \tfrac z2\big), \quad k_3 = \lambda y_n\big(1 + \tfrac z2 + \tfrac{z^2}{4}\big), \quad k_4 = \lambda y_n\big(1 + z + \tfrac{z^2}{2} + \tfrac{z^3}{4}\big).

Combine them with weights

1, 2, 2, 1

y_{n+1} = y_n + \tfrac h6(k_1 + 2k_2 + 2k_3 + k_4) = \Big(1 + z + \frac{z^2}{2} + \frac{z^3}{6} + \frac{z^4}{24}\Big)\,y_n.

The exact solution multiplies by

e^z = 1 + z + z^2/2 + z^3/6 + z^4/24 + z^5/120 + \dots

; RK4 matches it through

z^4

, so its local error is

O(h^5)

and its global error

O(h^4)

. This polynomial,

R(z)

, is RK4's amplification factor; lesson 9-3 uses it to decide when a fixed step is stable.

Simpson. If

f

depends on

t

only, then

k_2 = k_3 = f(t_n + h/2)

and the step is

\tfrac h6\big(f(t_n) + 4f(t_n + \tfrac h2) + f(t_n + h)\big)

— Simpson's rule (lesson 8-2) for

\int f\,dt

💡Butcher tableaux

Any explicit Runge–Kutta method is defined by coefficients

c_i

(where in the step stage

i

is evaluated),

a_{ij}

(how earlier stages build its input) and

b_i

(the final weights):

k_i = f(t_n + c_ih, y_n + h\sum_{j<i}a_{ij}k_j)

y_{n+1} = y_n + h\sum b_ik_i

. John Butcher's notation (1960s) writes them in a table and his theory of rooted trees gives the conditions for each order. The trade-off is unavoidable: order

p

needs at least

p

stages, and for

p \ge 5

strictly more (six stages for order 5).

📐Accuracy per evaluation

For

y' = y\cos t

y(0) = 1

(exact

e^{\sin t}

) at

t = 10

: with 100 steps, Euler (100 evaluations) is off by

9.2 \times 10^{-2}

, the midpoint method (200) by

5.8 \times 10^{-4}

, RK4 (400) by

1.6 \times 10^{-7}

. Doubling the steps divides the errors by 2, 4 and 16. For smooth problems high order wins decisively, even counting evaluations.

💡Embedded pairs and step-size control

An embedded pair computes

y_{n+1}

of order

p

and

\hat y_{n+1}

of order

p-1

from the same stages. Their difference

\text{err} = \|y_{n+1} - \hat y_{n+1}\|

estimates the local error of the lower-order one. Given a tolerance:

if $\text{err} \le \text{tol}$ , accept the step (and continue with the higher-order value — 'local extrapolation');
otherwise reject it and retry with a smaller $h$ ;
either way choose the next step as $h_{\text{new}} = 0.9\,h\,(\text{tol}/\text{err})^{1/p}$ , limited to a factor of about 0.2–5.

In solve_ivp, rtol and atol set the tolerance: the error in each component is kept below about

\text{atol} + \text{rtol}\cdot|y|

. They control the local error per step, not the global error — so always check results by tightening them.

Derivation

Where the step-size formula comes from

The estimate measures the local error of the order-

(p-1)

solution, which behaves like a constant times

h^p

\text{err} \approx C\,h^p

. The constant depends on the solution's derivatives, which change slowly from step to step, so assume it is about the same on the next step. To make the next error equal to the tolerance we want

C\,h_{\text{new}}^p = \text{tol}

. Divide by

\text{err} = Ch^p

\Big(\frac{h_{\text{new}}}{h}\Big)^p = \frac{\text{tol}}{\text{err}} \qquad\Longrightarrow\qquad h_{\text{new}} = h\,\Big(\frac{\text{tol}}{\text{err}}\Big)^{1/p}.

The factor 0.9 aims slightly below the tolerance so that the next step is rarely rejected, and the limits 0.2–5 stop one unusual estimate from changing

h

wildly. For Dormand–Prince 5(4),

p = 5

: an error 32 times too large shrinks the step by about a factor of 2.

Empirical evidence can never establish mathematical existence--nor can the mathematician's demand for existence be dismissed by the physicist as useless rigor. Only a mathematical existence proof can ensure that the mathematical description of a physical phenomenon is meaningful.

– Richard Courant, as quoted in The Parsimonious Universe (1996)

📐A Molniya orbit

A Molniya orbit — perigee 600 km, apogee 39,700 km, period 11.9 h — integrated for one revolution with RK45. At rtol = 1e-4 the spacecraft returns 462 km from its start; at 1e-6, 6.6 km; at 1e-9, 3 m, using 1,046 function evaluations. The controller takes steps of about 50 s near perigee, where the spacecraft moves at 10 km/s, and up to 765 s near apogee. A fixed step small enough for perigee would waste most of its work on the slow half of the orbit.

🛰Ballistics with drag and events

A 155 mm shell, 43 kg, fired at 400 m/s and 45°. In vacuum it would travel 16.3 km; with quadratic air drag (

C_d = 0.25

) solve_ivp gives 9.28 km, a 48.6 s flight and an impact at 249 m/s. An event function — height crossing zero while falling — stops the integration exactly at impact, located by root finding (Unit 2) on the solver's interpolant. With drag the best elevation is about 41°, not 45°. Exterior ballistics was where numerical ODE solving began on an industrial scale; the next lesson starts there.

💡Where you will use this

Orbital Mechanics: high-order adaptive Runge–Kutta (and RK8(7) pairs) for precise propagation. Rocket trajectories and Optimal Guidance: ascent and re-entry simulation with events at staging, burnout and impact. Attitude Control: simulation of controller and plant together. Chemical engineering: non-stiff reaction transients. In every case: set rtol/atol consciously, use dense_output or t_eval for output, and events for discontinuities.

📚Further reading

Ernst Hairer, Syvert Nørsett and Gerhard Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems (2nd ed., Springer, 1993), chapter II, 'Runge–Kutta and Extrapolation Methods' — order conditions, Butcher's trees and embedded pairs.
John Dormand and Peter Prince, 'A family of embedded Runge–Kutta formulae', Journal of Computational and Applied Mathematics 6(1), 19–26 (1980).
Lawrence Shampine and Mark Reichelt, 'The MATLAB ODE Suite', SIAM Journal on Scientific Computing 18(1), 1–22 (1997) — how a production solver chooses steps, estimates errors and handles events.

↓But there is a problem

Each Runge–Kutta step evaluates the right-hand side four, six or seven times — and then throws all of those evaluations away. The next step starts from scratch, although the slopes it just computed are excellent information about the solution's recent past. In 1883 the British mathematician Francis Bashforth was trying to compute the shapes of drops of liquid hanging under surface tension. The integration method he used had been devised for him by John Couch Adams — the astronomer who had predicted the existence of Neptune. Adams's idea was to remember.