DerivationODEsadvanced

Derivation: Heun's method

Full derivation of Heun: via the Taylor expansion to second order combined with the two-variable Taylor of f, and via the trapezoidal rule with an Euler prediction.

Route 1: Taylor to second order

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

    y(t+h)=y(t)+h y′(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 y′′y''. Differentiating y′(t)=f(t,y(t))y'(t)=f(t,y(t)) with the chain rule:

    y′′(t)=∂f∂t(t,y)+∂f∂y(t,y) y′(t)=ft+fy fy''(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+h ft+h f fy)+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+h ft+k fy+⋯f(t+h,\,y+k)=f+h\,f_t+k\,f_y+\cdots; choosing the increment k=h f(t,y)k=h\,f(t,y):

    f(t+h, y+hf(t,y))=f+h ft+h f fy+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 h⋅O(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 yk≈y(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+h f(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)