Heun's method

Averaging the initial slope and a predicted slope yields a second-order method. Full derivation via second-order Taylor and via the trapezoidal rule with an Euler prediction.

Averaging two slopes

Euler uses only the slope at the start of the subinterval, which is why it drifts as soon as the solution curves. Heun's method corrects that bias by averaging the initial slope with the slope at the point Euler would reach:

yk+1=yk+12k1+12k2,k1=hf(tk,yk)k2=hf(tk+1,yk+k1)y_{k+1}=y_k+\frac{1}{2}k_1+\frac{1}{2}k_2,\qquad \begin{aligned} k_1&=h\,f(t_k,y_k)\\ k_2&=h\,f(t_{k+1},\,y_k+k_1) \end{aligned}
Heun's method, order 2.
tₖtₖ₊₁ỹₖ₊₁predicted pointaverage slope
Heun first predicts with Euler and then corrects using the average of the initial slope and the predicted final slope.
Expand diagram

Heun first predicts with Euler and then corrects using the average of the initial slope and the predicted final slope.

Derivation

There are two natural routes to the formula: force the scheme to reproduce the Taylor expansion up to second order, or approximate the integral form of the IVP with the trapezoidal rule and predict the unknown value with Euler.

DerivationDerivation: Heun's methodView as its own page →

Route 1: Taylor to second order

  1. Expand the solution by Taylor one order further than in Euler:

    y(t+h)=y(t)+hy(t)+h22y(t)+O(h3)y(t+h)=y(t)+h\,y'(t)+\frac{h^2}{2}y''(t)+\mathcal{O}(h^3)
  2. We need yy''. Differentiating y(t)=f(t,y(t))y'(t)=f(t,y(t)) with the chain rule:

    y(t)=ft(t,y)+fy(t,y)y(t)=ft+fyfy''(t)=\frac{\partial f}{\partial t}(t,y)+\frac{\partial f}{\partial y}(t,y)\,y'(t)=f_t+f_y\,f
  3. Substituting y=fy'=f and y=ft+fyfy''=f_t+f_y f into the expansion and grouping half and half:

    y(t+h)=y(t)+h2f+h2(f+hft+hffy)+O(h3)y(t+h)=y(t)+\frac{h}{2}f+\frac{h}{2}\Bigl(f+h\,f_t+h\,f\,f_y\Bigr)+\mathcal{O}(h^3)
  4. The bracket matches a two-variable Taylor expansion of ff. In general, f(t+h,y+k)=f+hft+kfy+f(t+h,\,y+k)=f+h\,f_t+k\,f_y+\cdots; choosing the increment k=hf(t,y)k=h\,f(t,y):

    f(t+h,y+hf(t,y))=f+hft+hffy+O(h2)f\bigl(t+h,\,y+h f(t,y)\bigr)=f+h\,f_t+h\,f\,f_y+\mathcal{O}(h^2)
  5. Replacing the bracket by that evaluation (the hO(h2)h\cdot\mathcal{O}(h^2) error is absorbed into O(h3)\mathcal{O}(h^3)):

    y(t+h)=y(t)+h2f(t,y)+h2f(t+h,y+hf(t,y))+O(h3)y(t+h)=y(t)+\frac{h}{2}f(t,y)+\frac{h}{2}f\bigl(t+h,\,y+h f(t,y)\bigr)+\mathcal{O}(h^3)
  6. Evaluating at t=tkt=t_k with yky(tk)y_k\approx y(t_k) and discarding O(h3)\mathcal{O}(h^3) yields Heun's method, with local error O(h3)\mathcal{O}(h^3) and hence global error O(h2)\mathcal{O}(h^2):

    yk+1=yk+12hf(tk,yk)k1+12hf(tk+1,yk+k1)k2y_{k+1}=y_k+\frac{1}{2}\underbrace{h f(t_k,y_k)}_{k_1}+\frac{1}{2}\underbrace{h f(t_{k+1},\,y_k+k_1)}_{k_2}

Route 2: trapezoid with Euler prediction

tₖtₖ₊₁ỹₖ₊₁predicted pointaverage slope
Heun interprets the step as a trapezoid in the integral of slopes: one slope is known and the other is obtained with an Euler prediction.
Expand diagram

Heun interprets the step as a trapezoid in the integral of slopes: one slope is known and the other is obtained with an Euler prediction.

  1. Start from the exact integral form of the IVP on one subinterval:

    y(tk+1)=y(tk)+tktk+1f(τ,y(τ))dτy(t_{k+1})=y(t_k)+\int_{t_k}^{t_{k+1}} f\bigl(\tau,y(\tau)\bigr)\,d\tau
  2. Approximate that integral with the trapezoidal rule:

    y(tk+1)yk+h2(f(tk,yk)+f(tk+1,y(tk+1)))y(t_{k+1})\approx y_k+\frac{h}{2}\Bigl(f(t_k,y_k)+f\bigl(t_{k+1},y(t_{k+1})\bigr)\Bigr)
  3. The value y(tk+1)y(t_{k+1}) on the right-hand side is unknown. Instead of solving the implicit equation (that would be AM2, the implicit trapezoid), predict it with one Euler step:

    yˉk+1=yk+hf(tk,yk)\bar y_{k+1}=y_k+h\,f(t_k,y_k)
  4. Substituting the prediction into the trapezoid gives the same formula as via Taylor: Heun is a one-step predictor-corrector pair.

    yk+1=yk+h2(f(tk,yk)+f(tk+1,yˉk+1))y_{k+1}=y_k+\frac{h}{2}\Bigl(f(t_k,y_k)+f(t_{k+1},\bar y_{k+1})\Bigr)

Order and relation to other methods

Seen from the second derivation, Heun is the simplest predictor-corrector pair: it predicts with Euler and corrects with the trapezoidal rule. It is the same idea that, with multistep methods, produces the Adams-Bashforth-Moulton pairs. It is also the best-known member of the two-stage Runge-Kutta family, and the natural step towards the classical fourth-order Runge-Kutta.

ExampleOne Heun step

Take one Heun step with h=1h=1 for y=(30.1y)yy'=(3-0.1y)y, y(0)=10y(0)=10.

  1. Initial slope: k1=hf(0,10)=1(31)10=20k_1=h\,f(0,10)=1\cdot(3-1)\cdot 10=20.

  2. Predicted slope at Euler's arrival point, y0+k1=30y_0+k_1=30: k2=hf(1,30)=(33)30=0k_2=h\,f(1,30)=(3-3)\cdot 30=0.

  3. Average of both slopes:

    y1=10+1220+120=20y_1=10+\tfrac{1}{2}\cdot 20+\tfrac{1}{2}\cdot 0=20

The exact value is y(1)=27.2833y(1)=27.2833: with a single step the error is still large, but the table in the order exercise shows that as the mesh is refined Heun's error falls like h2h^2, much faster than Euler's.