Butcher Tableaux and the Order Conditions of Runge-Kutta Methods

A Runge-Kutta method replaces the single tangent evaluation of Euler with several evaluations inside one step, combined so that the extra work buys extra accuracy without extra past states. The combination is written as a table of coefficients, and the accuracy of the resulting method is decided by a set of algebraic identities among those coefficients. Understanding those identities explains why fourth order costs four stages, why no fifth-order method exists with four, and where the order barriers of Runge-Kutta theory come from.

The tableau as a data structure

An $s$-stage explicit or implicit scheme computes stage derivatives $k_i$ and then a weighted sum:

$$ k_i = f!\left(t_n + c_i h,; y_n + h\sum_{j=1}^{s} a_{ij} k_j\right), \qquad y_{n+1} = y_n + h\sum_{i=1}^{s} b_i k_i $$

The constants $a_{ij}$, $b_i$, $c_i$ are usually arranged as a Butcher tableau: a lower-triangular $a$ matrix with the row sums on the side as $c$, and the $b$ row underneath. For explicit methods $a_{ij} = 0$ for $j \ge i$, so each stage depends only on previously computed ones.

The classic fourth-order method, often attributed to Kutta, is worth writing out because its coefficients are the canonical example:

0   |
1/2 | 1/2
1/2 | 0    1/2
1   | 0    0    1
----+------------------
    | 1/6  1/3  1/3  1/6

The c column gives the time offset of each stage; the row sums of a equal the corresponding c for any consistent method, which is automatic for the explicit case above. The b row weights the stages in the final combination.

Two implementation consequences follow from the tableau’s shape. Explicit schemes need only one y work array plus a buffer for the current stage, because stage $i$ reads states formed from accumulated k values; a straightforward implementation keeps s history vectors of size $n$, giving $8sn$ bytes in double precision. Dense output between steps requires either a separate interpolating polynomial or a method with a continuous extension, since the stages only describe the interior of the step.

What “order $p$” means in terms of the coefficients

Consistency is the condition that the method reproduces a constant right-hand side. If $f$ does not depend on $y$, the exact solution is linear and every reasonable method gets it right; the algebraic statement is $\sum_i b_i = 1$. First order additionally requires the method to integrate a linear right-hand side exactly, which amounts to $\sum_i b_i c_i = 1/2$.

Each order imposes further polynomial identities. The systematic way to enumerate them is through rooted trees, which encode the partial derivatives of the exact solution. For $y’ = f$ the second derivative is $f’ f$, the third is $f” (f,f) + f’ f’ f$, and so on; each distinct way of composing derivatives corresponds to a rooted tree, and each tree contributes one order condition. Counting trees by order: order 1 has 1 tree, order 2 has 1, order 3 has 2, order 4 has 4, order 5 has 9, order 6 has 20, order 7 has 48. The elementary weights $\Phi_i$ assigned to each tree are polynomials in the tableau entries; the condition is that $\sum_j b_j \Phi_j(t) = 1/\gamma(t)$ where $\gamma(t)$ counts the symmetries of the tree.

For a fourth-order method with four stages, the conditions are:

sum(b_i)      = 1                order 1
sum(b_i c_i)  = 1/2              order 2
sum(b_i c_i^2) = 1/3             order 3
sum(b_i a_ij c_j) = 1/6          order 3
sum(b_i c_i^3) = 1/4             order 4
sum(b_i c_i a_ij c_j) = 1/8      order 4
sum(b_i a_ij c_j^2) = 1/12       order 4
sum(b_i a_ij a_jk c_k) = 1/24    order 4

That is eight equations in the tableau entries, counting the roots of the trees at each order. The tableau above satisfies all of them. Adding order raises the count combinatorially, and the number of stages has to grow at least as fast as the number of conditions, which is the crude reason high-order methods are expensive.

The barriers: no free stages

Two structural results constrain what tableaux can exist.

The first is about explicit methods. With $s$ stages, the maximum order is $s$ for $s \le 4$. For $s = 5$ the maximum order is still 4; the tableau has eleven or more free entries and the fifth-order conditions still cannot all be met. The reason is that the stage values themselves represent approximations with their own truncation errors, and the order conditions couple them in a way that cannot be satisfied with the available degrees of freedom. Reaching fifth order with an explicit method needs six stages. Six stages give at most fifth order; seventh order needs nine stages, eighth order needs eleven.

The second is about composition and the step. Every Runge-Kutta method, applied to a linear test problem, has a stability function $R(z) = 1 + z b^T(I - zA)^{-1}\mathbf{1}$. For an $s$-stage method, $R$ is a rational function whose numerator has degree at most $s$; to approximate $e^z$ to order $p$ the first $p+1$ Taylor coefficients must match, so $p \le s$ for implicit methods and $p \le s$ for explicit ones, with the extra stage shortfall for the explicit case as above.

There is a further price for high order in explicit methods: the stability interval on the negative real axis, where the scheme damps a decaying mode, shrinks as order rises. A lower-order method has a wider stability region per unit of work, which is why a low-order scheme with many cheap steps can beat a high-order scheme with few expensive ones for problems where the accuracy requirement is loose.

Where this connects to stiffness and to step control

The tableau framework makes explicit the two things an integrator must supply: a method with the required order and a mechanism for estimating the error. Embedding solves both at once. A method is embedded when a second set of weights $\hat b$ of lower order is available for the same stages, so the difference between the two combinations gives a local error estimate at the cost of the stage values already computed. The difference in weights multiplies the k values accumulated in the same pass, so the estimate is nearly free. Families of embedded pairs, indexed by the orders $(p, \hat p)$ they provide, are the basis for adaptive stepping described in step-size control for embedded pairs.

Order alone does not make a method suitable for a stiff problem. High-order explicit methods have bounded stability regions and are hopeless when the Jacobian’s eigenvalues are large and negative. Implicit tableaux with $a_{ij} \ne 0$ for $j \ge i$ are the alternative, but their stability functions are rational approximations to $e^z$ with a pole structure that determines whether they damp the infinite modes. The analysis of that damping is the content of the Dahlquist barrier and related results, discussed in A-stability and the Dahlquist barrier.

Reading a tableau before trusting it

The practical diagnostic when choosing a method is to check three properties of the tableau directly, without running anything. Consistency sums: $\sum b_i = 1$ and $\sum b_i c_i = 1/2$ catch transcription errors and mistyped coefficients, which are a real source of wrong results when a tableau is copied by hand. Stage consistency: for a consistent method $\sum_{j} a_{ij} = c_i$ for every stage $i$, which is a bookkeeping check on the interior rows. Zero-pattern: for an explicit method, $a_{ij}$ must vanish for $j \ge i$. A single nonzero on or above the diagonal turns the method implicit, and the implementation must then solve a system that it was not written for.

A fourth property matters for a class of methods used to preserve invariants: whether the method is symplectic or, more loosely, whether it has a nontrivial conservation law. For a separable Hamiltonian system, no explicit Runge-Kutta method of the form above is symplectic; the order conditions for the invariant conflict with the order conditions for accuracy. That conflict is why Hamiltonian integration uses different machinery, examined in symplectic integrators and energy drift.

For most work the conclusion is deflationary: order 4 or order 5 with an embedded error estimate covers the overwhelming majority of non-stiff problems, the coefficients are standard, and the real decisions are the step controller, the tolerance, and the handling of the right-hand side’s cost. Higher-order methods pay off when the right-hand side is expensive relative to the arithmetic of the method itself, which is the case when each evaluation costs a mesh solve or a large sparse matrix-vector product.