4  Higher Order Methods

Lecture 4

Published

February 2, 2026

Overview
  1. Error estimates for one-step methods
  2. Higher order Taylor methods
  3. Intro to Runge Kutta methods
Note

I encourage you to play around with the juptyer notebook for this lecture - you can copy the code with the this notebook button on the side of this page.

Warning

These notes are mainly a record of what we discussed and are not a substitute for attending the lectures and reading books! If anything is unclear/wrong, let me know and I will update the notes.

Tip

This lecture is mostly based on Burden, Faires, and Burden (2015) section 5.2, 5.3

Recall that we are considering IVPs: seek u :[0,T] \to \mathbb R such that

\begin{align} &u(0) = u_0 \nonumber\\ &u'(t) = f\big( t, u(t) \big). \tag{IVP} \end{align}

We will suppose that f is continuous on [0,T] \times \mathbb R and Lipschitz in the second argument and so the problem is well-posed and there exists a unique solution to (\text{IVP}). Recall Euler’s method: we approximate u on the equispaced mesh \{t_j = j h\}_{j=0}^n where h := \frac{T}{n} is the mesh size with the following difference equation

\begin{align} u_{j+1} = u_j + h f( t_j, u_j ). \tag{Euler} \end{align}

We saw last week that, if \|u''\|_{L^\infty} \leq M, then we obtain the bound | u(t_j) - u_j | \leq \frac{h M}{2L} \left( e^{Lt_j} - 1 \right). That is, the error decays with rate O(h) = O(n^{-1}) as n \to \infty.

Exercise 4.1 Prove that Euler’s method applied to (\text{IVP}) converges with rate O(h) in each of the following cases (these are the numerical experiments we saw last week):

  • (i) u_0 > 0, f(t, u) = \mu u, \mu > 0,
  • (ii) u_0 = 1, f(t,u) = u( 1 - \frac{u}{20}), T = 5,
  • (iii) u_0 = 0, f(t, u) = t^2 - u \sin t, T = 20,

Write down the corresponding error estimates. Compare your error bounds with the true errors that we observed last week. How pessimistic are the error bounds?

A natural question one may ask is “how can we do better than O(h) convergence?”:

4.1 General One-Step Methods

We now consider the following one-step method:

\begin{align} u_{j+1} = u_{j} + h \phi( t_j, u_j; h) \tag{1-step} \end{align}

for some choice of \phi.

function oneStep( u0, ϕ, T, n )
    h = T/n     
    t = 0:h:T    
    u = zeros(n+1)
    u[1] = u0
    for j = 1:n
        u[j+1] = u[j] + h * ϕ( t[j], u[j], h )
    end
    return u
end 
oneStep (generic function with 1 method)
Note

When \phi(t, u; h) := f(t, u), (\text{1-step}) is simply Euler’s method.

We first define the local truncation error for general one-step methods:

Definition 4.1 (Local truncation error) The local truncation error of the method (\text{1-step}) is given by

\begin{align} \tau_{j+1}(h) := \frac{u(t_{j+1}) - u(t_j)}{h} - \phi(t_j, u(t_j);h). \end{align}

We say (\text{1-step}) is consistent if

\begin{align} \max_j \tau_{j+1}(h) \to 0 \end{align}

as h\to 0.

It turns out that bounds on the local truncation error lead to a global error bound:

Theorem 4.1 Suppose that

  • (i) |\tau_{j+1}(h)| \leq Ch^p,
  • (ii) \left| \frac{\partial \phi}{\partial u} \right| \leq L

for all (t,u) \in[0,T] \times \mathbb R. Then, we have

\begin{align} | u(t_j) - u_j | \leq \frac{Ch^p}{L} \left( e^{L t_j} - 1 \right). \end{align}

Definition 4.2 (Order of accuracy) If \tau_{j+1}(h) = O(h^p) for some p, we say (\text{1-step}) has order of accuracy p.

Note

This is a more general version of the theorem we saw last week for Euler’s method. In that case we assumed that \|u''\|_{L^\infty} \leq M and showed |\tau_{j+1}(h)| \leq \frac{h M}{2}, giving order of accuracy 1.

Proof. The proof is exactly the same as in the case of Euler.

Define e_{j} := u(t_j) - u_j and note that

\begin{align} e_{j+1} &= u(t_{j+1}) - \left[ u_j + h \phi( t_j, u_j; h ) \right] \nonumber\\ % &= \left[ u(t_j) - u_j \right] \nonumber\\ % &\qquad + \left[ u(t_{j+1}) - u(t_j) - h \phi( t_j, u(t_j); h ) \right] \nonumber\\ % &\qquad + h \left[ \phi( t_j, u(t_j); h ) - \phi( t_j, u_j; h ) \right] \nonumber\\ % &= e_j + h \tau_{j+1}(h) + h \left[ \phi( t_j, u(t_j); h ) - \phi( t_j, u_j; h ) \right]. \end{align}

Using the fact that | \phi( t_j, u(t_j); h ) - \phi( t_j, u_j; h ) | \leq L | u(t_j) - u_j | = L|e_j| and |\tau_{j+_1}(h)| \leq Ch^p, we have

\begin{align} |e_{j+1}| \leq (1 + hL)|e_j| + C h^{p+1}. \end{align}

Applying A1 [Section B, Exercise 4 and 3], we conclude

\begin{align} |e_{j+1}| &\leq (1 + hL)^{j+1} \left( |e_0| + \frac{Ch^{p}}{L} \right) - \frac{Ch^{p}}{L} \nonumber\\ % &\leq \frac{Ch^{p}}{L} \left( (1 + hL)^{j+1} - 1 \right)\nonumber\\ % &\leq \frac{Ch^{p}}{L} \left( \big( e^{hL} \big)^{j+1} - 1 \right) % = \frac{Ch^{p}}{L} \left( e^{Lt_{j+1}} - 1 \right). \end{align}

Question: How to choose \phi to improve the order of accuracy?

4.2 Higher order Taylor methods

\begin{align} \frac{u(t_{j+1}) - u(t_j)}{h} % &= u'(t_j) + \frac{h}{2} u''(t_j) + \frac{h^2}{6} u'''(t_j) + \frac{h^3}{24} u^{(4)}(t_j) + \dots \end{align}

Burden, Faires, and Burden (2015) section 5.3.

Exercise 4.2 Write down the Taylor methods of order 1 and 2.

Example 4.1 Let’s solve

\begin{align} u'(t) &= u(t) - t^2 + 1 \nonumber\\ u(0) &= 0 \end{align}

numerically. You can check that this problem has a unique solution given by u(t) = (t+1)^2 - e^t. We use this example to analyse the error made in the Taylor method of order 1 (Euler’s method) and 2:

u0, T, n = 0., 5., 20
f(t,u) = u - t^2 + 1
f_t(t,u) = -2*t
f_u(t,u) = 1
ϕ(t, u, h) = f(t,u) + (h/2)*( f_t(t,u) + f_u(t,u) * f(t,u) )

u1 = oneStep( u0, (t,u,h)->f(t,u), T, n )
u2 = oneStep( u0, ϕ, T, n )

u_exact(t) = (t+1)^2 - exp(t)

plot(u_exact, 0, T, 
    label="exact",
    xlabel=L"t", ylabel=L"u(t)")

scatter!( 0:T/n:T, u1, 
    label="Euler (n=$n)")

scatter!( 0:T/n:T, u2, 
   label="Taylor order 2 (n=$n)" )

Here’s the error as a function of n:

N = 500
errs1 = zeros(N)
errs2 = zeros(N)
for n ∈ 1:N
    u1 = oneStep( u0, (t,u,h)->f(t,u), T, n )
    u2 = oneStep( u0, ϕ, T, n )
    mesh = 0:T/n:T
    errs1[n] = maximum( @. abs( u1 - u_exact(mesh) ) )
    errs2[n] = maximum( @. abs( u2 - u_exact(mesh) ) )
end

scatter( errs1, xaxis=:log, yaxis=:log, xlabel=L"n", label="Euler's method", lw=3, title="Error" )
plot!( 500*(1:N).^(-1), label=L"n^{-1}", linestyle=:dash )

scatter!( errs2, xaxis=:log, yaxis=:log, xlabel=L"n", label="Taylor order 2", lw=3 )
plot!( 1000*(1:N).^(-2), label=L"n^{-2}", linestyle=:dash )

However, the formulae for the higher order derivatives of f become more complicated and we don’t want to compute the derivatives of f explicitly.

4.3 Runge-Kutta Methods

Taylor theorem in multiple variables, Runge-Kutta of order 2: Burden, Faires, and Burden (2015) section 5.4

Midpoint method is often written as a multi-stage method:

\begin{align} k_1 &:= h f(t_j, u_j) \nonumber\\ k_2 &:= h f(t_j + \tfrac12 h, u_j + \tfrac12 k_1) \nonumber\\ u_{j+1} &:= u_j + k_2 \tag{Midpoint} \end{align}

Example 4.2 We return to Example 4.1 but now also implement the midpoint method

n = 20
ϕ2(t, u, h) = f(t+h/2,u+(h/2)*f(t,u))

u1 = oneStep( u0, (t,u,h)->f(t,u), T, n )
u2 = oneStep( u0, ϕ, T, n )
u3 = oneStep( u0, ϕ2, T, n )

plot(u_exact, 0, T, 
    label="exact",
    xlabel=L"t", ylabel=L"u(t)")

scatter!( 0:T/n:T, u1, 
    label="Euler (n=$n)")

scatter!( 0:T/n:T, u2, 
   label="Taylor order 2 (n=$n)" )

scatter!( 0:T/n:T, u3, 
   label="Midpoint (n=$n)" )

Here’s a plot of the errors:

N = 500
errs3 = zeros(N)
for n ∈ 1:N
    u3 = oneStep( u0, ϕ2, T, n )
    mesh = 0:T/n:T
    errs3[n] = maximum( @. abs( u3 - u_exact(mesh) ) )
end

scatter( errs1, xaxis=:log, yaxis=:log, xlabel=L"n", label="Euler's method", lw=3, title="Error" )
plot!( 500*(1:N).^(-1), label=L"n^{-1}", linestyle=:dash )

scatter!( errs2, xaxis=:log, yaxis=:log, xlabel=L"n", label="Taylor order 2", lw=3 )
scatter!( errs3, xaxis=:log, yaxis=:log, xlabel=L"n", label="Midpoint", lw=3 )
plot!( 1000*(1:N).^(-2), label=L"n^{-2}", linestyle=:dash )

Here, we can see the order of accuracy is indeed p=2.

Definition 4.3 A general s-stage Runge-Kutta method is of the form:

\begin{align} k_1 &= h f(t_j, u_j) \nonumber\\ k_2 &= h f(t_j + c_1 h, u_j + a_{11} k_1) \nonumber \\ k_3 &= h f(t_j + c_2 h, u_j + a_{21} k_1 + a_{22} k_2) \nonumber\\ \vdots& \nonumber\\ k_s &= h f(t_j + c_{s-1} h, u_j + a_{s-1,1} k_1 + \dots + a_{s-1,s-1} k_{s-1}) \nonumber\\ % u_{j+1} &= u_{j} + b_1 k_1 + b_2 k_2 + \dots + b_s k_s. \nonumber \end{align}

This Runge-Kutta method is conveniently written as the Butcher tableau:

\begin{array} {c|cccc} 0\\ c_1 & a_{11} \\ c_2 & a_{21} & a_{22} \\ \vdots & \vdots& \vdots & \ddots& \\ c_{s-1} & a_{s-1,1}& a_{s-1,2} & \cdots & a_{s-1,s-1} \\ \hline & b_1 & b_2 & \dots & b_{s-1} & b_s \end{array}

Exercise 4.3 Write down the Butcher tableau for (i) Euler’s method, (ii) Midpoint method, and the following 2-stage, order 2 methods:

(iii): \begin{align} u_{j+1} &= u_j + \frac{h}{2} \left[ f(t_j, u_j) + f\big( t_{j+1}, u_j + h f(t_j, u_j) \big) \right] \end{align}

(iv): \begin{align} u_{j+1} &= u_j + \frac{h}{4} \left[ f(t_j, u_j) + 3 f\big( t_{j} + \tfrac{2}{3} h, u_j + \tfrac{2}{3} h f(t_j, u_j) \big) \right] \end{align}

Exercise 4.4 Show that the Modified Euler method has order of accuracy 2.

TipTO-DO
  • Assignment 2: Due next Wednesday,
  • Read: Burden, Faires, and Burden (2015) sections 5.1 - 5.4 - please come to the problem class on Wednesday with any questions you have,
  • (Sign-up to office hours if you would like)

Next: Problem class

Burden, Richard L., Douglas J. Faires, and Annette M. Burden. 2015. Numerical Analysis. 10th ed. CENGAGE Learning.