Euler’s method is first-order accurate: halving the step size halves the global error. This chapter asks a natural question: can we do better without a fundamentally different strategy? The answer is yes — and the improvement costs exactly one additional evaluation of \(f\) per step. By sampling the slope at a carefully chosen second point inside the interval \([x_n, x_{n+1}]\), we can cancel the leading error term and achieve second-order accuracy. The surprising result is that there is not one second-order method but an entire family of them, parameterized by a single free constant. This chapter derives that family systematically and shows two of its most important members.
5.1 The Core Idea: Estimate the Slope at a Second Point
Euler’s method uses only the slope at the left endpoint of each interval:
\[y_{n+1} = y_n + h \cdot f(x_n, y_n)\]
The truncation error is \(O(h^2)\) per step because the method matches only the first two terms of the Taylor series: the constant term and the \(h \cdot y'(x_n)\) term. The \(h^2\) term and everything above it are discarded.
The key insight for improvement: instead of using only the slope at \(x_n\), what if we took a preliminary half-step to estimate where the solution will be at \(x_n + h/2\), then used the slope at that midpoint for the actual update? Geometrically, the slope at the midpoint of the interval is a better approximation to the average slope over \([x_n, x_{n+1}]\) than the slope at the left endpoint alone.
This intuition can be made precise. The two-variable Taylor series (established in Chapter 3) lets us expand \(f\) at any nearby point and compare the result term by term with the exact Taylor series for \(y(x_n + h)\). When the terms match through \(h^2\), the method is second-order accurate. The question is: what constraints must the coefficients satisfy to force that match?
5.2 The General Two-Stage RK Scheme
A two-stage Runge–Kutta method uses two slope evaluations per step. Following the unifying structure from Section 4.1, the update formula is:
The unknown coefficients are \(b_1\), \(b_2\), \(c_2\), and \(a_{21}\). There are four of them but, as we will see, only three independent equations. That underdetermination is precisely what produces the “family” of methods.
The subscript convention for the coefficients follows the Butcher tableau introduced in Section 5.10: \(b_i\) are the weights, \(c_i\) are the stage positions along the \(x\)-axis, and \(a_{ij}\) are the coupling coefficients that determine how each stage depends on earlier stages.
5.3 Expanding the Second Stage with the Two-Variable Taylor Series
To find what constraints the order conditions impose on our four coefficients, we need to expand \(k_2\) around \((x_n, y_n)\). This is exactly the two-variable Taylor expansion from Chapter 3, applied to \(f\) at the point \((x_n + c_2 h,\ y_n + a_{21} h k_1)\):
\[k_2 = f(x_n, y_n) + c_2 \cdot h \cdot f_x + a_{21} \cdot h \cdot k_1 \cdot f_y + O(h^2)\]
Since \(k_1 = f(x_n, y_n)\), the expansion simplifies to:
\[k_2 = f + c_2 \cdot h \cdot f_x + a_{21} \cdot h \cdot f \cdot f_y + O(h^2)\]
where \(f\), \(f_x\), and \(f_y\) are all evaluated at \((x_n, y_n)\) and the \(O(h^2)\) terms involve second-order partial derivatives.
Substituting this expansion into the update formula gives the numerical increment as a power series in \(h\). We will use SymPy to perform this expansion symbolically and then compare it term by term with the exact Taylor series for \(y(x_n + h)\).
First, let’s set up the symbolic framework. We treat \(f\) as an abstract function and let SymPy expand everything. Following the recipe from Chapter 3, we write the two-variable Taylor expansion of \(f\) at a displaced point as a direct sum of partial-derivative terms:
import sympy as spx, y, h, c2, a21, b1, b2 = sp.symbols('x y h c2 a21 b1 b2')f = sp.Function('f')def taylor2(dx, dy, order=2):"""Two-variable Taylor expansion of f at (x+dx, y+dy) around (x, y).""" f0 = f(x, y) expr = f0if order >=1: expr += dx * f0.diff(x) + dy * f0.diff(y)if order >=2: expr += (dx**2/2* f0.diff(x, 2)+ dx*dy * f0.diff(x).diff(y)+ dy**2/2* f0.diff(y, 2))return sp.expand(expr)# Manual collect-by-h. SymPy's sp.collect can fail on expressions with# multivariate Derivatives, so we group by h^k by extracting coefficients# directly.def collect_h(expr, sym, n_max=8): expr = sp.expand(expr)returnsum(expr.coeff(sym, k) * sym**k for k inrange(n_max +1))k1 = f(x, y)k2 = taylor2(c2*h, a21*h*k1, order=2)collect_h(k2, h)
The output shows \(k_2\) as \(f + h\,(c_2\, f_x + a_{21}\, f\, f_y) + O(h^2)\) (in SymPy’s Derivative notation for the partial derivatives). Now we form the full two-stage increment \(b_1 k_1 + b_2 k_2\) and compare it to the exact Taylor series for \(y(x + h)\).
5.4 Matching Coefficients: The Order Conditions for RK2
For the two-stage method to be second-order accurate, the numerical increment \(h\,(b_1 k_1 + b_2 k_2)\) must match the Taylor expansion of \(y(x_n + h)\) through terms of order \(h^2\). We compare the two expansions term by term.
The exact Taylor expansion of \(y(x_n + h)\) uses the chain rule. Since \(y'(x) = f(x, y(x))\), the second derivative is:
Four unknowns, three equations: the system is underdetermined by one degree of freedom. Choosing any value of \(b_2\) (other than 0) uniquely determines the remaining three:
This is the “RK2 family”. Every choice of \(b_2 \in (0, 1]\) produces a distinct but equally valid second-order Runge–Kutta method. SymPy confirms this:
The output is {b1: 1 - b2, c2: 1/(2*b2), a21: 1/(2*b2)}. The single free parameter \(b_2\) indexes the entire family. Two canonical choices are \(b_2 = 1\) (the Midpoint Method) and \(b_2 = 1/2\) (Heun’s Method). Both are genuinely different algorithms but both are second-order accurate.
5.6 The Midpoint Method: One Choice of Solution
Setting \(b_2 = 1\) gives \(b_1 = 0\), \(c_2 = 1/2\), \(a_{21} = 1/2\), placing the second stage exactly at the midpoint:
The method takes a half-step (\(h/2\)) to estimate the solution at the midpoint, evaluates \(f\) there, then uses that midpoint slope to advance the full step. The first slope \(k_1\) is used only to locate the midpoint — it does not appear in the final update. This is why \(b_1 = 0\).
5.7 Heun’s Method (Improved Euler): Another Choice
Setting \(b_2 = 1/2\) gives \(b_1 = 1/2\), \(c_2 = 1\), \(a_{21} = 1\). The second stage is now at \(x_n + h\) (the right endpoint of the interval), and the update uses the average of the slopes at the left and right endpoints:
\[k_1 = f(x_n,\ y_n)\]
\[k_2 = f(x_n + h,\ y_n + h \cdot k_1)\]
\[y_{n+1} = y_n + \frac{h}{2}\,(k_1 + k_2)\]
Heun’s method has a predictor–corrector interpretation: \(k_1\) is an Euler “predictor” step to estimate \(y_{n+1}\), then \(k_2\) re-evaluates \(f\) at that predicted point to “correct” the slope. The final update averages the two slopes. This interpretation is not just a mnemonic — it is the exact geometric content of setting \(c_2 = 1\) in the general scheme.
def heun_step(fx, x, y, h): k1 = fx(x, y) k2 = fx(x + h, y + h*k1)return y + (h/2)*(k1 + k2)
Three methods compared against the exact \(e^x\) (gray dashed).
import pandas as pdpd.DataFrame({"Method": ["Euler", "Midpoint", "Heun"],"Error at x=2": [ys_e[-1] - np.exp(2), ys_m[-1] - np.exp(2), ys_h[-1] - np.exp(2)]})
Method
Error at x=2
0
Euler
-1.428592
1
Midpoint
-0.126809
2
Heun
-0.126809
With \(h = 0.25\), Euler’s error should be roughly four times larger than Midpoint’s or Heun’s — consistent with Euler being \(O(h)\) and the others \(O(h^2)\).
A method is second-order if halving \(h\) reduces the global error by a factor of \(4\) (since \((h/2)^2 = h^2/4\)). We verify this by computing the error at \(x = 2\) for a sequence of step sizes and checking the ratio of successive errors.
def ratios(errs): out = ["---"]for i inrange(1, len(errs)): out.append(f"{errs[i]/errs[i-1]:.4f}")return outpd.DataFrame({"h": h_values,"Euler err": [f"{e:.4e}"for e in eur_errors],"ratio (E)": ratios(eur_errors),"Midpoint err": [f"{e:.4e}"for e in mid_errors],"ratio (M)": ratios(mid_errors),"Heun err": [f"{e:.4e}"for e in heun_errors],"ratio (H)": ratios(heun_errors),})
h
Euler err
ratio (E)
Midpoint err
ratio (M)
Heun err
ratio (H)
0
0.50000
2.3266e+00
---
4.1616e-01
---
4.1616e-01
---
1
0.25000
1.4286e+00
0.6140
1.2681e-01
0.3047
1.2681e-01
0.3047
2
0.12500
8.0581e-01
0.5641
3.4973e-02
0.2758
3.4973e-02
0.2758
3
0.06250
4.3039e-01
0.5341
9.1757e-03
0.2624
9.1757e-03
0.2624
4
0.03125
2.2278e-01
0.5176
2.3492e-03
0.2560
2.3492e-03
0.2560
The ratio columns confirm the order: Euler’s ratios approach \(0.5\) (first-order, since halving \(h\) halves the error), while Midpoint and Heun’s ratios approach \(0.25\) (second-order, since halving \(h\) quarters the error). Both RK2 methods achieve the same asymptotic convergence rate despite their different formulas.
5.10 The Butcher Tableau: A Compact Notation for RK Methods
The general two-stage RK scheme involves four coefficients: \(b_1\), \(b_2\), \(c_2\), and \(a_{21}\). As we move to three-stage (Chapter 6) and four-stage (Chapter 7) methods, the number of coefficients grows rapidly. Writing out every formula in full becomes unwieldy. The Butcher tableau is a compact notation that organizes all coefficients of any explicit RK method in a single structured array.
The tableau for a two-stage method looks like this:
The left column holds the \(c_i\) values (stage positions), the main body holds the \(a_{ij}\) coupling coefficients, and the bottom row holds the \(b_i\) weights. For explicit methods, the \(a\) matrix is strictly lower triangular: each stage only depends on stages that came before it.
The tableau encodes the entire algorithm. Given the tableau, you can reconstruct the formulas for \(k_1\), \(k_2\), and \(y_{n+1}\) mechanically. Given the formulas, you can fill in the tableau mechanically. The two representations carry identical information; the tableau is simply more compact.
This is pure notation — no new mathematics is introduced here. The same Taylor expansion logic from Sections 5.3 and 5.4 determines whether any tableau gives a valid second-order method. In Chapters 6 and 7 we will use the tableau to organize the growing list of coefficients as we move to three and four stages, and the pattern in the tableau will make the structure of higher-order methods visually apparent.
