In Chapter 5 we discovered that demanding second-order accuracy does not fix a unique method but instead produces a one-parameter family. This chapter asks the natural follow-up: what happens when we demand third-order accuracy? The same structure reappears — the conditions can still be satisfied in multiple ways — but the algebra is substantially harder. More importantly, examining how the number of constraints grows at each order reveals why fourth-order Runge–Kutta is the last method derivable by hand without combinatorial machinery.
A three-stage Runge–Kutta method evaluates \(f\) at three points per step. Following the unifying structure from Section 4.1, the update formula is:
\[y_{n+1} = y_n + h \cdot (b_1 \cdot k_1 + b_2 \cdot k_2 + b_3 \cdot k_3)\]
where the three stage slopes are:
\[k_1 = f(x_n,\ y_n)\]
\[k_2 = f(x_n + c_2 \cdot h,\ y_n + a_{21} \cdot h \cdot k_1)\]
\[k_3 = f(x_n + c_3 \cdot h,\ y_n + h \cdot (a_{31} \cdot k_1 + a_{32} \cdot k_2))\]
The eight free coefficients fall into three groups: \(b_1\), \(b_2\), \(b_3\) are the weights in the update formula; \(c_2\), \(c_3\) are the stage positions along the \(x\)-axis (with \(c_1 = 0\) always); and \(a_{21}\), \(a_{31}\), \(a_{32}\) are the coupling coefficients that determine how each stage depends on previous ones. For an explicit method the coupling matrix is strictly lower triangular: \(k_1\) has no coupling, \(k_2\) depends only on \(k_1\), and \(k_3\) depends on both \(k_1\) and \(k_2\).
Two of the eight coefficients are not independently free. The consistency conditions \(c_2 = a_{21}\) and \(c_3 = a_{31} + a_{32}\) require that each stage position equal the row sum of its coupling coefficients. These conditions place each evaluation point at the \(x\)-coordinate it nominally claims, and they reduce the effective free parameter count from eight to six.
Compare this to the two-stage scheme from Chapter 5, which had four parameters and one degree of freedom after applying order conditions. Adding a third stage doubles the parameter count, but also adds new constraints. The order conditions will narrow the six effective free parameters to just two.
For the three-stage scheme to be third-order accurate, the numerical increment \(h \cdot (b_1 k_1 + b_2 k_2 + b_3 k_3)\) must match the exact Taylor expansion of \(y(x_n + h)\) through and including the \(h^3\) term. The derivation follows the same path as Chapter 5: expand \(k_2\) and \(k_3\) using the two-variable Taylor series from Chapter 3, form the full increment, and compare it term by term with the exact series. The key difference is that \(k_3\) depends on \(k_2\), and this coupling generates a fourth independent condition that has no analog in the two-stage case.
Matching the coefficient of \(h\) gives the first order condition:
\[b_1 + b_2 + b_3 = 1\]
Matching the coefficient of \(h^2\) yields one condition on the weighted stage positions:
\[b_2 \cdot c_2 + b_3 \cdot c_3 = \tfrac{1}{2}\]
At order \(h^3\), the exact Taylor derivative \(y''' = f_{xx} + 2 f \cdot f_{xy} + f_x \cdot f_y + f \cdot (f_y)^2 + f^2 \cdot f_{yy}\) contains five terms, but when the numerical increment is expanded and compared, these five terms collapse into only two independent constraints. The first constraint comes from matching the \(f_{xx}\), \(f \cdot f_{xy}\), and \(f^2 \cdot f_{yy}\) terms (all of which share the same coefficient in both the exact and numerical expansions). The second comes from matching the \(f_x \cdot f_y\) and \(f \cdot (f_y)^2\) terms, which involve the coupling coefficient \(a_{32}\).
Condition 3 (from the \(f_{xx}\), \(f \cdot f_{xy}\), and \(f^2 \cdot f_{yy}\) terms):
\[b_2 \cdot c_2^2 + b_3 \cdot c_3^2 = \tfrac{1}{3}\]
Condition 4 (from the \(f_x \cdot f_y\) and \(f \cdot (f_y)^2\) terms):
\[b_3 \cdot a_{32} \cdot c_2 = \tfrac{1}{6}\]
Condition 4 involves only \(a_{32}\) and \(c_2\) — not \(a_{31}\). The reason is that \(a_{31}\) tracks how \(k_3\) depends on \(k_1\), which is evaluated at the base point \((x_n, y_n)\) and carries no displacement information. The cross-derivative structure arises exclusively from how \(k_3\) inherits the \(x\)-displacement of \(k_2\) through the coupling \(a_{32}\): the term \(a_{32} \cdot c_2\) measures how far \(k_2\) was displaced in \(x\), and this displacement is what produces the mixed second-order partial derivative structure at order \(h^3\).
All four conditions are independent — satisfying the first two does not force the last two. The Python demonstration below constructs a method that satisfies conditions 1 and 2 exactly while conditions 3 and 4 both fail. Note that \(f(x, y) = y\) (used in later sections) has \(f_{xx} = f_{xy} = f_{yy} = 0\), which would hide conditions 2 and 3; the test function below uses \(f = x + y\) so all four conditions are visible.
from fractions import Fraction
b1t, b2t, b3t = Fraction(1, 4), Fraction(1, 2), Fraction(1, 4)
c2t, c3t = Fraction(1, 2), Fraction(1, 1)
a21t, a31t, a32t = Fraction(1, 2), Fraction(0, 1), Fraction(1, 1)
[c2t == a21t, c3t == a31t + a32t][True, True][b1t + b2t + b3t == 1,
b2t*c2t + b3t*c3t == Fraction(1, 2),
b2t*c2t**2 + b3t*c3t**2 == Fraction(1, 3),
b3t*a32t*c2t == Fraction(1, 6)][True, True, False, False]The first list shows consistency holds; the second list shows the first two order conditions are met but the last two fail.
The phrase “solution space shrinks” requires care. Adding stages always adds more unknowns, so the degrees of freedom do not literally decrease. The meaningful measure is the ratio of order conditions to stages: this ratio grows steeply, and when it exceeds the effective free parameters per stage, achieving the next order requires more stages than the order number. That is the Butcher barrier — encountered in Chapter 9.
The table below shows the condition count, minimum stages required, and the ratio:
import pandas as pd
cond_counts = [1, 2, 4, 8, 17, 37]
min_stages = [1, 2, 3, 4, 6, 7]
orders = [1, 2, 3, 4, 5, 6]
ratios = [round(c / s, 3) for c, s in zip(cond_counts, min_stages)]
pd.DataFrame({
"Order p": orders,
"Conditions": cond_counts,
"Min stages s": min_stages,
"Conds/s": ratios,
})| Order p | Conditions | Min stages s | Conds/s | |
|---|---|---|---|---|
| 0 | 1 | 1 | 1 | 1.000 |
| 1 | 2 | 2 | 2 | 1.000 |
| 2 | 3 | 4 | 3 | 1.333 |
| 3 | 4 | 8 | 4 | 2.000 |
| 4 | 5 | 17 | 6 | 2.833 |
| 5 | 6 | 37 | 7 | 5.286 |
Reading the table: the ratio of conditions to minimum stages grows from \(1.0\) at order 2, to \(1.33\) at order 3, to \(2.0\) at order 4, to \(2.83\) at order 5. The jump between orders 4 and 5 is especially sharp: a five-stage method provides only about ten effective free parameters after applying consistency, but order 5 requires 17 conditions. This shortfall forces the use of six stages to achieve fifth-order accuracy, even though five stages would seem sufficient by naive counting. Chapter 9 explains why.
For the methods in this book: RK3 (three stages) has 8 total parameters minus 2 consistency conditions minus 4 order conditions = 2 degrees of freedom. RK4 (four stages) will have 13 minus 3 minus 8 = 2 degrees of freedom. The absolute count stays at 2, but the fraction of parameters consumed by constraints grows: \(6/8 = 75\%\) for RK3, \(11/13 = 85\%\) for RK4. The sense in which the solution “space” shrinks is that it occupies a smaller and smaller fraction of the total parameter space as the order rises.
We have four order conditions and six effective free parameters (after applying consistency). Treating \(c_2\) and \(c_3\) as the two free parameters, conditions 2 and 3 form a \(2 \times 2\) linear system in \(b_2\) and \(b_3\) (for given \(c_2\) and \(c_3\)). Condition 1 then gives \(b_1 = 1 - b_2 - b_3\), and condition 4 gives \(a_{32} = 1/(6 \cdot b_3 \cdot c_2)\). The remaining coupling coefficients follow from consistency: \(a_{21} = c_2\) and \(a_{31} = c_3 - a_{32}\).
SymPy solves all four conditions simultaneously:
import sympy as sp
b1, b2, b3, c2, c3, a32 = sp.symbols('b1 b2 b3 c2 c3 a32')
sp.solve(
[b1 + b2 + b3 - 1,
b2*c2 + b3*c3 - sp.Rational(1, 2),
b2*c2**2 + b3*c3**2 - sp.Rational(1, 3),
b3*a32*c2 - sp.Rational(1, 6)],
[b1, b2, b3, a32])[((6*c2*c3 - 3*c2 - 3*c3 + 2)/(6*c2*c3),
-(3*c3 - 2)/(6*c2*(c2 - c3)),
(3*c2 - 2)/(6*c3*(c2 - c3)),
c3*(c2 - c3)/(c2*(3*c2 - 2)))]The output gives \(b_2\), \(b_3\), \(a_{32}\) as rational functions of \(c_2\) and \(c_3\), and \(b_1 = 1 - b_2 - b_3\). The expressions are valid whenever \(c_2\) and \(c_3\) are both nonzero and unequal to each other (\(c_2 = c_3\) would make the \(2 \times 2\) linear system singular, because both stages would lie at the same \(x\)-position). Every valid choice of \(c_2\), \(c_3\) with these restrictions yields a distinct but equally valid third-order method.
sp.solve(
[b1 + b2 + b3 - 1,
b2*sp.Rational(1, 2) + b3*1 - sp.Rational(1, 2),
b2*sp.Rational(1, 2)**2 + b3*1**2 - sp.Rational(1, 3),
b3*a32*sp.Rational(1, 2) - sp.Rational(1, 6)],
[b1, b2, b3, a32])[(1/6, 2/3, 1/6, 2)]The output confirms \(b_1 = 1/6\), \(b_2 = 2/3\), \(b_3 = 1/6\), and \(a_{32} = 2\). Applying consistency: \(a_{21} = c_2 = 1/2\) and \(a_{31} = c_3 - a_{32} = 1 - 2 = -1\). These are exactly the coefficients of Kutta’s third-order method from 1901.
Wilhelm Kutta derived this method in his 1901 paper as the natural three-stage analog of Simpson’s rule. Setting \(c_2 = 1/2\) and \(c_3 = 1\) places the three evaluation points at the left endpoint, the midpoint, and the right endpoint of each step interval — the same three nodes used by Simpson’s rule for numerical integration. The update weights \(1/6\), \(2/3\), \(1/6\) are exactly Simpson’s weights.
The full specification is:
\[k_1 = f(x_n,\ y_n)\]
\[k_2 = f\!\left(x_n + \tfrac{h}{2},\ y_n + \tfrac{h}{2} \cdot k_1\right)\]
\[k_3 = f\!\left(x_n + h,\ y_n + h \cdot (-k_1 + 2 \cdot k_2)\right)\]
\[y_{n+1} = y_n + \frac{h}{6}\,(k_1 + 4 \cdot k_2 + k_3)\]
Several features deserve comment. First, the \(a_{31} = -1\) coupling means that the \(y\)-argument of \(k_3\) is \(y_n - h \cdot k_1 + 2 h \cdot k_2\). This unusual negative coefficient is not a flaw: it is a direct consequence of forcing the third stage to \(x_n + h\) while simultaneously satisfying all four order conditions. Second, the update formula \(\tfrac{h}{6}(k_1 + 4 k_2 + k_3)\) matches Simpson’s \(1/3\) rule: the midpoint slope \(k_2\) receives four times the weight of the two endpoint slopes, which is exactly the weighting that makes numerical integration exact for polynomials up to degree three.
from fractions import Fraction
b1v, b2v, b3v = Fraction(1, 6), Fraction(2, 3), Fraction(1, 6)
c2v, c3v = Fraction(1, 2), Fraction(1, 1)
a21v, a31v, a32v = Fraction(1, 2), Fraction(-1, 1), Fraction(2, 1)
[c2v == a21v, c3v == a31v + a32v][True, True][b1v + b2v + b3v == 1,
b2v*c2v + b3v*c3v == Fraction(1, 2),
b2v*c2v**2 + b3v*c3v**2 == Fraction(1, 3),
b3v*a32v*c2v == Fraction(1, 6)][True, True, True, True]All six should return True. The method is third-order accurate and is still widely used today as a lightweight alternative to RK4 in applications where one fewer function evaluation per step matters.
We implement Kutta’s third-order method and compare it against Euler’s method and the Midpoint Method on the test problem \(dy/dx = y\), \(y(0) = 1\), whose exact solution is \(e^x\).
import numpy as np
def rk3_step(fx, x, y, h):
k1 = fx(x, y)
k2 = fx(x + h/2, y + (h/2)*k1)
k3 = fx(x + h, y + h*(-k1 + 2*k2))
return y + (h/6)*(k1 + 4*k2 + k3)
def rk3_method(fx, x0, y0, x_end, h):
xs = [x0]; ys = [y0]
while xs[-1] < x_end - h/2:
ys.append(rk3_step(fx, xs[-1], ys[-1], h))
xs.append(xs[-1] + h)
return np.array(xs), np.array(ys)def f(x, y):
return y
def euler_step(fx, x, y, h):
return y + h * fx(x, y)
def euler_method(fx, x0, y0, x_end, h):
xs = [x0]; ys = [y0]
while xs[-1] < x_end - h/2:
ys.append(euler_step(fx, xs[-1], ys[-1], h))
xs.append(xs[-1] + h)
return np.array(xs), np.array(ys)
def midpoint_step(fx, x, y, h):
k1 = fx(x, y)
k2 = fx(x + h/2, y + (h/2)*k1)
return y + h*k2
def midpoint_method(fx, x0, y0, x_end, h):
xs = [x0]; ys = [y0]
while xs[-1] < x_end - h/2:
ys.append(midpoint_step(fx, xs[-1], ys[-1], h))
xs.append(xs[-1] + h)
return np.array(xs), np.array(ys)import matplotlib.pyplot as plt
h_cmp = 0.25
xs_e, ys_e = euler_method(f, 0, 1, 2, h_cmp)
xs_m, ys_m = midpoint_method(f, 0, 1, 2, h_cmp)
xs_3, ys_3 = rk3_method(f, 0, 1, 2, h_cmp)
xc = np.linspace(0, 2, 200)
fig, ax = plt.subplots(figsize=(6.5, 4.2))
ax.plot(xc, np.exp(xc), color='gray', linestyle='--', label='Exact: $e^x$')
ax.plot(xs_e, ys_e, color='red', marker='o', markersize=4, label='Euler')
ax.plot(xs_m, ys_m, color='blue', marker='s', markersize=4, label='Midpoint')
ax.plot(xs_3, ys_3, color='purple', marker='^', markersize=4, label='RK3')
ax.set_xlabel('x'); ax.set_ylabel('y')
ax.set_title('Euler (red), Midpoint (blue), RK3 (purple), h = 0.25')
ax.legend()
fig.tight_layout()
fig.savefig("06_RK3/rk3-comparison.pdf")
fig.savefig("06_RK3/rk3-comparison.svg")
plt.close(fig)
import pandas as pd
h_values = [0.5, 0.25, 0.125, 0.0625, 0.03125]
eur_errors = [abs(euler_method(f, 0, 1, 2, hv)[1][-1] - np.exp(2))
for hv in h_values]
mid_errors = [abs(midpoint_method(f, 0, 1, 2, hv)[1][-1] - np.exp(2))
for hv in h_values]
rk3_errors = [abs(rk3_method(f, 0, 1, 2, hv)[1][-1] - np.exp(2))
for hv in h_values]
def ratios(errs):
out = ["---"]
for i in range(1, len(errs)):
out.append(f"{errs[i]/errs[i-1]:.4f}")
return out
pd.DataFrame({
"h": h_values,
"Euler": [f"{e:.4e}" for e in eur_errors],
"ratio (E)": ratios(eur_errors),
"Midpoint": [f"{e:.4e}" for e in mid_errors],
"ratio (M)": ratios(mid_errors),
"RK3": [f"{e:.4e}" for e in rk3_errors],
"ratio (R3)": ratios(rk3_errors),
})| h | Euler | ratio (E) | Midpoint | ratio (M) | RK3 | ratio (R3) | |
|---|---|---|---|---|---|---|---|
| 0 | 0.50000 | 2.3266e+00 | --- | 4.1616e-01 | --- | 5.1635e-02 | --- |
| 1 | 0.25000 | 1.4286e+00 | 0.6140 | 1.2681e-01 | 0.3047 | 7.8801e-03 | 0.1526 |
| 2 | 0.12500 | 8.0581e-01 | 0.5641 | 3.4973e-02 | 0.2758 | 1.0884e-03 | 0.1381 |
| 3 | 0.06250 | 4.3039e-01 | 0.5341 | 9.1757e-03 | 0.2624 | 1.4301e-04 | 0.1314 |
| 4 | 0.03125 | 2.2278e-01 | 0.5176 | 2.3492e-03 | 0.2560 | 1.8328e-05 | 0.1282 |
The ratio columns confirm the order of each method. Euler’s ratios approach \(0.5\) (first-order), Midpoint’s ratios approach \(0.25\) (second-order), and RK3’s ratios approach \(0.125\) (third-order). Halving the step size reduces the RK3 error by a factor of \(2^3 = 8\).
The Butcher tableau introduced in Section 5.10 organizes all coefficients of an RK method in a single array. For a general three-stage method the tableau has the form:
\[\begin{array}{c|ccc} 0 & & & \\ c_2 & a_{21} & & \\ c_3 & a_{31} & a_{32} & \\ \hline & b_1 & b_2 & b_3 \end{array}\]
For Kutta’s third-order method, substituting \(c_2 = 1/2\), \(c_3 = 1\), \(a_{21} = 1/2\), \(a_{31} = -1\), \(a_{32} = 2\), \(b_1 = b_3 = 1/6\), \(b_2 = 2/3\):
\[\begin{array}{c|ccc} 0 & & & \\ \tfrac{1}{2} & \tfrac{1}{2} & & \\ 1 & -1 & 2 & \\ \hline & \tfrac{1}{6} & \tfrac{2}{3} & \tfrac{1}{6} \end{array}\]
Three patterns visible in the tableau will continue into RK4. First, the \(c\) column gains one entry per stage. Second, the lower-triangular coupling body grows: RK2 had one entry (\(a_{21}\)), RK3 has three (\(a_{21}\), \(a_{31}\), \(a_{32}\)), and RK4 will have six. Third, the weight row gains one entry per stage. The negative \(a_{31} = -1\) is unusual but not exceptional: it is the necessary consequence of placing the third stage at \(x_n + h\) while satisfying all four order conditions.
The growth of the coupling body is the combinatorial source of the condition-count explosion at higher orders. In an explicit \(s\)-stage method the coupling body has \(s(s-1)/2\) entries: \(1\) for RK2, \(3\) for RK3, \(6\) for RK4. Each entry is a free parameter, but each also participates in the order conditions in new combinations. Chapter 9 connects this to the theory of rooted trees, which gives a precise count of how many independent conditions arise at each order.
A four-stage scheme adds five new coefficients compared to three stages: \(b_4\), \(c_4\), \(a_{41}\), \(a_{42}\), \(a_{43}\). The total parameter count rises from 8 to 13. Three consistency conditions (\(c_2 = a_{21}\), \(c_3 = a_{31} + a_{32}\), \(c_4 = a_{41} + a_{42} + a_{43}\)) reduce this to ten effective free parameters.
At order 4 there are eight order conditions in total. The four conditions inherited from order 3 (conditions 1 through 4 above) still apply. The new conditions at \(h^4\) come from matching four structurally independent derivative groups in the exact \(y''''\) expansion, and they involve the full coupling body \(a_{31}\), \(a_{32}\), \(a_{41}\), \(a_{42}\), \(a_{43}\). The system has ten free parameters and eleven constraints (three consistency plus eight order), leaving exactly two degrees of freedom.
import pandas as pd
stage_count = [2, 3, 4]
total_params = [s + (s - 1) + s*(s - 1)//2 for s in stage_count]
consistency = [s - 1 for s in stage_count]
order_conds = [3, 4, 8]
dof = [tp - co - oc for tp, co, oc
in zip(total_params, consistency, order_conds)]
pd.DataFrame({
"Stages": stage_count,
"Total params": total_params,
"Consistency": consistency,
"Order conds": order_conds,
"DOF": dof,
})| Stages | Total params | Consistency | Order conds | DOF | |
|---|---|---|---|---|---|
| 0 | 2 | 4 | 1 | 3 | 0 |
| 1 | 3 | 8 | 2 | 4 | 2 |
| 2 | 4 | 13 | 3 | 8 | 2 |
The table confirms: 1 degree of freedom for RK2, 2 for RK3, and 2 for RK4. The two degrees of freedom in RK4 correspond to two genuine choices in the classical derivation: one is related to the time-symmetry of the method (choosing between the classical RK4 and the “3/8-rule” variant), and the other is an internal gauge. Chapter 7 derives all eight order conditions, identifies both free choices, and shows how the classical coefficients \(1/6, 1/3, 1/3, 1/6\) emerge as the canonical solution.