4 Runge-Kutta Methods
The Euler method uses a single evaluation of \bm{f} per step and is limited to first order; the Taylor methods of the previous chapter reach any order, but at the price of the derivatives of \bm{f}. Runge–Kutta methods occupy the middle ground: they obtain high order using only evaluations of \bm{f}, placed at several intermediate points inside the step (see @butcher2008numerical, @hairer1993solving).
The presentation here goes from the general to the particular. The general method is defined first, with no restriction on its coefficients; the familiar explicit schemes then appear as the special case in which the coefficient matrix is strictly lower triangular, a structure that reduces the computation of a step to a sequence of function evaluations. Taking the opposite route — explicit first, implicit as an afterthought — hides the fact that the order conditions, the Butcher tableau and the whole theory are the same in both cases, and that only the cost of one step changes.
Before introducing Runge–Kutta schemes, recall the standard initial-value problem for an explicit ODE,
\begin{cases} \bm{x}'(t)=\bm{f}(\bm{x}(t),t),\\ \bm{x}(0)=\bm{x}_0. \end{cases}
More generally, an ODE may be given in fully implicit form,
\begin{cases} \bm{F}(\bm{x}(t),\bm{x}'(t),t)=\bm{0},\\ \bm{x}(0)=\bm{x}_0, \end{cases}
where the equation need not be solved explicitly for \bm{x}'. The explicit form is the special case \bm{F}(\bm{x},\bm{v},t)=\bm{v}-\bm{f}(\bm{x},t). To describe an ODE rather than a DAE, \partial\bm{F}/\partial\bm{v} must be locally nonsingular along the solution, so that the implicit-function theorem locally determines \bm{v}=\bm{x}'. Runge–Kutta methods for this general form will be written below using the same stage increments, with the stage derivatives entering \bm{F} directly.
4.1 The General Runge–Kutta Method
An s-stage Runge–Kutta method advances the solution of \bm{y}'=\bm{f}(\bm{y},t) from t_k to t_{k+1}=t_k+\Delta t through s stages \bm{K}_1,\ldots,\bm{K}_s:
\begin{aligned} \bm{z}_{k+1} &= \bm{z}_k + \sum_{i=1}^{s}b_i\bm{K}_i \\ &\left\{\begin{aligned} \bm{K}_1 &= \Delta t\bm{f}\Big(\bm{z}_k + \displaystyle\sum_{j=1}^{s}a_{1j}\bm{K}_j,t_k+\Delta t c_1\Big) \\[0.25em] \bm{K}_2 &= \Delta t\bm{f}\Big(\bm{z}_k + \displaystyle\sum_{j=1}^{s}a_{2j}\bm{K}_j,t_k+\Delta t c_2\Big) \\[0.25em] \vdots \\[0.25em] \bm{K}_s &= \Delta t\bm{f}\Big(\bm{z}_k + \displaystyle\sum_{j=1}^{s}a_{sj}\bm{K}_j,t_k+\Delta t c_s\Big) \end{aligned} \right. \end{aligned} \tag{4.1}
The method is completely described by the s^2+2s coefficients a_{ij}, b_i, c_i, collected in the Butcher tableau
\def\arraystretch{1.5} \begin{array} {c|cccc} c_1 & a_{11} & a_{12} & \ldots & a_{1s} \\ c_2 & a_{21} & a_{22} & \ldots & a_{2s} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ c_s & a_{s1} & a_{s2} & \ldots & a_{ss} \\ \hline & b_1 & b_2 & \ldots & b_s \end{array} \qquad\text{or}\qquad \begin{array}{c|c} \bm{c} & \bm{A}\\ \hline & \bm{b}^T \end{array}
The quantities \bm{K}_i/\Delta t are approximations of \bm{f} at the internal times t_k+c_i\Delta t, and the final update is a weighted average of them. The abscissae are not free: they must satisfy the stage-consistency relation
c_i=\sum_{j=1}^{s}a_{ij}, \qquad i=1,\ldots,s, \qquad\text{that is}\qquad \bm{A}\bm{1}=\bm{c},
so that the i-th stage really is an approximation at time t_k+c_i\Delta t; when 0\le c_i\le1 the stage times lie inside the current step. This relation is assumed throughout and is revisited in 2.
In Equation 4.1 the stages appear on both sides: each \bm{K}_i depends on all the others. For a problem of dimension d, one step therefore requires the solution of a nonlinear system of s\cdot d equations, not the evaluation of a formula. Whether that system is expensive, cheap, or absent altogether depends entirely on the shape of \bm{A} — which is the subject of the next section.
Runge–Kutta stages for a fully implicit ODE
For the fully implicit problem, the same tableau and stage increments give
\begin{aligned} \bm{z}_{k+1} &= \bm{z}_k+\sum_{i=1}^s b_i\bm{K}_i,\\ \bm{F}\!\left(\bm{z}_k+\sum_{j=1}^s a_{ij}\bm{K}_j, \frac{\bm{K}_i}{\Delta t}, t_k+c_i\Delta t\right)&=\bm{0},\qquad i=1,\ldots,s. \end{aligned} \tag{4.2}
Here \bm{K}_i/\Delta t approximates the derivative at stage i. Unlike the explicit-ODE formula, even a strictly lower-triangular tableau does not in general give an explicit algorithm: once the earlier stages are known, the stage equation still has to be solved for \bm{K}_i because the derivative is an argument of \bm{F}. A lower-triangular (DIRK) tableau makes these nonlinear stage solves sequential; a full tableau couples all stages and requires a coupled nonlinear solve. Thus fully implicit ODEs generally call for a diagonally implicit or fully implicit RK procedure, rather than an explicit RK procedure in the usual computational sense.
4.2 Constructing Runge-Kutta Methods
Everything needed to build a Runge–Kutta method and to use it in practice: the shape of the coefficient matrix, how a step is actually solved, the conditions that fix the order, and the embedded pairs that drive step-size selection.
The Structure of \bm{A}
Three shapes of the coefficient matrix are used in practice, and they correspond to three very different costs per step.
Fully implicit methods
\bm{A} is a full matrix: all stages are coupled and the sd\times sd system must be solved as a whole. This is the most expensive option, and it buys the most: an s-stage method can reach order 2s (the Gauss–Legendre family), and these methods have the best stability properties available — several of them are A-stable or L-stable for every s. They are the methods of choice when the problem is very stiff, when high accuracy is required at large step sizes, or when a structural property such as symplecticity must be preserved.
Diagonally implicit methods (DIRK)
\bm{A} is lower triangular, a_{ij}=0 for j>i:
\def\arraystretch{1.5} \begin{array}{c|cccc} c_1 & a_{11} & 0 & \ldots & 0 \\ c_2 & a_{21} & a_{22} & \ldots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ c_s & a_{s1} & a_{s2} & \ldots & a_{ss} \\ \hline & b_1 & b_2 & \ldots & b_s \end{array}
The i-th stage then involves only \bm{K}_1,\ldots,\bm{K}_i, so the stages are computed one at a time:
\bm{K}_i = \Delta t\bm{f}\Big(\bm{z}_k+\underbrace{\sum_{j<i}a_{ij}\bm{K}_j}_{\text{already known}} +a_{ii}\bm{K}_i,t_k+\Delta t c_i\Big), \qquad i=1,\ldots,s.
Each stage is still implicit — \bm{K}_i appears on both sides through the diagonal entry a_{ii} — but the system to be solved has dimension d, not sd. Since the cost of a direct solve grows like the cube of the dimension, replacing one sd\times sd factorization by s factorizations of size d\times d saves a factor of about s^2.
Two refinements matter in practice.
- SDIRK (singly diagonally implicit): all the diagonal entries are equal, a_{11}=\cdots=a_{ss}=\gamma. Then every stage requires the same Newton matrix \bm{I}-\Delta t\gamma\bm{J}_{\bm{f}}, which can be factorized once per step — or even reused across several steps — instead of once per stage.
- ESDIRK (explicit first stage): a_{11}=0, so the first stage is explicit and free, while the remaining ones are singly diagonally implicit. This also makes the method FSAL (Section 4.2.5), saving one further evaluation.
DIRK and SDIRK methods are the standard compromise in production stiff solvers: they keep A-stability and L-stability, which an explicit method can never have, while paying only a few d\times d solves per step. Their price is the order barrier — a DIRK method of order p generally needs more stages than a fully implicit one — and a slightly larger error constant.
Explicit methods
\bm{A} is strictly lower triangular, a_{ij}=0 for j\ge i — the diagonal vanishes as well:
\begin{array} {c|ccccc} 0 & 0 & 0 & 0 & 0 & 0 \\ c_2 & a_{21} & 0 & 0 & 0 & 0 \\ c_3 & a_{31} & a_{32} & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & 0 & 0 \\ c_s & a_{s1} & a_{s2} & \ldots & a_{s,s-1} & 0 \\ \hline & b_1 & b_2 & \ldots & b_{s-1} & b_s \end{array}
Now the right-hand side of the i-th stage involves only stages already computed, and Equation 4.1 collapses from a nonlinear system to an explicit recipe:
\begin{aligned} \bm{z}_{k+1} &= \bm{z}_k + \sum_{i=1}^{s} b_i \bm{K}_i, \\ &\left\{\begin{aligned} \bm{K}_1 &= \Delta t \bm{f}\Big(\bm{z}_k,t_k\Big), \\[0.5em] \bm{K}_2 &= \Delta t \bm{f}\Big(\bm{z}_k + a_{21} \bm{K}_1,t_k + \Delta t c_2\Big), \\[0.5em] \bm{K}_3 &= \Delta t \bm{f}\Big(\bm{z}_k + a_{31} \bm{K}_1 + a_{32} \bm{K}_2,t_k + \Delta t c_3\Big), \\[0.5em] \vdots & \\[0.5em] \bm{K}_s &= \Delta t \bm{f}\Big(\bm{z}_k + \sum_{j=1}^{s-1} a_{sj} \bm{K}_j,t_k + \Delta t c_s\Big). \end{aligned} \right. \end{aligned} \tag{4.3}
One step costs exactly s evaluations of \bm{f}, with no linear algebra, no Jacobian and no iteration. This is why explicit methods dominate non-stiff computation. The simplification is not free: as shown in Section 4.2.4.1, the strict triangularity makes \bm{A} nilpotent and caps the attainable order at s, and the stability function becomes a polynomial, so the absolute-stability region is necessarily bounded.
Summary of the three structures
| structure | \bm{A} | cost of one step | max order | A-stable? |
|---|---|---|---|---|
| explicit | strictly lower triangular | s evaluations of \bm{f} | s | never |
| DIRK / SDIRK | lower triangular | s systems d\times d | below 2s | possible |
| fully implicit | full | one system sd\times sd | 2s | possible |
Everything that follows — order conditions, embedded pairs, stability functions — applies verbatim to all three. Only the work needed to evaluate one step changes.
Solving the Stage Equations
For a fully implicit or diagonally implicit method, Equation 4.1 must actually be solved at every step. The two candidates are the methods of the chapter on nonlinear equations: fixed-point iteration and Newton’s method.
Fixed-point iteration, and why it is not enough
The simplest option is to read Equation 4.1 directly as a fixed-point iteration,
\bm{K}_i^{(m+1)} = \Delta t\bm{f}\Bigl(\bm{z}_k+\sum_{j=1}^s a_{ij}\bm{K}_j^{(m)},\;t_k+c_i\Delta t\Bigr), \qquad m=0,1,2,\ldots
started for instance from \bm{K}_i^{(0)}=\Delta t\bm{f}(\bm{z}_k,t_k). If \bm{f} is Lipschitz in \bm{y} with constant L, the map is a contraction as soon as
\boxed{\; \Delta t\,L\max_{i}\sum_{j=1}^{s}|a_{ij}|<1, \;} \tag{4.4}
and the iteration then converges linearly with a rate equal to the left-hand side of Equation 4.4.
Condition Equation 4.4 imposes \Delta t\lesssim1/L. But implicit methods exist precisely for stiff problems, where L is huge: the contraction condition would reintroduce exactly the step-size restriction that implicitness was meant to remove. On a stiff problem a fixed-point solver cancels the benefit of the implicit formula, and Newton’s method becomes mandatory.
Solving the coupled stage equations with Newton’s method
For an implicit Runge–Kutta method the stages are coupled. It is convenient to work with stage states
\bm{Y}_i = \bm{z}_n +\Delta t\sum_{j=1}^s a_{ij}\bm{f}(\bm{Y}_j,t_n+c_j\Delta t), \qquad i=1,\ldots,s.
If each state has dimension d, stack all stages into the vector
\bm{Y} = \begin{pmatrix} \bm{Y}_1\\ \vdots\\ \bm{Y}_s \end{pmatrix} \in\mathbb R^{sd}.
Define the nonlinear residual blocks
\bm{G}_i(\bm{Y}) = \bm{Y}_i-\bm{z}_n -\Delta t\sum_{j=1}^s a_{ij} \bm{f}(\bm{Y}_j,t_n+c_j\Delta t).
The stage problem is the nonlinear system
\boxed{\bm{G}(\bm{Y})=\bm{0}},
of dimension sd. Newton’s method computes \bm{S}^{(m)} from
\bm{J}_{\bm{G}}(\bm{Y}^{(m)})\bm{S}^{(m)} =-\bm{G}(\bm{Y}^{(m)}),
and sets
\bm{Y}^{(m+1)}=\bm{Y}^{(m)}+\bm{S}^{(m)}.
The (i,j) block of the Newton matrix is
\boxed{ \left[\bm{J}_{\bm{G}}\right]_{ij} = \delta_{ij}\bm{I} -\Delta t a_{ij}\bm{J}_{\bm{f}}(\bm{Y}_j,t_n+c_j\Delta t) }
with d\times d blocks. Hence a fully implicit s-stage RK method generally requires the solution of a coupled linear system of dimension sd at every Newton iteration.
For large systems, forming and factorizing this full matrix can dominate the cost. Practical solvers may therefore use simplified Newton iterations, reuse Jacobians and factorizations, exploit Kronecker-product structure, or use iterative linear solvers with preconditioning. The structural alternative is to give up the full \bm{A} altogether and use a DIRK method (Section 4.2.1), whose Newton matrix is block lower triangular.
Unlike the fixed-point iteration, Newton’s method converges quadratically and requires no condition of the form Equation 4.4: the step size \Delta t stays free, which is the whole reason for choosing an implicit method.
- Simplified Newton. Freeze \bm{J}_{\bm{f}} at (t_k,\bm{z}_k) and reuse the same factorization across all iterations, and often across several consecutive steps. Quadratic convergence is lost, but so is the dominant cost, the \mathcal O\bigl((sd)^3\bigr) factorization.
- DIRK and SDIRK methods. With \bm{A} lower triangular the Newton matrix is block lower triangular and the stages are solved one at a time: s systems of size d\times d instead of one of size sd\times sd (Section 4.2.1). If in addition all the diagonal entries are equal, the same factorization serves every stage.
- Predictor. A good starting guess \bm{K}^{(0)}, extrapolated from the previous step, usually reduces the iteration count to two or three.
Order Conditions: the Two-Stage Case
The conditions listed in the next section are usually presented as a table to be looked up. They are, however, nothing more than the result of matching two Taylor expansions, and carrying out the match once — for the simplest nontrivial case — explains both their form and their number.
Two-Variable Taylor Expansion with Remainder
Everything rests on expanding f around (y,t) in both arguments. It is worth stating the expansion with a quantitative remainder, because the order conditions and the error constants both come from there.
Proposition 4.1 (Taylor expansion with remainder) Let f\in C^{1} in a neighbourhood of (y,t). Then
f(y+\beta,t+\alpha) = f(y,t)+\beta f_y(y,t)+\alpha f_t(y,t)+r_1(\alpha,\beta),
and if the gradient (f_y,f_t) is Lipschitz with constant M_1 in that neighbourhood,
\boxed{\; \bigl|r_1(\alpha,\beta)\bigr|\le\frac{M_1}{2}\bigl(\alpha^2+\beta^2\bigr). \;} \tag{4.5}
If f\in C^{2} the second-order expansion holds,
f(y+\beta,t+\alpha) = f+\beta f_y+\alpha f_t +\frac12\bigl(\beta^2f_{yy}+2\alpha\beta f_{ty}+\alpha^2f_{tt}\bigr) +r_2(\alpha,\beta),
with |r_2|=\mathcal O\bigl((\alpha^2+\beta^2)^{3/2}\bigr), and |r_2|\le\frac{M_2}{6}(\alpha^2+\beta^2)^{3/2} if the Hessian is Lipschitz with constant M_2.
In a Runge–Kutta method the expansion is always applied with \alpha=c_i\Delta t and \beta=\mathcal O(\Delta t): every increment is proportional to \Delta t, so the remainder Equation 4.5 is \mathcal O(\Delta t^2) and, at second order, \mathcal O(\Delta t^3). This is what allows the expansions to be truncated at the desired order and the terms to be collected in powers of \Delta t.
A method of order p requires f\in C^{p}: the order conditions involve the partial derivatives of f up to order p-1, and the exact solution must be C^{p+1}. If f is merely Lipschitz the method still converges — the convergence theorem Theorem 2.1 uses only the Lipschitz property — but with order 1. High order is paid for in regularity, not only in stages.
The Local Truncation Error of a Two-Stage Method
Consider the scalar problem y'=f(y,t) and the two-stage increment
\Phi(t,y,\Delta t)=y+A\Delta t f(y,t)+B\Delta t f\bigl(y+D\Delta t f(y,t),\;t+C\Delta t\bigr),
with four parameters A,B,C,D to be determined. The local truncation error is the difference between the exact solution and one step of the method started on it,
\tau(\Delta t)=y(t+\Delta t)-\Phi\bigl(t,y(t),\Delta t\bigr).
For the exact solution the differential equation supplies the first two derivatives,
y(t+\Delta t) = y(t)+\underbrace{f}_{y'}\Delta t +\underbrace{\bigl(f_t+f_yf\bigr)}_{y''}\frac{\Delta t^2}{2} +\mathcal O(\Delta t^3),
all quantities being evaluated at (y(t),t). For the method, Proposition 4.1 with \alpha=C\Delta t and \beta=D\Delta t f gives
f\bigl(y+D\Delta t f,t+C\Delta t\bigr) = f+C\Delta t f_t+D\Delta t f_yf+\mathcal O(\Delta t^2).
Subtracting and collecting powers of \Delta t,
\boxed{ \begin{aligned} \tau(\Delta t) &= \Delta t\,f\,(1-A-B) \\ &\quad +\frac{\Delta t^2}{2}\,f_t\,(1-2BC) +\frac{\Delta t^2}{2}\,f_yf\,(1-2BD) +\mathcal O(\Delta t^3). \end{aligned} } \tag{4.6}
Reading the Conditions Off the Expansion
The three quantities f, f_t and f_yf are independent: no relation holds between them for a general right-hand side, so each of their coefficients must vanish separately. This is the mechanism that generates the order conditions, and the reason their number grows with the order — at order p one coefficient appears for each elementary differential of order p, that is, for each rooted tree with p nodes.
Requiring \tau(\Delta t)=\mathcal O(\Delta t^3), hence a globally second-order method, gives
\left\{ \begin{aligned} 1-A-B &= 0,\\ 1-2BC &= 0,\\ 1-2BD &= 0, \end{aligned} \right. \qquad\Longrightarrow\qquad A=1-B,\quad C=D=\frac1{2B}.
Three conditions and four unknowns leave a one-parameter family of second-order two-stage methods,
y_{k+1} = y_k+(1-B)\Delta t f(y_k,t_k) +B\Delta t f\!\left(y_k+\frac{\Delta t}{2B}f(y_k,t_k),\,t_k+\frac{\Delta t}{2B}\right),
which in the standard notation b_1=1-B, b_2=B, c_1=0, c_2=a_{21}=\frac1{2B} reads
\begin{aligned} y_{k+1} &= y_k+b_1\bm{K}_1+b_2\bm{K}_2,\\ \bm{K}_1 &= \Delta t f(y_k,t_k),\\ \bm{K}_2 &= \Delta t f\bigl(y_k+a_{21}\bm{K}_1,\;t_k+c_2\Delta t\bigr). \end{aligned}
The familiar members of the family correspond to particular choices of the free parameter:
\begin{array}{c|c|l} B & c_2 & \text{method}\\ \hline 1/2 & 1 & \text{Heun (explicit trapezoidal)}\\ 1 & 1/2 & \text{explicit midpoint}\\ 3/4 & 2/3 & \text{Ralston (minimal error constant)} \end{array}
The remaining freedom is not accuracy but the size of the leading error term: Ralston’s choice minimises the coefficient of \Delta t^3 in \tau, which is why it appears as the second-order member of embedded pairs.
Beyond Two Stages
Repeating the construction with s stages produces s^2+2s coefficients a_{ij},b_i,c_i and, at each order, one algebraic equation per elementary differential. The bookkeeping becomes unmanageable by hand beyond order four, and the systematic tool is Butcher’s rooted-tree formalism. The resulting conditions in vector form are collected next; Section 7.1 applies them to concrete tableaux.
Order Conditions: the General Case
The derivation of 1 produced three conditions for a two-stage method of order two. The same procedure, carried out with s stages and pushed to order p, produces a finite list of algebraic identities among the coefficients a_{ij}, b_i, c_i. This section collects that list.
Two facts organise the whole section.
- The conditions are cumulative. A method has order p when every condition listed for orders 1,\ldots,p holds; it has order exactly p when in addition at least one condition of order p+1 fails. Nothing can be skipped: all the conditions of a level must hold before the level is attained.
- Each condition corresponds to one elementary differential of the right-hand side — one way of combining f and its derivatives — and therefore to one rooted tree. This is why the number of conditions grows the way it does, and it is the content of Section 4.2.4.1.
Preliminary: reduction to autonomous form
The conditions take their cleanest form for an autonomous system \bm{y}'=\bm{f}(\bm{y}). This is no loss of generality: the non-autonomous problem \bm{y}'=\bm{f}(\bm{y},t) becomes autonomous by appending the time variable,
\widetilde{\bm{y}}=\begin{pmatrix}t\\ \bm{y}\end{pmatrix}, \qquad \widetilde{\bm{f}}(\widetilde{\bm{y}}) =\begin{pmatrix}1\\ \bm{f}(\bm{y},t)\end{pmatrix}.
A Runge–Kutta method is invariant under this transformation — it produces the same numerical solution for the original and the augmented system — precisely when the stage-consistency relation below holds. For such methods the autonomous and non-autonomous order conditions coincide, and one may work with whichever form is more convenient.
Stage consistency
Before any order condition comes the relation defining the stage abscissae,
\boxed{\; {\color{magenta} c_i=\sum_{j=1}^{s}a_{ij}, \qquad i=1,\ldots,s, \qquad\text{that is}\qquad \bm{A}\bm{1}=\bm{c}}, \;} \tag{4.7}
with \bm{1}=(1,\ldots,1)^T. It states that the i-th stage is a consistent approximation of the solution at time t_k+c_i\Delta t. It is assumed throughout, it is what makes the method invariant under the autonomisation above, and a tableau that violates it has almost certainly been mistranscribed.
Order 1
In addition to Equation 4.7, one condition:
{\color{blue}\sum_{i=1}^{s}b_i=1.}
The increment must reproduce the first term \Delta t f of the Taylor expansion. Every usable method satisfies it; a tableau that does not is inconsistent and does not converge at all.
Order 2
One further condition:
{\color{DarkRed}\sum_{i=1}^{s}b_ic_i=\frac12.}
This is the weight \tfrac12 \Delta t^2 y'' of the expansion, and it is where the abscissae first enter.
Order 3
Two further conditions:
{\color{DarkGreen} \sum_{i=1}^{s}b_ic_i^2=\frac13, \qquad \sum_{i,j=1}^{s}b_ia_{ij}c_j=\frac16.}
The two correspond to the two distinct elementary differentials of order three, f_{yy}(f,f) and f_yf_yf: the first is seen by the abscissae alone, the second only through the internal coupling \bm{A}. An explicit one-stage method can never satisfy the second, and this is the structural reason why more stages are needed as the order grows.
Order 4
Four further conditions:
{\color{teal} \begin{aligned} \sum_{i=1}^{s}b_ic_i^3 &= \frac14, & \sum_{i,j=1}^{s}b_ic_ia_{ij}c_j &= \frac18, \\[0.4em] \sum_{i,j=1}^{s}b_ia_{ij}c_j^2 &= \frac1{12}, & \sum_{i,j,k=1}^{s}b_ia_{ij}a_{jk}c_k &= \frac1{24}. \end{aligned}}
Note the last one: it involves a double application of \bm{A}, so an explicit method needs at least three nonzero sub-diagonal couplings, hence at least four stages, to satisfy it.
Vector form
Writing the sums as inner products makes the conditions far easier to check by hand. With \circ the componentwise (Hadamard) product, set
\boxed{\; \bm{d}=\bm{A}\bm{c}, \qquad \bm{e}=\bm{A}^T\bm{b}, \qquad \bm{c}^{\circ2}=\bm{c}\circ\bm{c}, \qquad \bm{c}^{\circ3}=\bm{c}\circ\bm{c}\circ\bm{c}. \;} \tag{4.8}
Then the eight conditions through order four read
\color{teal} \left. \begin{aligned} &\color{DarkGreen}\underbrace{ \color{DarkRed}\overbrace{ \color{blue}\overbrace{ \underbrace{\color{magenta}\bm{A}\bm{1}=\bm{c}}_{\text{stage consistency}},\quad \color{blue}\bm{b}\cdot\bm{1}=1 }^{\textstyle\text{order }1},\quad \color{DarkRed}\bm{b}\cdot\bm{c}=\dfrac12 }^{\textstyle\text{order }2},\quad \color{DarkGreen}\bm{b}\cdot\bm{c}^{\circ2}=\dfrac13,\quad \bm{b}\cdot\bm{d}=\dfrac16 }_{\textstyle\text{order }3} \\[1.2em] &\color{teal}\quad \bm{b}\cdot\bm{c}^{\circ3}=\dfrac14,\quad (\bm{b}\circ\bm{c})\cdot\bm{d}=\dfrac18,\quad \bm{e}\cdot\bm{c}^{\circ2}=\dfrac1{12},\quad \bm{e}\cdot\bm{d}=\dfrac1{24} \end{aligned} \right\} {\color{teal}\text{order }4} \tag{4.9}
The colours mark the levels and are used consistently in this section: \color{magenta}\text{stage consistency}, \color{blue}\text{order }1, \color{DarkRed}\text{order }2, \color{DarkGreen}\text{order }3, \color{teal}\text{order }4, \color{Purple}\text{order }5.
The two auxiliary vectors do all the work: \bm{d}=\bm{A}\bm{c} carries the internal coupling seen from the stages, and \bm{e}=\bm{A}^T\bm{b} the same coupling seen from the weights. Both are computed once and reused. Section 7.1 applies this table to a dozen tableaux, and Section 7.3 turns it into a few lines of code.
Order 5, for reference
Nine further conditions, listed here so that the fourth-order claims of the tableaux below can actually be sharpened to exactly four:
{\color{Purple} \begin{aligned} \sum_i b_ic_i^4 &= \frac15, & \sum_i b_ic_i^2 d_i &= \frac1{10}, & \sum_i b_i d_i^2 &= \frac1{20}, \\[0.4em] \sum_i b_ic_i\bigl(\bm{A}\bm{c}^{\circ2}\bigr)_i &= \frac1{15}, & \sum_i b_i\bigl(\bm{A}\bm{c}^{\circ3}\bigr)_i &= \frac1{20}, & \sum_i b_ic_i\bigl(\bm{A}\bm{d}\bigr)_i &= \frac1{30}, \\[0.4em] \sum_i b_i\bigl(\bm{A}(\bm{c}\circ\bm{d})\bigr)_i &= \frac1{40}, & \sum_i b_i\bigl(\bm{A}^2\bm{c}^{\circ2}\bigr)_i &= \frac1{60}, & \sum_i b_i\bigl(\bm{A}^2\bm{d}\bigr)_i &= \frac1{120}. \end{aligned}} \tag{4.10}
For the classical RK4 tableau the very first of these already fails, \bm{b}\cdot\bm{c}^{\circ4}=\tfrac5{24}\ne\tfrac15, so its order is exactly four. All nine hold for the three-stage Radau IIA method (order five) and for the three-stage Gauss–Legendre method (order six).
Higher orders: the rooted-tree rule
Writing conditions one by one stops being practical at order five. Butcher’s observation is that the terms of the Taylor expansion of the exact solution and of the numerical one are both indexed by rooted trees, and that matching them term by term gives one condition per tree.
In bracket notation a tree is built from the single node \bullet: [\tau_1\cdots\tau_m] denotes the tree whose root has the subtrees \tau_1,\ldots,\tau_m. To each tree \tau with \rho(\tau) nodes one associates
- the elementary differential F(\tau), a combination of f and its derivatives, which appears in the expansion of the exact solution;
- the elementary weight \Phi(\tau), a polynomial in a_{ij},b_i,c_i, which appears in the expansion of the numerical solution;
- the density \gamma(\tau), an integer computed from the tree alone.
Theorem 4.1 (Butcher’s order conditions) A Runge–Kutta method has order p if and only if
\boxed{\; \Phi(\tau)=\frac{1}{\gamma(\tau)} \quad\text{for every rooted tree }\tau\text{ with }\rho(\tau)\le p. \;}
The first eight trees reproduce exactly the conditions already listed:
\begin{array}{c|c|c|c|c} \rho & \tau & F(\tau) & \Phi(\tau) & \gamma(\tau)\\ \hline 1 & \bullet & f & {\color{blue}\bm{b}\cdot\bm{1}} & 1\\[0.3em] 2 & [\bullet] & f_yf & {\color{DarkRed}\bm{b}\cdot\bm{c}} & 2\\[0.3em] 3 & [\bullet^2] & f_{yy}(f,f) & {\color{DarkGreen}\bm{b}\cdot\bm{c}^{\circ2}} & 3\\[0.3em] 3 & [[\bullet]] & f_yf_yf & {\color{DarkGreen}\bm{b}\cdot\bm{d}} & 6\\[0.3em] 4 & [\bullet^3] & f_{yyy}(f,f,f) & {\color{teal}\bm{b}\cdot\bm{c}^{\circ3}} & 4\\[0.3em] 4 & [\bullet[\bullet]] & f_{yy}(f,f_yf) & {\color{teal}(\bm{b}\circ\bm{c})\cdot\bm{d}} & 8\\[0.3em] 4 & [[\bullet^2]] & f_yf_{yy}(f,f) & {\color{teal}\bm{e}\cdot\bm{c}^{\circ2}} & 12\\[0.3em] 4 & [[[\bullet]]] & f_yf_yf_yf & {\color{teal}\bm{e}\cdot\bm{d}} & 24 \end{array}
Two structural facts follow at once and are worth keeping in mind.
- The right-hand side 1/\gamma(\tau) depends only on the tree, never on the method: the conditions are the same for explicit and implicit schemes, and only the solvability of the system differs.
- A tree of height m — a chain of m nested brackets — produces a condition involving \bm{A}^{m-1}. For an explicit method \bm{A} is strictly lower triangular and nilpotent, so \bm{A}^{m-1}=\bm 0 once m-1\ge s: the corresponding condition becomes 0=1/\gamma(\tau) and is unsatisfiable. This is the elementary reason why an explicit s-stage method cannot exceed order s.
How many conditions?
The number of rooted trees with p nodes grows quickly, and so does the cumulative number of conditions that must hold for order p:
| p | new conditions | total through order p |
|---|---|---|
| 1 | 1 | 1 |
| 2 | 1 | 2 |
| 3 | 2 | 4 |
| 4 | 4 | 8 |
| 5 | 9 | 17 |
| 6 | 20 | 37 |
| 7 | 48 | 85 |
| 8 | 115 | 200 |
| 9 | 286 | 486 |
| 10 | 719 | 1205 |
Constructing a method of order eight by solving two hundred polynomial equations by hand is out of the question; this is why high-order methods are built using the structural shortcuts of Section 4.2.4.1.1 rather than by brute force.
Attainable order: the Butcher barriers
For explicit methods the number of stages limits the order:
| order p | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
|---|---|---|---|---|---|---|---|---|
| minimum stages s | 1 | 2 | 3 | 4 | 6 | 7 | 9 | 11 |
Up to order four an extra stage buys an extra order; from order five on it does not, and the gap widens. This is the practical reason why RK4 occupies the position it does, and why high-order explicit codes (Dormand–Prince 8(7), for instance) use a dozen stages or more.
For implicit methods the bound is p\le2s, attained by the Gauss–Legendre family; Radau methods reach 2s-1 and Lobatto methods 2s-2, trading one or two orders for the stiff behaviour that the position of the abscissae gives them.
Simplifying assumptions
High-order implicit methods are not obtained by solving the order conditions directly, but through three families of structural identities due to Butcher. For integers p,\eta,\xi:
\begin{aligned} B(p):\quad & \sum_{i=1}^{s}b_ic_i^{\,k-1}=\frac1k, && k=1,\ldots,p;\\[0.4em] C(\eta):\quad & \sum_{j=1}^{s}a_{ij}c_j^{\,k-1}=\frac{c_i^{\,k}}{k}, && i=1,\ldots,s,\quad k=1,\ldots,\eta;\\[0.4em] D(\xi):\quad & \sum_{i=1}^{s}b_ic_i^{\,k-1}a_{ij}=\frac{b_j}{k}\bigl(1-c_j^{\,k}\bigr), && j=1,\ldots,s,\quad k=1,\ldots,\xi. \end{aligned} \tag{4.11}
Their meaning is concrete: B(p) says that the weights b_i integrate polynomials of degree <p exactly on [0,1] — it is a quadrature condition; C(\eta) says that each stage integrates polynomials of degree <\eta exactly on [0,c_i]; D(\xi) is the adjoint statement.
Theorem 4.2 (Butcher’s theorem) If a Runge–Kutta method satisfies B(p), C(\eta) and D(\xi) with
p\le\eta+\xi+1 \qquad\text{and}\qquad p\le2\eta+2,
then it has order p.
This reduces the construction of a high-order implicit method to a quadrature problem. Choosing the c_i as the nodes of a Gauss, Radau or Lobatto quadrature rule on [0,1] fixes B(p); the linear systems C(\eta) and D(\xi) then determine \bm{A} and \bm{b}. The classical families are exactly the three choices of nodes:
| family | abscissae | B | C | D | order |
|---|---|---|---|---|---|
| Gauss–Legendre | interior | 2s | s | D(s) | 2s |
| Radau IA | c_1=0 | 2s-1 | s-1 | D(s) | 2s-1 |
| Radau IIA | c_s=1 | 2s-1 | s | D(s-1) | 2s-1 |
| Lobatto IIIA | c_1=0,\ c_s=1 | 2s-2 | s | D(s-2) | 2s-2 |
| Lobatto IIIB | c_1=0,\ c_s=1 | 2s-2 | s-2 | D(s) | 2s-2 |
| Lobatto IIIC | c_1=0,\ c_s=1 | 2s-2 | s-1 | D(s-1) | 2s-2 |
The three-stage members of each family, whose tableaux appear below, illustrate the pattern: Gauss–Legendre reaches order six with three stages, Radau IIA order five, the three Lobatto variants order four. The Lobatto variants differ from one another not in order but in stability: IIIA and IIIB share the same stability function, with R(\infty)=1, so they are A-stable but not L-stable, whereas IIIC has R(\infty)=0 and damps the stiff modes completely. Stability is the subject of the next chapter.
- To check a given tableau, use the vector form Equation 4.9 and work level by level; Section 7.1 does this for a dozen methods and Section 7.3 automates it.
- To construct an explicit method of low order, write the conditions as a polynomial system in the free coefficients and solve it; the two-stage case is done in full in 1.
- To construct an implicit method of high order, never touch the order conditions: choose the quadrature nodes and use Theorem 4.2.
Embedded Runge–Kutta Pairs and Adaptive Step Size
In practical ODE solvers the step size is rarely kept fixed. A standard strategy is to use an embedded Runge–Kutta pair: the same stages \bm{K}_i are combined with two different sets of weights, \bm{b} and \widehat{\bm{b}}, producing two approximations of different orders,
\bm{z}_{k+1} =\bm{z}_k+\sum_{i=1}^s b_i\bm{K}_i, \qquad \widehat{\bm{z}}_{k+1} =\bm{z}_k+\sum_{i=1}^s \widehat b_i\bm{K}_i.
Their difference
\boxed{ \bm{e}_{k+1}^{\rm est} = \bm{z}_{k+1}-\widehat{\bm{z}}_{k+1} = \sum_{i=1}^s(b_i-\widehat b_i)\bm{K}_i }
is an inexpensive estimate of the local error because no additional evaluations of \bm{f} are required.
For a vector problem, the error is normally scaled componentwise using absolute and relative tolerances. One common norm is
\mathrm{err} = \left[ \frac1d\sum_{j=1}^d \left( \frac{e_j^{\rm est}} {\mathrm{atol}_j+\mathrm{rtol}_j\max(|z_{k,j}|,|z_{k+1,j}|)} \right)^2 \right]^{1/2}.
A step is accepted when \mathrm{err}\le1. If the estimator behaves like C\Delta t^q, a basic step-size controller is
\boxed{ \Delta t_{\rm new} = \eta \Delta t\,\mathrm{err}^{-1/q} }
with a safety factor 0<\eta<1 and practical upper/lower bounds on the change factor. For a common embedded pair of orders p and p-1, the difference is typically O(\Delta t^p), so q=p.
Adaptive step-size control is one of the main reasons modern explicit RK methods such as Bogacki–Shampine and Dormand–Prince are substantially more efficient than fixed-step RK4 in general-purpose software.
Example: Bogacki–Shampine 3(2) pair
A concrete embedded pair is the Bogacki–Shampine method, used in MATLAB’s ode23. Its four stages are determined by
\begin{array}{c|cccc} 0&0&0&0&0\\ \frac12&\frac12&0&0&0\\ \frac34&0&\frac34&0&0\\ 1&\frac29&\frac13&\frac49&0\\ \hline &\frac29&\frac13&\frac49&0\\ &\frac7{24}&\frac14&\frac13&\frac18 \end{array}
The first weight vector gives an order-three approximation,
\bm{z}_{n+1}^{[3]} = \bm{z}_n +\Delta t\sum_{i=1}^4 b_i\bm{k}_i,
whereas the second gives an order-two approximation,
\bm{z}_{n+1}^{[2]} = \bm{z}_n +\Delta t\sum_{i=1}^4 \widehat b_i\bm{k}_i.
The difference
\bm{e}_{n+1}^{\rm est} = \bm{z}_{n+1}^{[3]}-\bm{z}_{n+1}^{[2]}
is O(\Delta t^3) and is used to control the local error. An important implementation feature is FSAL (first same as last): after an accepted step, the last stage
\bm{k}_4=\bm{f}(\bm{z}_{n+1}^{[3]},t_{n+1})
can be reused as the first stage of the next step. Thus, after the first accepted step, only three new right-hand-side evaluations are normally required.
A common controller has the form
\Delta t_{\rm new} = \eta \Delta t\,\mathrm{err}^{-1/3},
because the embedded difference is of order \Delta t^3. In robust implementations this factor is clipped to a prescribed interval to prevent abrupt changes of step size.
The stability function R(q) of a Runge–Kutta method, and the absolute-stability region it defines, are derived from the Butcher tableau in Section 6.4.5.
4.3 Examples of Runge–Kutta Methods
The tableaux below are collected here for reference. Each one can be checked against the conditions of Equation 4.9; the verification is carried out explicitly for a dozen of them in Section 7.1, and automated in Section 7.3.
Explicit methods
Explicit Euler.
The Explicit Euler Method Represented Using a Butcher Tableau:
\def\arraystretch{1.5} \begin{array}{c|c} 0 & 0 \\ \hline & 1 \end{array} \tag{4.12}
Midpoint method.
The Midpoint Method, also known as the Collatz Method, is a variant of the explicit Euler Method. This approach is classified as a second-order method and consists of two stages, allowing for improved accuracy in numerical integration.
\def\arraystretch{1.5} \begin{array}{c|cc} 0 & 0 & 0 \\ 1/2 & 1/2 & 0 \\ \hline & 0 & 1 \end{array} \tag{4.13}
Heun’s method.
Heun’s Method, often referred to as the explicit trapezoid rule, enhances the standard explicit Euler Method. This technique is a second-order method that employs two stages to achieve greater accuracy in numerical integration.
\def\arraystretch{1.5} \begin{array}{c|cc} 0 & 0 & 0 \\ 1 & 1 & 0 \\ \hline & 1/2 & 1/2 \end{array} \tag{4.14}
Ralston’s method.
Ralston’s method is a second-order technique with two stages, specifically designed to achieve a minimal local error bound.
\def\arraystretch{1.5} \begin{array}{c|cc} 0 & 0 & 0 \\ 2/3 & 2/3 & 0 \\ \hline & 1/4 & 3/4 \end{array} \tag{4.15}
Generic second-order method.
\def\arraystretch{1.5} \begin{array}{c|cc} 0 & 0 & 0 \\ \alpha & \alpha & 0 \\ \hline\\[-1em] & 1-\frac{1}{2\alpha} & \frac{1}{2\alpha} \end{array} \qquad \alpha\neq0. \tag{4.16}
If one additionally wants the internal stage time to lie in the current step, one restricts 0<\alpha\le1. The endpoint \alpha=1 gives Heun’s method and must therefore not be excluded.
Runge–Kutta’s third-order method (RK3).
\def\arraystretch{1.5} \begin{array} {c|ccc} 0 & 0 & 0 & 0 \\ 1/2 & 1/2 & 0 & 0 \\ 1 & -1 & 2 & 0 \\ \hline & 1/6 & 2/3 & 1/6 \end{array} \tag{4.17}
Heun’s third-order method.
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1/3 & 1/3 & 0 & 0 \\ 2/3 & 0 & 2/3 & 0 \\ \hline & 1/4 & 0 & 3/4 \end{array} \tag{4.18}
Ralston’s third-order method.
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1/2 & 1/2 & 0 & 0 \\ 3/4 & 0 & 3/4 & 0 \\ \hline & 2/9 & 1/3 & 4/9 \end{array} \tag{4.19}
Strong stability preserving Runge–Kutta method (SSPRK3).
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1/2 & 1/4 & 1/4 & 0 \\ \hline & 1/6 & 1/6 & 2/3 \end{array} \tag{4.20}
Generic third-order method.
A one-parameter family of explicit three-stage, third-order methods with c_2=\alpha and c_3=1 is
\def\arraystretch{1.7} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \alpha & \alpha & 0 & 0 \\ 1 & \dfrac{3\alpha^2-3\alpha+1}{\alpha(3\alpha-2)} & \dfrac{\alpha-1}{\alpha(3\alpha-2)} & 0 \\[0.6em] \hline & \dfrac{3\alpha-1}{6\alpha} & \dfrac{1}{6\alpha(1-\alpha)} & \dfrac{2-3\alpha}{6(1-\alpha)} \end{array} \tag{4.21}
with
\alpha\neq0,\qquad \alpha\neq\frac23,\qquad \alpha\neq1.
Substitution into the third-order conditions gives
\bm{b}^T\bm{1}=1,\qquad \bm{b}^T\bm{c}=\frac12,\qquad \bm{b}^T\bm{c}^{\circ2}=\frac13,\qquad \bm{b}^T\bm{A}\bm{c}=\frac16.
These identities provide a direct check that the family has order three.
Classical Runge–Kutta method (RK4).
The classical fourth-order Runge-Kutta:
\def\arraystretch{1.5} \begin{array} {c|cccc} 0 & 0 & 0 & 0 & 0 \\ 1/2 & 1/2 & 0 & 0 & 0 \\ 1/2 & 0 & 1/2 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ \hline & 1/6 & 2/6 & 2/6 & 1/6 \end{array} \tag{4.22}
Ralston’s fourth-order method.
It is a four-stage, fourth-order Runge–Kutta method whose coefficients were chosen to reduce the leading truncation-error terms.
\def\arraystretch{1.5} \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0 \\ 2/5 & 2/5 & 0 & 0 & 0 \\ c_3 & a_{3,1} & a_{3,2} & 0 & 0 \\ 1 & a_{4,1} & a_{4,2} & a_{4,3} & 0 \\ \hline\\ & b_1 & b_2 & b_3 & b_4 \end{array} \tag{4.23}
where
\begin{aligned} b_1 & = \frac{263}{1812}+\frac{2}{151}\sqrt{5}\\ b_2 & = \frac{125}{3828}-\frac{250}{957}\sqrt{5} \\ b_3 & = \frac{3426304}{5924787}+\frac{553984}{1974929}\sqrt{5} \\ b_4 & = \frac{10}{41}-\frac{4}{123}\sqrt{5} \\ c_3 & = \frac{7}{8}-\frac{3}{16}\sqrt{5} \\ a_{3,1} & = \frac{357}{256}\sqrt{5}-\frac{2889}{1024} \\ a_{4,1} & = \frac{1047}{3020}\sqrt{5}-\frac{673}{1208} \\ a_{3,2} &= \frac{3785}{1024}-\frac{405}{256}\sqrt{5} \\ a_{4,2} &= -\frac{975}{2552}-\frac{1523}{1276}\sqrt{5} \\ a_{4,3} &= \frac{93408}{48169}+\frac{203968}{240845}\sqrt{5} \end{aligned}
Fourth-order 3/8-rule method.
The 3/8 rule is a four-stage, fourth-order explicit Runge–Kutta method.
\def\arraystretch{1.5} \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0 \\ 1/3 & 1/3 & 0 & 0 & 0 \\ 2/3 & -1/3 & 1 & 0 & 0 \\ 1 & 1 & -1 & 1 & 0 \\ \hline & 1/8 & 3/8 & 3/8 & 1/8 \end{array}
Implicit methods
The one- and two-stage schemes below are the implicit counterparts of the explicit methods above; the Gauss–Legendre, Lobatto and Radau families that follow are the high-order methods built from quadrature nodes through the simplifying assumptions of Section 4.2.4.1.1.
Basic implicit methods
Implicit Euler method.
The Implicit Euler Method Represented Using a Butcher Tableau:
\def\arraystretch{1.5} \begin{array}{c|c} 1 & 1 \\ \hline & 1 \end{array} \tag{4.24}
Implicit midpoint method.
The implicit midpoint rule is a one-stage, second-order implicit Runge–Kutta method. It should not be confused with the explicit midpoint (Collatz) method, which has two stages.
\def\arraystretch{1.5} \begin{array}{c|c} 1/2 & 1/2 \\ \hline & 1 \end{array} \tag{4.25}
Crank–Nicolson method.
The Crank-Nicolson method corresponds to the implicit trapezoidal rule.
\def\arraystretch{1.5} \begin{array}{c|cc} 0 & 0 & 0 \\ 1 & 1/2 & 1/2 \\ \hline & 1/2 & 1/2 \end{array} \tag{4.26}
Methods from quadrature nodes
Gauss–Legendre methods.
These methods are based on the points of Gauss-Legendre quadrature.
Two-stage Gauss–Legendre method (order 4).
\def\arraystretch{3} \begin{array}{c|cc} \dfrac{1}{2}-\dfrac{\sqrt3}{6} & \dfrac{1}{4} & \dfrac{1}{4}-\dfrac{\sqrt3}{6} \\ \dfrac{1}{2}+\dfrac{\sqrt3}{6} & \dfrac{1}{4}+\dfrac{\sqrt3}{6} & \dfrac{1}{4} \\ \hline & \dfrac{1}{2} & \dfrac{1}{2} \end{array} \tag{4.27}
Gauss–Legendre sixth-order method.
\def\arraystretch{3} \begin{array}{c|ccc} \dfrac{1}{2} - \dfrac{\sqrt{15}}{10} & \dfrac{5}{36} & \dfrac{2}{9}- \dfrac{\sqrt{15}}{15} & \dfrac{5}{36} - \dfrac{\sqrt{15}}{30} \\ \dfrac{1}{2} & \dfrac{5}{36} + \dfrac{\sqrt{15}}{24} & \dfrac{2}{9} & \dfrac{5}{36} - \dfrac{\sqrt{15}}{24} \\ \dfrac{1}{2} + \dfrac{\sqrt{15}}{10} & \dfrac{5}{36} + \dfrac{\sqrt{15}}{30} & \dfrac{2}{9} + \dfrac{\sqrt{15}}{15} & \dfrac{5}{36} \\[0.5em] \hline & \dfrac{5}{18} & \dfrac{4}{9} & \dfrac{5}{18} \end{array} \tag{4.28}
Lobatto methods.
Lobatto IIIA fourth-order method.
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1/2 & 5/24& 1/3 & -1/24 \\ 1 & 1/6 & 2/3 & 1/6 \\ \hline & 1/6 & 2/3 & 1/6 \end{array} \tag{4.29}
Lobatto IIIB fourth-order method.
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 1/6 & -1/6 & 0 \\ 1/2 & 1/6 & 1/3 & 0 \\ 1 & 1/6 & 5/6 & 0 \\ \hline & 1/6 & 2/3 & 1/6 \end{array} \tag{4.30}
Lobatto IIIC fourth-order method.
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 1/6 & -1/3 & 1/6 \\ 1/2 & 1/6 & 5/12 & -1/12 \\ 1 & 1/6 & 2/3 & 1/6 \\ \hline & 1/6 & 2/3 & 1/6 \end{array} \tag{4.31}
Lobatto IIIC^\star fourth-order method.
\def\arraystretch{1.5} \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ 1/2 & 1/4 & 1/4 & 0 \\ 1 & 0 & 1 & 0 \\ \hline & 1/6 & 2/3 & 1/6 \end{array} \tag{4.32}
Radau methods.
Radau IA fifth-order method.
\def\arraystretch{3} \begin{array}{c|ccc} 0 & \dfrac{1}{9} & \dfrac{-1 - \sqrt{6}}{18} & \dfrac{-1 + \sqrt{6}}{18} \\ \dfrac{3}{5} - \dfrac{\sqrt{6}}{10} & \dfrac{1}{9} & \dfrac{11}{45} + \dfrac{7\sqrt{6}}{360} & \dfrac{11}{45} - \dfrac{43\sqrt{6}}{360}\\ \dfrac{3}{5} + \dfrac{\sqrt{6}}{10} & \dfrac{1}{9} & \dfrac{11}{45} + \dfrac{43\sqrt{6}}{360} & \dfrac{11}{45} - \dfrac{7\sqrt{6}}{360} \\ \hline & \dfrac{1}{9} & \dfrac{4}{9} + \dfrac{\sqrt{6}}{36} & \dfrac{4}{9} - \dfrac{\sqrt{6}}{36} \end{array} \tag{4.33}
Radau IIA fifth-order method.
\def\arraystretch{3} \begin{array}{c|ccc} \dfrac{2}{5} - \dfrac{\sqrt{6}}{10} & \dfrac{11}{45} - \frac{7\sqrt{6}}{360} & \dfrac{37}{225} - \dfrac{169\sqrt{6}}{1800} & -\dfrac{2}{225} + \dfrac{\sqrt{6}}{75} \\ \dfrac{2}{5} + \dfrac{\sqrt{6}}{10} & \dfrac{37}{225} + \dfrac{169\sqrt{6}}{1800} & \dfrac{11}{45} + \dfrac{7\sqrt{6}}{360} & -\dfrac{2}{225} - \dfrac{\sqrt{6}}{75}\\ 1 & \dfrac{4}{9} - \dfrac{\sqrt{6}}{36} & \dfrac{4}{9} + \dfrac{\sqrt{6}}{36} & \dfrac{1}{9} \\ \hline & \dfrac{4}{9} - \dfrac{\sqrt{6}}{36} & \dfrac{4}{9} + \dfrac{\sqrt{6}}{36} & \dfrac{1}{9} \end{array} \tag{4.34}
Summary
- A Runge–Kutta method replaces the derivatives of \bm{f} required by a Taylor method with several evaluations of \bm{f} at intermediate points. The Butcher tableau (\bm{A},\bm{b},\bm{c}) contains all the information about the scheme, and the stage-consistency relation \bm{A}\bm{1}=\bm{c} is presupposed everywhere.
- In general a step is a nonlinear system of dimension sd. The shape of \bm{A} decides the cost: full means one sd\times sd system and up to order 2s; lower triangular (DIRK, SDIRK) means s systems of size d, with a single Newton matrix per step in the SDIRK case; strictly lower triangular means no system at all, s evaluations of \bm{f}, but order at most s and a bounded stability region.
- Fixed-point iteration on the stages is unusable on stiff problems because its contraction condition \Delta t L\max_i\sum_j|a_{ij}|<1 reinstates exactly the step restriction that implicitness was meant to remove; Newton’s method, in simplified or DIRK form, is what makes implicitness worthwhile.
- The order conditions follow from matching the Taylor expansion of the exact solution with that of the increment function. The two-stage case
- shows the mechanism in full; the general case
- gives one condition per rooted tree, \Phi(\tau)=1/\gamma(\tau), eight of them through order four and seventeen through order five.
- Explicit s-stage methods have order at most s, and from order five on they need strictly more stages than the order; implicit s-stage methods reach at most 2s. High-order implicit methods are not built by solving the order conditions but by choosing quadrature nodes and applying the simplifying assumptions B(p), C(\eta), D(\xi).
- Embedded pairs provide a local error estimate at no extra cost and are the basis of adaptive step-size control; the FSAL property makes the extra stage almost free.