7 Checking the Order of a RK Method

Author

Affiliation

Enrico Bertolazzi

University of Trento, Department of Industrial Engineering

The order conditions collected in Section 4.2.4 are usually quoted and rarely used. This chapter does the opposite: it applies them, one tableau at a time, to the methods introduced in the previous chapters. The point of the exercise is twofold — to make the conditions concrete, and to show that determining the order of a method is a finite, purely algebraic computation that needs no Taylor expansion and no numerical experiment.

7.1 The Conditions in Vector Form

The conditions themselves, their derivation and the rooted-tree rule behind them are in Section 4.2.4. This section only recalls the notation and then puts it to work.

Let the method have s stages, with Butcher tableau

\begin{array}{c|c} \bm{c} & \bm{A}\\ \hline & \bm{b}^T \end{array}

and let \bm{1}=(1,\ldots,1)^T\in\mathbb R^s. With \circ the componentwise (Hadamard) product, write \bm{c}^{\circ2}=\bm{c}\circ\bm{c}, \bm{c}^{\circ3}=\bm{c}\circ\bm{c}\circ\bm{c}, and introduce the two auxiliary vectors of Equation 4.8,

\bm{d}=\bm{A}\bm{c}, \qquad \bm{e}=\bm{A}^T\bm{b}.

The eight conditions through order four are then those of Equation 4.9, preceded by the stage-consistency relation \bm{A}\bm{1}=\bm{c} of Equation 4.7:

\begin{array}{r|l} \text{order} & \text{new conditions}\\ \hline 0 & \bm{A}\bm{1}=\bm{c}\quad\text{(stage consistency)}\\[0.3em] 1 & \bm{b}\cdot\bm{1}=1\\[0.3em] 2 & \bm{b}\cdot\bm{c}=\dfrac12\\[0.6em] 3 & \bm{b}\cdot\bm{c}^{\circ2}=\dfrac13, \qquad \bm{b}\cdot\bm{d}=\dfrac16\\[0.6em] 4 & \bm{b}\cdot\bm{c}^{\circ3}=\dfrac14, \qquad (\bm{b}\circ\bm{c})\cdot\bm{d}=\dfrac18,\\[0.4em] & \bm{e}\cdot\bm{c}^{\circ2}=\dfrac1{12}, \qquad \bm{e}\cdot\bm{d}=\dfrac1{24} \end{array}

The conditions are cumulative: a method has order p when all the conditions listed for orders 1,\ldots,p hold and at least one condition of order p+1 fails. The row labelled 0 is not an order condition by itself, but every formula below presupposes it and it is worth verifying first because a tableau that violates it has been transcribed incorrectly.

How to organise the computation by hand

Build one table per method, adding rows as needed:

write \bm{b} and check that its entries sum to 1;
write \bm{c} (the row sums of \bm{A}), then \bm{b}\circ\bm{c} and sum;
form \bm{c}^{\circ2}, \bm{c}^{\circ3} by squaring and cubing entries;
form \bm{d}=\bm{A}\bm{c} and \bm{e}=\bm{A}^T\bm{b} once, and reuse them.

Stop as soon as a condition fails: the order is then one less than the level of the failing condition. Everything is rational arithmetic, so the checks are exact.

7.2 Worked Verifications

The nine tableaux below are checked one by one with the conditions in vector form. Each verification is self-contained: read the one you need.

Explicit Euler

\def\arraystretch{1.4} \begin{array}{c|c} 0 & 0 \\ \hline & 1 \end{array} \qquad \bm{A}=(0),\quad \bm{b}=(1),\quad \bm{c}=(0).

Stage consistency: \bm{A}\bm{1}=(0)=\bm{c}. Then

\begin{array}{l|c|l} \text{condition} & \text{value} & \\ \hline \bm{b}\cdot\bm{1}=1 & 1 & \text{OK}\\ \bm{b}\cdot\bm{c}=\frac12 & 0 & \text{fails} \end{array}

Order 1. The failure is immediate: a single stage evaluated at the left endpoint cannot see the curvature of the solution.

Explicit Midpoint

\def\arraystretch{1.4} \begin{array}{c|cc} 0 & 0 & 0 \\ 1/2 & 1/2 & 0 \\ \hline & 0 & 1 \end{array} \qquad \bm{c}=\begin{pmatrix}0\\ 1/2\end{pmatrix},\quad \bm{b}=\begin{pmatrix}0\\ 1\end{pmatrix},\quad \bm{d}=\bm{A}\bm{c}=\begin{pmatrix}0\\ 0\end{pmatrix}.

\begin{array}{c||cc||cc} \hline \bm{b} & 0 & 1 & 1 & \text{OK}\\ \hline\hline \bm{c} & 0 & 1/2 & & \\ \hline \bm{b}\circ\bm{c} & 0 & 1/2 & 1/2 & \text{OK}\\ \hline\hline \bm{c}^{\circ2} & 0 & 1/4 & & \\ \hline \bm{b}\circ\bm{c}^{\circ2} & 0 & 1/4 & 1/4 & \text{fails }(\neq 1/3)\\ \hline \end{array}

Order 2. Once a third-order condition fails there is no need to test \bm{b}\cdot\bm{d}=\frac16: the order is already determined.

Ralston’s Second-Order Method

\def\arraystretch{1.4} \begin{array}{c|cc} 0 & 0 & 0 \\ 2/3 & 2/3 & 0 \\ \hline & 1/4 & 3/4 \end{array} \qquad \bm{c}=\begin{pmatrix}0\\ 2/3\end{pmatrix},\quad \bm{d}=\bm{A}\bm{c}=\begin{pmatrix}0\\ 0\end{pmatrix}.

\begin{array}{c||cc||cc} \hline \bm{b} & 1/4 & 3/4 & 1 & \text{OK}\\ \hline\hline \bm{c} & 0 & 2/3 & & \\ \hline \bm{b}\circ\bm{c} & 0 & 1/2 & 1/2 & \text{OK}\\ \hline\hline \bm{c}^{\circ2} & 0 & 4/9 & & \\ \hline \bm{b}\circ\bm{c}^{\circ2} & 0 & 1/3 & 1/3 & \text{OK}\\ \hline \bm{d} & 0 & 0 & & \\ \hline \bm{b}\circ\bm{d} & 0 & 0 & 0 & \text{fails }(\neq 1/6)\\ \hline \end{array}

Order 2. This example is instructive: Ralston’s method satisfies \bm{b}\cdot\bm{c}^{\circ2}=\frac13, one of the two third-order conditions, and fails the other. Satisfying some conditions of a given level buys nothing; the choice c_2=2/3 minimises the leading error constant among second-order two-stage methods, but it does not raise the order.

Third-Order Runge–Kutta (RK3)

\def\arraystretch{1.4} \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} \qquad \bm{c}=\begin{pmatrix}0\\ 1/2\\ 1\end{pmatrix},\quad \bm{d}=\bm{A}\bm{c}=\begin{pmatrix}0\\ 0\\ 1\end{pmatrix}.

\begin{array}{c||ccc||cc} \hline \bm{b} & 1/6 & 2/3 & 1/6 & 1 & \text{OK}\\ \hline\hline \bm{c} & 0 & 1/2 & 1 & & \\ \hline \bm{b}\circ\bm{c} & 0 & 1/3 & 1/6 & 1/2 & \text{OK}\\ \hline\hline \bm{c}^{\circ2} & 0 & 1/4 & 1 & & \\ \hline \bm{b}\circ\bm{c}^{\circ2} & 0 & 1/6 & 1/6 & 1/3 & \text{OK}\\ \hline \bm{d} & 0 & 0 & 1 & & \\ \hline \bm{b}\circ\bm{d} & 0 & 0 & 1/6 & 1/6 & \text{OK}\\ \hline\hline \bm{c}^{\circ3} & 0 & 1/8 & 1 & & \\ \hline \bm{b}\circ\bm{c}^{\circ3} & 0 & 1/12& 1/6 & 1/4 & \text{OK}\\ \hline (\bm{b}\circ\bm{c})\circ\bm{d} & 0 & 0 & 1/6 & 1/6 & \text{fails }(\neq 1/8)\\ \hline \end{array}

Order 3. Here the method passes the first fourth-order condition \bm{b}\cdot\bm{c}^{\circ3}=\frac14 and fails the second. Testing the conditions in the order listed in Equation 4.9 is therefore not an early-exit strategy: all conditions of a level must be checked before the level can be declared satisfied.

Fourth-Order Runge–Kutta: the 3/8 Rule

\def\arraystretch{1.4} \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}

Stage consistency gives \bm{c}=(0,\tfrac13,\tfrac23,1)^T, and

\bm{d}=\bm{A}\bm{c}=\begin{pmatrix}0\\ 0\\ 1/3\\ 1/3\end{pmatrix}, \qquad \bm{e}=\bm{A}^T\bm{b}=\begin{pmatrix}1/8\\ 1/4\\ 1/8\\ 0\end{pmatrix}.

\begin{array}{c||cccc||cc} \hline \bm{b} & 1/8 & 3/8 & 3/8 & 1/8 & 1 & \text{OK}\\ \hline\hline \bm{c} & 0 & 1/3 & 2/3 & 1 & & \\ \hline \bm{b}\circ\bm{c} & 0 & 1/8 & 1/4 & 1/8 & 1/2 & \text{OK}\\ \hline\hline \bm{c}^{\circ2} & 0 & 1/9 & 4/9 & 1 & & \\ \hline \bm{b}\circ\bm{c}^{\circ2} & 0 & 1/24 & 1/6 & 1/8 & 1/3 & \text{OK}\\ \hline \bm{d} & 0 & 0 & 1/3 & 1/3 & & \\ \hline \bm{b}\circ\bm{d} & 0 & 0 & 1/8 & 1/24 & 1/6 & \text{OK}\\ \hline\hline \bm{c}^{\circ3} & 0 & 1/27 & 8/27 & 1 & & \\ \hline \bm{b}\circ\bm{c}^{\circ3} & 0 & 1/72 & 1/9 & 1/8 & 1/4 & \text{OK}\\ \hline (\bm{b}\circ\bm{c})\circ\bm{d} & 0 & 0 & 1/12 & 1/24 & 1/8 & \text{OK}\\ \hline \bm{e} & 1/8 & 1/4 & 1/8 & 0 & & \\ \hline \bm{e}\circ\bm{c}^{\circ2} & 0 & 1/36 & 1/18 & 0 & 1/12 & \text{OK}\\ \hline \bm{e}\circ\bm{d} & 0 & 0 & 1/24 & 0 & 1/24 & \text{OK}\\ \hline \end{array}

Order at least 4. All eight conditions hold. Establishing that the order is exactly four would require the nine fifth-order conditions Equation 4.10; for a four-stage explicit method this is unnecessary, since the attainable order of an explicit s-stage method is at most s (Table 4.3).

Implicit Midpoint

\def\arraystretch{1.4} \begin{array}{c|c} 1/2 & 1/2 \\ \hline & 1 \end{array} \qquad \bm{A}=(1/2),\quad \bm{b}=(1),\quad \bm{c}=(1/2),\quad \bm{d}=\bm{A}\bm{c}=(1/4).

\begin{array}{l|c|l} \text{condition} & \text{value} & \\ \hline \bm{b}\cdot\bm{1}=1 & 1 & \text{OK}\\ \bm{b}\cdot\bm{c}=\frac12 & 1/2 & \text{OK}\\ \bm{b}\cdot\bm{c}^{\circ2}=\frac13 & 1/4 & \text{fails}\\ \bm{b}\cdot\bm{d}=\frac16 & 1/4 & \text{fails} \end{array}

Order 2. A single implicit stage already attains the maximum order 2s=2 of a one-stage method; no one-stage method of order three exists.

Crank–Nicolson (Trapezoidal Rule)

\def\arraystretch{1.4} \begin{array}{c|cc} 0 & 0 & 0 \\ 1 & 1/2 & 1/2 \\ \hline & 1/2 & 1/2 \end{array} \qquad \bm{c}=\begin{pmatrix}0\\ 1\end{pmatrix},\quad \bm{d}=\bm{A}\bm{c}=\begin{pmatrix}0\\ 1/2\end{pmatrix}.

\begin{array}{c||cc||cc} \hline \bm{b} & 1/2 & 1/2 & 1 & \text{OK}\\ \hline\hline \bm{c} & 0 & 1 & & \\ \hline \bm{b}\circ\bm{c} & 0 & 1/2 & 1/2 & \text{OK}\\ \hline\hline \bm{c}^{\circ2} & 0 & 1 & & \\ \hline \bm{b}\circ\bm{c}^{\circ2} & 0 & 1/2 & 1/2 & \text{fails }(\neq 1/3)\\ \hline \end{array}

Order 2. The trapezoidal rule is only a diagonally implicit two-stage method; its interest lies entirely in its stability (A-stable, but not L-stable), not in its order.

Gauss–Legendre of Order Four

\def\arraystretch{2.2} \begin{array}{c|cc} \dfrac12-\dfrac{\sqrt3}{6} & \dfrac14 & \dfrac14-\dfrac{\sqrt3}{6} \\ \dfrac12+\dfrac{\sqrt3}{6} & \dfrac14+\dfrac{\sqrt3}{6} & \dfrac14 \\ \hline & \dfrac12 & \dfrac12 \end{array}

Writing \gamma=\sqrt3/6 so that \bm{c}=(\tfrac12-\gamma,\ \tfrac12+\gamma)^T,

\bm{d}=\bm{A}\bm{c} =\begin{pmatrix}\tfrac16-\tfrac{\sqrt3}{12}\\[0.3em] \tfrac16+\tfrac{\sqrt3}{12}\end{pmatrix}, \qquad \bm{e}=\bm{A}^T\bm{b} =\begin{pmatrix}\tfrac14+\tfrac{\sqrt3}{12}\\[0.3em] \tfrac14-\tfrac{\sqrt3}{12}\end{pmatrix}.

The checks use repeatedly c_1+c_2=1 and c_1c_2=\tfrac14-\gamma^2=\tfrac16:

\begin{array}{l|l|l} \text{condition} & \text{value} & \\ \hline \bm{b}\cdot\bm{1}=1 & \tfrac12+\tfrac12=1 & \text{OK}\\[0.3em] \bm{b}\cdot\bm{c}=\tfrac12 & \tfrac12(c_1+c_2)=\tfrac12 & \text{OK}\\[0.3em] \bm{b}\cdot\bm{c}^{\circ2}=\tfrac13 & \tfrac12\bigl((c_1+c_2)^2-2c_1c_2\bigr)=\tfrac12\bigl(1-\tfrac13\bigr)=\tfrac13 & \text{OK}\\[0.3em] \bm{b}\cdot\bm{d}=\tfrac16 & \tfrac12\cdot\tfrac13=\tfrac16 & \text{OK}\\[0.3em] \bm{b}\cdot\bm{c}^{\circ3}=\tfrac14 & \tfrac14 & \text{OK}\\[0.3em] (\bm{b}\circ\bm{c})\cdot\bm{d}=\tfrac18 & \tfrac18 & \text{OK}\\[0.3em] \bm{e}\cdot\bm{c}^{\circ2}=\tfrac1{12} & \tfrac1{12} & \text{OK}\\[0.3em] \bm{e}\cdot\bm{d}=\tfrac1{24} & \tfrac1{24} & \text{OK} \end{array}

Order at least 4 — and in fact exactly 4=2s, the maximum attainable by an s-stage implicit method, which is the defining property of the Gauss–Legendre family. The surds cancel in every condition, as they must: the conditions are symmetric functions of c_1,c_2, whose elementary symmetric polynomials c_1+c_2=1 and c_1c_2=\tfrac16 are rational.

Lobatto IIIB of Order Four

\def\arraystretch{1.4} \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}

The row sums are 0,\tfrac12,1, so stage consistency holds, and

\bm{d}=\bm{A}\bm{c}=\begin{pmatrix}-1/12\\ 1/6\\ 5/12\end{pmatrix}, \qquad \bm{e}=\bm{A}^T\bm{b}=\begin{pmatrix}1/6\\ 1/3\\ 0\end{pmatrix}.

\begin{array}{c||ccc||cc} \hline \bm{b} & 1/6 & 2/3 & 1/6 & 1 & \text{OK}\\ \hline\hline \bm{c} & 0 & 1/2 & 1 & & \\ \hline \bm{b}\circ\bm{c} & 0 & 1/3 & 1/6 & 1/2 & \text{OK}\\ \hline\hline \bm{c}^{\circ2} & 0 & 1/4 & 1 & & \\ \hline \bm{b}\circ\bm{c}^{\circ2} & 0 & 1/6 & 1/6 & 1/3 & \text{OK}\\ \hline \bm{d} & -1/12 & 1/6 & 5/12 & & \\ \hline \bm{b}\circ\bm{d} & -1/72 & 1/9 & 5/72 & 1/6 & \text{OK}\\ \hline\hline \bm{c}^{\circ3} & 0 & 1/8 & 1 & & \\ \hline \bm{b}\circ\bm{c}^{\circ3} & 0 & 1/12 & 1/6 & 1/4 & \text{OK}\\ \hline (\bm{b}\circ\bm{c})\circ\bm{d} & 0 & 1/18 & 5/72 & 1/8 & \text{OK}\\ \hline \bm{e} & 1/6 & 1/3 & 0 & & \\ \hline \bm{e}\circ\bm{c}^{\circ2} & 0 & 1/12 & 0 & 1/12 & \text{OK}\\ \hline \bm{e}\circ\bm{d} & -1/72 & 1/18 & 0 & 1/24 & \text{OK}\\ \hline \end{array}

Order at least 4. Two features are worth noticing. The last column of \bm{A} vanishes, so the third stage does not feed back into the others: the method is not stiffly accurate in the IIIA sense. And \bm{d} has a negative entry, a reminder that the order conditions are algebraic identities with no sign structure to exploit.

7.3 Automated Verification

Checking a tableau by hand is useful once; after that it should be automated. The conditions of Equation 4.9 are a handful of inner products, so a few lines suffice to determine the order of any tableau — including the ones just examined, which serve as a regression test of the computations above.

rk_order <- function(A, b, tol = 1e-12) {
  A <- as.matrix(A); b <- as.numeric(b)
  s <- nrow(A); one <- rep(1, s)
  cc <- as.numeric(A %*% one)          # stage consistency defines c
  d  <- as.numeric(A %*% cc)
  e  <- as.numeric(t(A) %*% b)
  ok <- function(x, target) abs(x - target) < tol
  lev <- list(
    `1` = c(sum(b) - 1),
    `2` = c(sum(b * cc) - 1/2),
    `3` = c(sum(b * cc^2) - 1/3, sum(b * d) - 1/6),
    `4` = c(sum(b * cc^3) - 1/4, sum(b * cc * d) - 1/8,
            sum(e * cc^2) - 1/12, sum(e * d) - 1/24)
  )
  p <- 0
  for (k in names(lev)) {
    if (all(abs(lev[[k]]) < tol)) p <- as.integer(k) else break
  }
  p
}

g <- sqrt(3) / 6
tableaux <- list(
  "Explicit Euler"     = list(A = matrix(0, 1, 1), b = 1),
  "Explicit midpoint"  = list(A = rbind(c(0,0), c(1/2,0)), b = c(0,1)),
  "Ralston 2"          = list(A = rbind(c(0,0), c(2/3,0)), b = c(1/4,3/4)),
  "RK3"                = list(A = rbind(c(0,0,0), c(1/2,0,0), c(-1,2,0)),
                              b = c(1/6,2/3,1/6)),
  "RK4 (classical)"    = list(A = rbind(c(0,0,0,0), c(1/2,0,0,0),
                                        c(0,1/2,0,0), c(0,0,1,0)),
                              b = c(1/6,1/3,1/3,1/6)),
  "RK4 (3/8 rule)"     = list(A = rbind(c(0,0,0,0), c(1/3,0,0,0),
                                        c(-1/3,1,0,0), c(1,-1,1,0)),
                              b = c(1/8,3/8,3/8,1/8)),
  "Implicit midpoint"  = list(A = matrix(1/2, 1, 1), b = 1),
  "Crank--Nicolson"    = list(A = rbind(c(0,0), c(1/2,1/2)), b = c(1/2,1/2)),
  "Gauss--Legendre 4"  = list(A = rbind(c(1/4, 1/4 - g), c(1/4 + g, 1/4)),
                              b = c(1/2,1/2)),
  "Lobatto IIIA 4"     = list(A = rbind(c(0,0,0), c(5/24,1/3,-1/24),
                                        c(1/6,2/3,1/6)),
                              b = c(1/6,2/3,1/6)),
  "Lobatto IIIB 4"     = list(A = rbind(c(1/6,-1/6,0), c(1/6,1/3,0),
                                        c(1/6,5/6,0)),
                              b = c(1/6,2/3,1/6))
)

res <- data.frame(
  Method = names(tableaux),
  Stages = vapply(tableaux, function(m) nrow(as.matrix(m$A)), integer(1)),
  Order  = vapply(tableaux, function(m) rk_order(m$A, m$b), integer(1)),
  row.names = NULL
)
knitr::kable(res, align = "lrr",
             caption = "Order determined from the conditions of @eq-order-conditions-table (capped at 4).")

Order determined from the conditions of Equation 4.9 (capped at 4).
Method	Stages	Order
Explicit Euler	1	1
Explicit midpoint	2	2
Ralston 2	2	2
RK3	3	3
RK4 (classical)	4	4
RK4 (3/8 rule)	4	4
Implicit midpoint	1	2
Crank–Nicolson	2	2
Gauss–Legendre 4	2	4
Lobatto IIIA 4	3	4
Lobatto IIIB 4	3	4

The column is capped at four because only the conditions up to that level have been implemented; for the two Lobatto methods and for Gauss–Legendre the true order is four, while a fifth-order test would be needed to confirm it rather than assume it.

What this check does and does not prove

The conditions verify the order of consistency of the tableau, that is, the rate at which the local truncation error vanishes. They say nothing about:

stability — a high-order tableau may have a tiny absolute stability region and be useless on a stiff problem (compare the 3/8 rule with Gauss–Legendre);
the error constant — two methods of the same order can differ by an order of magnitude in actual error, which is exactly what distinguishes Ralston’s method from the explicit midpoint method;
implementation cost — an implicit tableau of order four requires solving a nonlinear system of dimension sd at every step.

Order is a necessary property, never a sufficient criterion for choosing a method.

Exercises

Verify that the classical RK4 tableau satisfies all eight conditions of Equation 4.9, and identify which ones fail at order five.
Show that no explicit two-stage method can have order three, by writing the general tableau with c_2=\alpha and deriving a contradiction between \bm{b}\cdot\bm{c}^{\circ2}=\tfrac13 and \bm{b}\cdot\bm{d}=\tfrac16.
Heun’s third-order method has \bm{c}=(0,\tfrac13,\tfrac23)^T and \bm{b}=(\tfrac14,0,\tfrac34)^T. Determine \bm{A} from the order conditions and confirm the order with rk_order.
For the Lobatto IIIA method of order four, verify that the last row of \bm{A} equals \bm{b}^T (the method is stiffly accurate) and explain why this implies R(\infty)=0.
Modify rk_order to include the nine fifth-order conditions Equation 4.10 and recompute the table. Which of the methods listed actually attain order five?

Summary

The order of a Runge–Kutta method follows from a finite set of algebraic identities among \bm{A}, \bm{b} and \bm{c} — no Taylor expansion and no numerical experiment is required.
Two auxiliary vectors, \bm{d}=\bm{A}\bm{c} and \bm{e}=\bm{A}^T\bm{b}, express all the conditions through order four.
The conditions are cumulative and must be checked level by level: satisfying some conditions of a level, as Ralston’s method and RK3 do, does not raise the order.
Explicit s-stage methods have order at most s, and from order five on they need strictly more stages than the order (Table 4.3); implicit s-stage methods reach at most 2s, attained by the Gauss–Legendre family.