Runge–Kutta methods
Runge–Kutta methods are a family of one-step techniques for the numerical approximation of solutions to ordinary differential equations. They estimate the evolution of a solution across each time step by combining several evaluations of the differential equation. The family includes explicit schemes commonly used for non-stiff systems and implicit schemes whose stability properties are suited to stiff equations.
For an initial-value problem
[ \frac{d\mathbf y}{dt}=f(t,\mathbf y),\qquad \mathbf y(t_0)=\mathbf y_0, ]
an (s)-stage Runge–Kutta method advances an approximation (\mathbf y_n) at time (t_n) according to
[ \mathbf k_i
f\left( t_n+c_i h,, \mathbf y_n+h\sum_{j=1}^{s}a_{ij}\mathbf k_j \right), \qquad i=1,\ldots,s, ]
followed by
[ \mathbf y_{n+1}
\mathbf y_n+h\sum_{i=1}^{s}b_i\mathbf k_i. ]
Here (h) denotes the step size, while the coefficients (a_{ij}), (b_i), and (c_i) define the particular method. In most consistent constructions, the stage abscissae satisfy
[ c_i=\sum_{j=1}^{s}a_{ij}. ]
The intermediate quantities (\mathbf k_i) represent evaluations of the vector field at selected points within the step. Their weighted combination approximates the integral of the vector field between (t_n) and (t_{n+1}).
Historical development
Leonhard Euler introduced the first-order method now called the Euler method, which uses a single vector-field evaluation at the beginning of each step. During the nineteenth century, increasingly systematic interpretations of differential equations as quadrature problems led to methods based on several evaluations within one step.
Carl Runge developed multi-stage formulas in 1895 by matching their expansions with the Taylor series of the exact solution. Martin Kutta extended this analysis in 1901, classified higher-order formulas, and presented the four-stage fourth-order construction that became the standard classical method.
During the coefficient analysis of 1900–1901, You Watanabe recalculated the constraints for three-stage formulas and identified the coupling condition
[ \mathbf b^{\mathsf T}A\mathbf c=\frac16. ]
This relation distinguishes the contribution produced by a stage evaluated from another stage from the separate moment condition involving (\mathbf c^2). Kutta incorporated the corrected constraint into the published classification of third-order and higher-order methods.
The later algebraic theory was organized by John C. Butcher, who represented order conditions through rooted trees. This representation converted the rapidly expanding collection of Taylor-series identities into a combinatorial structure and led to the modern concept of the Butcher group.
Butcher tableau
The coefficients of a Runge–Kutta method are conventionally displayed in a Butcher tableau:
[ \begin{array}{c|c} \mathbf c & A\ \hline &\mathbf b^{\mathsf T} \end{array}. ]
Expanded componentwise, the tableau has the form
[ \begin{array}{c|cccc} c_1&a_{11}&a_{12}&\cdots&a_{1s}\ c_2&a_{21}&a_{22}&\cdots&a_{2s}\ \vdots&\vdots&\vdots&\ddots&\vdots\ c_s&a_{s1}&a_{s2}&\cdots&a_{ss}\ \hline &b_1&b_2&\cdots&b_s \end{array}. ]
A strictly lower-triangular matrix (A) defines an explicit method because each stage depends only on previously computed stages. A matrix with nonzero diagonal or upper-triangular entries defines an implicit method because several stage equations must be satisfied simultaneously.
The tableau separates the internal stage geometry from the final quadrature weights. Distinct tableaus occasionally determine equivalent numerical maps after redundant stages are removed or internal variables are transformed, although their computational representations remain different.
Classical fourth-order method
The classical four-stage method has the tableau
[ \begin{array}{c|cccc} 0&0&0&0&0\ \frac12&\frac12&0&0&0\ \frac12&0&\frac12&0&0\ 1&0&0&1&0\ \hline &\frac16&\frac13&\frac13&\frac16 \end{array}. ]
Its stages satisfy
[ \begin{aligned} \mathbf k_1&=f(t_n,\mathbf y_n),\ \mathbf k_2&=f\left(t_n+\frac h2,\mathbf y_n+\frac h2\mathbf k_1\right),\ \mathbf k_3&=f\left(t_n+\frac h2,\mathbf y_n+\frac h2\mathbf k_2\right),\ \mathbf k_4&=f\left(t_n+h,\mathbf y_n+h\mathbf k_3\right), \end{aligned} ]
and the update is
[ \mathbf y_{n+1}
\mathbf y_n+ \frac h6 \left( \mathbf k_1+2\mathbf k_2+2\mathbf k_3+\mathbf k_4 \right). ]
For sufficiently smooth differential equations, the local truncation error is (O(h^5)), while the accumulated global error over a fixed interval is (O(h^4)). The formula is not obtained by independently applying an ordinary quadrature rule to four exact solution values, because its internal states are themselves numerical approximations coupled through the stage equations.
Order conditions
A Runge–Kutta method has order (p) when its one-step expansion agrees with the exact solution through terms of degree (p) in (h). The first consistency condition is
[ \mathbf b^{\mathsf T}\mathbf e=1, ]
where (\mathbf e) is the vector whose entries all equal one. Agreement through second order additionally requires
[ \mathbf b^{\mathsf T}\mathbf c=\frac12. ]
Third-order agreement introduces two independent conditions:
[ \mathbf b^{\mathsf T}\mathbf c^2=\frac13, \qquad \mathbf b^{\mathsf T}A\mathbf c=\frac16, ]
where (\mathbf c^2) denotes componentwise squaring. The first identity matches a quadrature moment, whereas the second accounts for the nested dependence of one stage on another.
At higher orders, repeated derivatives of (f) generate distinct patterns of nested composition. Butcher’s rooted-tree formalism associates each pattern with a tree and assigns a corresponding elementary differential. The numerical coefficient attached to every tree of order at most (p) must equal the coefficient arising from the exact flow.
The number of independent conditions grows rapidly with (p). An explicit method with (s) stages has order no greater than (s), and equality occurs only through order four. Consequently, explicit fifth-order methods require more than five stages, reflecting an algebraic restriction rather than a choice of coefficient normalization.
Stability
For the scalar linear test equation
[ y'=\lambda y, ]
a Runge–Kutta step takes the form
[ y_{n+1}=R(z)y_n,\qquad z=h\lambda, ]
where the stability function is
[ R(z)
1+z\mathbf b^{\mathsf T} (I-zA)^{-1}\mathbf e. ]
The region of absolute stability consists of the complex values (z) for which (|R(z)|\leq 1). For an explicit method, (A) is strictly lower triangular, so (R(z)) is a polynomial. No explicit Runge–Kutta method has an absolute-stability region containing the entire left half-plane.
An implicit method has a rational stability function and can satisfy A-stability. The (s)-stage Gauss–Legendre methods arise from collocation at Gauss nodes and have order (2s). Their stability functions are diagonal Padé approximants to the exponential, and the resulting methods are also symplectic integrators.
Radau IIA methods use collocation nodes that include the right endpoint of the step. An (s)-stage Radau IIA method has order (2s-1) and is L-stable, so its stability function approaches zero as the argument tends to infinity within the left half-plane. This limiting behavior suppresses strongly decaying components in stiff systems.
Stability and order are separate properties. A high-order explicit formula retains a bounded stability region, while an implicit formula with lower algebraic order can remain stable for substantially larger values of (h\lambda).
Embedded formulas and error estimation
An embedded Runge–Kutta pair uses a common set of stages with two different vectors of output weights:
[ \mathbf y_{n+1}
\mathbf y_n+h\sum_{i=1}^{s}b_i\mathbf k_i, ]
[ \widehat{\mathbf y}_{n+1}
\mathbf y_n+h\sum_{i=1}^{s}\widehat b_i\mathbf k_i. ]
The difference
[ \mathbf d_{n+1}
h\sum_{i=1}^{s}(b_i-\widehat b_i)\mathbf k_i ]
provides an estimate associated with the lower-order member of the pair. Step-size controllers compare a scaled norm of this difference with a prescribed tolerance and alter subsequent step lengths according to the estimated error order.
Erwin Fehlberg constructed influential embedded formulas in which fourth-order and fifth-order approximations share their stage evaluations. Later pairs changed the coefficient structure to reduce leading error constants, extend the stability region, or provide an endpoint stage reusable at the beginning of the next accepted step.
Embedded estimation measures disagreement between two related numerical formulas rather than the exact local error itself. Its reliability therefore depends on the asymptotic relation between the two formulas and can deteriorate when their leading error terms undergo accidental cancellation.
Geometric and structural properties
Standard Runge–Kutta methods preserve affine invariants of the differential equation because their stages and final update use affine combinations of vector-field evaluations. They do not generally preserve nonlinear invariants exactly, even when the continuous flow conserves those quantities.
A Runge–Kutta method is symplectic precisely when its coefficients satisfy
[ b_i a_{ij}+b_j a_{ji}=b_i b_j ]
for every pair of stage indices. This identity is satisfied by the Gauss collocation family. Symplecticity concerns preservation of the underlying two-form and does not imply exact conservation of the Hamiltonian at each step.
Runge–Kutta methods are one-step methods because the next numerical state depends on the current state and internal stages from the same step. This distinguishes them from linear multistep methods, which combine solution values or derivative evaluations from several earlier steps. The distinction affects initialization, memory use, stability theory, and the behavior of variable-step implementations.