This chapter completes the derivation of Runge–Kutta methods begun in Chapters 5 and 6. In those chapters we derived the order conditions for two-stage and three-stage schemes, observed that the number of constraints grows with each additional order, and established that four-stage methods inherit all prior conditions while adding four new ones at order \(h^4\). Here we carry out the complete derivation: we expand all four stages symbolically in Python with SymPy, extract the eight order conditions, solve the system, and confirm that the classical coefficients \(1/6, 1/3, 1/3, 1/6\) satisfy all constraints exactly.
Every step uses only the two-variable Taylor series from Chapter 3 and the algebraic pipeline introduced in Chapters 5 and 6. No facts are borrowed from outside that framework. At the end of this chapter we will have derived the classical RK4 method from first principles, and we will understand why four stages is the last order at which this derivation is tractable by hand.
7.1 Setting Up the Four-Stage Scheme
Following the unifying structure first introduced in Section 4.1, a four-stage explicit Runge–Kutta method advances the solution by one step using a weighted combination of four slope estimates:
The scheme has thirteen coefficients in total: four weights (\(b_1\), \(b_2\), \(b_3\), \(b_4\)), three stage positions (\(c_2\), \(c_3\), \(c_4\) with \(c_1 = 0\) by convention), and six coupling coefficients (\(a_{21}\), \(a_{31}\), \(a_{32}\), \(a_{41}\), \(a_{42}\), \(a_{43}\)). The stage positions are not free: the consistency conditions \(c_2 = a_{21}\), \(c_3 = a_{31} + a_{32}\), \(c_4 = a_{41} + a_{42} + a_{43}\) reduce the effective count to ten. As shown in the table below, eight order conditions then leave exactly two degrees of freedom.
import pandas as pdstage_count = [2, 3, 4]total_params = [s + (s -1) + s*(s -1)//2for s in stage_count]consistency = [s -1for s in stage_count]order_conds_count = [3, 4, 8]dof = [tp - co - oc for tp, co, ocinzip(total_params, consistency, order_conds_count)]pd.DataFrame({"Stages": stage_count,"Total params": total_params,"Consistency": consistency,"Order conds": order_conds_count,"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
7.2 Letting SymPy Expand Each Stage Symbolically
The central computation in this chapter is expanding \(k_2\), \(k_3\), and \(k_4\) as power series in \(h\) around the base point \((x_n, y_n)\). Each stage evaluates \(f\) at a displaced point, and the chain rule turns that evaluation into a sum of partial derivatives of \(f\) weighted by the displacements. The complication is that \(k_3\) depends on \(k_2\), and \(k_4\) depends on \(k_2\) and \(k_3\), so the expansions are nested. The code below defines a helper fT that encodes the two-variable Taylor expansion from Chapter 3 to third order, then applies it to each stage in sequence, truncating at \(h^3\) after each substitution to control growth.
import sympy as spxn, yn, h = sp.symbols('xn yn h')c2, c3, c4 = sp.symbols('c2 c3 c4')a21, a31, a32, a41, a42, a43 = sp.symbols('a21 a31 a32 a41 a42 a43')b1, b2, b3, b4 = sp.symbols('b1 b2 b3 b4')f = sp.Function('f')def fT(a, b):"""Two-variable Taylor expansion of f at (xn+a, yn+b) around (xn, yn), to third order in the displacements.""" f0 = f(xn, yn)return (f0+ a*f0.diff(xn) + b*f0.diff(yn)+ (a**2/2)*f0.diff(xn, 2)+ a*b*f0.diff(xn).diff(yn)+ (b**2/2)*f0.diff(yn, 2)+ (a**3/6)*f0.diff(xn, 3)+ (a**2*b/2)*f0.diff(xn, 2).diff(yn)+ (a*b**2/2)*f0.diff(xn).diff(yn, 2)+ (b**3/6)*f0.diff(yn, 3))def truncate_h(expr, n):"""Keep only terms up to and including h^n.""" expr = sp.expand(expr)returnsum(expr.coeff(h, k) * h**k for k inrange(n +1))k1 = f(xn, yn)k2 = truncate_h(fT(c2*h, a21*h*k1), 3)k3 = truncate_h(fT(c3*h, h*(a31*k1 + a32*k2)), 3)k4 = truncate_h(fT(c4*h, h*(a41*k1 + a42*k2 + a43*k3)), 3)# truncate_h already groups by h, so k2 displays as a polynomial in h.k2
The \(h^0\) term of every \(k_i\) is \(f(x_n, y_n)\), confirming that all four stages reduce to the same base evaluation as \(h\) tends to zero. The \(h^1\) term of \(k_2\) is \(c_2 f_x + a_{21} f_0 f_y\) — the directional derivative of \(f\) in the \((c_2, a_{21} f_0)\) direction. The \(h^2\) and \(h^3\) terms of \(k_3\) and \(k_4\) inherit coupling coefficients \(a_{32}\), \(a_{42}\), \(a_{43}\) through the nested substitutions. This is the mechanism by which later stages “know” about earlier ones.
7.3 Extracting Coefficients by Power of \(h\) — A Symbolic Pipeline
With \(k_1\) through \(k_4\) expanded in \(h\), the full numerical increment \(h\,(b_1 k_1 + b_2 k_2 + b_3 k_3 + b_4 k_4)\) is a power series in \(h\) whose coefficients involve the \(b\), \(c\), and \(a\) parameters together with partial derivatives of \(f\). For a fourth-order method every coefficient of \(h^p\) in that series must equal the corresponding coefficient from the exact Taylor expansion of \(y(x_n + h) - y_n\). The pipeline below forms the increment, truncates it at fifth order in \(h\), and extracts the symbolic coefficient at each power. Matching those coefficients to the exact derivatives produces the eight order conditions.
The \(h^1\) coefficient equals \((b_1 + b_2 + b_3 + b_4)\, f(x_n, y_n)\). For this to match the exact first-order term \(f(x_n, y_n)\) we need \(b_1 + b_2 + b_3 + b_4 = 1\), which is the first order condition. The \(h^2\) coefficient decomposes into terms proportional to the partial derivatives \(f_x(x_n, y_n)\) and \(f_y(x_n, y_n)\); matching their coefficients to the exact \((1/2) y''\) produces the second condition. The \(h^3\) and \(h^4\) levels each yield two and four independent conditions respectively, giving eight in total. The next section records all eight.
7.4 The Eight Order Conditions — Grouped by Power of \(h\)
Matching the numerical increment to the exact Taylor expansion of \(y\) term by term produces eight conditions on the ten effective free parameters. The conditions are grouped below by the power of \(h\) they arise from. Conditions 1 and 2 are identical to those for RK2 and RK3 — they are inherited at every higher order. Conditions 3 and 4 first appear at RK3. Conditions 5 through 8 are new at RK4.
Order \(h\) — one condition (the weights must sum to one):
\[b_1 + b_2 + b_3 + b_4 = 1\]
Order \(h^2\) — one condition (the weighted stage positions must equal \(1/2\)):
The eight order conditions involve ten effective free parameters (after the three consistency substitutions eliminate \(a_{21}\), \(a_{31}\), and \(a_{41}\)). Eight conditions in ten unknowns leave exactly two degrees of freedom. The next section shows how two structural choices reduce the system to a unique solution.
7.5 Solving the System in SymPy
The eight order conditions plus three consistency conditions form a system of eleven equations in thirteen unknowns, leaving two degrees of freedom. The classical RK4 method corresponds to a specific pair of structural choices: place the two intermediate stages at the interval midpoint (\(c_2 = c_3 = 1/2\)) and the final stage at the endpoint (\(c_4 = 1\)), and zero out the coupling from stage 1 into stage 3 and from stage 2 into stage 4 (\(a_{31} = 0\), \(a_{42} = 0\)). These choices create the symmetric “ladder” structure that makes RK4 elegant and geometrically transparent.
SymPy returns a single solution: \(a_{21} = 1/2\), \(a_{32} = 1/2\), \(a_{41} = 0\), \(a_{43} = 1\), \(b_1 = b_4 = 1/6\), \(b_2 = b_3 = 1/3\). These are the classical RK4 coefficients. The two free structural choices (midpoint staging and the zero coupling entries) fix the last two degrees of freedom, leaving no ambiguity.
All eight conditions evaluate to True. The classical RK4 method is not the only fourth-order method — any other choice of \(c_2\), \(c_3\) (subject to certain constraints) generates a valid fourth-order scheme — but the classical coefficients have the simplest algebraic form and an elegant geometric interpretation that the following sections develop.
7.6 The Classical RK4 Coefficients
The SymPy solve call returns a unique solution: \(a_{21} = 1/2\), \(a_{32} = 1/2\), \(a_{41} = 0\), \(a_{43} = 1\), \(b_1 = b_4 = 1/6\), \(b_2 = b_3 = 1/3\). Together with the fixed parameters \(c_2 = c_3 = 1/2\), \(c_4 = 1\), \(a_{31} = 0\), \(a_{42} = 0\), these define the classical Runge–Kutta 4th-order method. Because \(a_{31} = 0\) the third stage does not depend on the first stage — only on the second. Because \(a_{41} = 0\) and \(a_{42} = 0\) the fourth stage depends only on the third.
The classical RK4 formula has a striking symmetry: \(k_2\) and \(k_3\) are structurally identical in that both evaluate \(f\) at the midpoint \(x_n + h/2\), yet they use different \(y\)-arguments. This “leapfrog” structure is a direct consequence of the zero coupling entries \(a_{31} = 0\) and \(a_{42} = 0\) that we chose to fix the degrees of freedom.
7.7 Symbolic Verification
Having derived the classical coefficients algebraically, we verify them computationally in two steps. First, all eight order conditions evaluate to True when the classical values are substituted. Second, a convergence table on the test problem \(dy/dx = y\), \(y(0) = 1\) shows that the global error at \(x = 2\) scales as \(h^4\): halving \(h\) reduces the error by a factor of \(16\).
import pandas as pddef rk4_run(h): x_n =0.0; y_n =1.0 n_steps =round(2/ h)for _ inrange(n_steps): y_n = rk4_step(f_test, x_n, y_n, h) x_n += hreturnabs(y_n - exact_solution(2))h_values = [0.4, 0.2, 0.1, 0.05, 0.025]errors = [rk4_run(h) for h in h_values]pd.DataFrame({"h": h_values,"error at x=2": [f"{e:.4e}"for e in errors],})
h
error at x=2
0
0.400
2.2624e-03
1
0.200
1.6686e-04
2
0.100
1.1332e-05
3
0.050
7.3830e-07
4
0.025
4.7114e-08
Each time \(h\) is halved the error decreases by approximately a factor of \(16 = 2^4\), confirming fourth-order convergence. This is the defining property of a fourth-order method: one extra digit of accuracy for roughly the same number of function evaluations per unit interval.
7.8 Geometric Interpretation
The four \(k\)-values correspond to four evaluations of \(f\) at locations spread across \([x_n, x_n + h]\). Stage positions \(c_1 = 0\), \(c_2 = c_3 = 1/2\), \(c_4 = 1\) place \(k_1\) at the left endpoint, \(k_2\) and \(k_3\) at the midpoint, and \(k_4\) at the right endpoint. The \(y\)-coordinate at each sampling point is a linear combination of previous \(k\)-values, so the four points trace a path that adapts to the local curvature of the solution.
Four sampling locations for RK4 along the solution curve \(y = e^x\). Red is \(k_1\), blue is \(k_2\), green is \(k_3\), orange is \(k_4\).
The two midpoint points (blue and green) are distinct because \(k_3\) uses a \(y\)-value corrected by \(k_2\): \(k_3\) lands above \(k_2\) because \(k_2\) overshot the curve, and \(k_3\) corrects upward toward the true solution.
7.9 Understanding the Weights: A Simpson’s Rule Connection
The weights \(b_1 = b_4 = 1/6\) and \(b_2 = b_3 = 1/3\) are not arbitrary. Simpson’s composite rule for integrating over \([x_n, x_n + h]\) assigns weight \(1/6\) to each endpoint and \(2/3\) to the midpoint. Because RK4 evaluates \(f\) twice at the midpoint (\(k_2\) and \(k_3\)), each midpoint evaluation receives half the Simpson midpoint weight: \((2/3)/2 = 1/3\). Advancing a solution of \(dy/dx = f(x, y)\) over one step is equivalent to numerically integrating \(f\) along the solution curve, and the RK4 weights are the Simpson quadrature weights for that integral.
The combined midpoint weight is \(1/3 + 1/3 = 2/3\), matching the Simpson coefficient for the midpoint. The two endpoint weights \(1/6\) each match the Simpson endpoint coefficients. No other four-point quadrature rule with these stage positions achieves fourth-order accuracy, so the classical RK4 weights are uniquely determined by the Simpson correspondence.
7.10 Why RK4 Is the Last Method You Can Derive by Hand
The number of order conditions grows rapidly with the target order. At order \(p = 4\) there are \(8\) conditions, which is exactly what we solved in this chapter using SymPy. At order \(p = 5\) the count jumps to \(17\), and the minimum number of stages required increases to \(6\) — a fifth-order method cannot be achieved with just five stages. This is Butcher’s fundamental result from 1963. Beyond order 4 no derivation by hand is feasible.
Each order condition corresponds to one rooted tree of the appropriate order. The number of independent derivative structures at order \(p\) equals the number of rooted trees with \(p\) nodes — a combinatorial count that grows super-exponentially. Chapter 9 explores this connection and shows how SymPy can still derive methods at orders 5 and 6 by solving the expanded systems symbolically.
7.11 Butcher Tableau for Classical RK4
The Butcher tableau encodes all coefficients of a Runge–Kutta method in a compact matrix form. For a four-stage method the tableau has four \(c\)-values (stage positions), a \(4 \times 4\) coefficient matrix \(A\) whose lower-triangular entries drive the stage evaluations, and a row of four \(b\)-weights. The general four-stage tableau in symbolic form is:
Reading the tableau row by row recovers the \(k\)-stage formulas. Row 1 says \(k_1\) evaluates at \((x_n, y_n)\) with no prior stages. Row 2 says \(k_2\) uses a half-step with coefficient \(1/2\) from \(k_1\). Row 3 says \(k_3\) uses a half-step from \(k_2\) only (the zero in position \(a_{31}\) means \(k_1\) does not appear). Row 4 says \(k_4\) uses a full step from \(k_3\) only. The \(b\)-row averages the four slopes with Simpson weights \(1/6, 1/3, 1/3, 1/6\).
This completes the derivation of the classical Runge–Kutta 4th-order method. Starting from the general four-stage scheme in Section 7.1, we used the two-variable Taylor expansion from Chapter 3 to expand all four stages symbolically, extracted eight order conditions by matching Taylor coefficients, applied two structural constraints to fix the remaining degrees of freedom, and arrived at the unique classical RK4 coefficients. The following chapter places this result in the broader context of convergence analysis and applications.
