Lecture 7: Zero-stability of Runge-Kutta and multistep methods
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.
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.
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 one-step methods:
\begin{align} u_{j+1} = u_j + h \phi(t_j, u_j; h) \tag{1-step} \end{align}
where t_j = j \frac{T}{n} = j h defines the mesh. We saw last week that, if the local truncation error is \tau_{j+1} = O(h^p) as h\to 0 (i.e. if the method is consistent with order of accuracy p) and |\frac{\partial \phi}{\partial u}| \leq L on [0,T] \times \mathbb R (i.e. the method is zero-stable), then we have the error bound
\begin{align} |u(t_j) - u_j| \leq \frac{Ch^p}{L} ( e^{Lt_j} - 1) \end{align}
(i.e. the method is convergent). This result applies to Runge-Kutta methods:
function RK( u0, f, T, n, A, b, c )
h, t = T/n, 0:T/n:T
u = fill(float(u0), n+1)
s, = size(A)
s = s+1
k = fill(float(u0), s)
for j = 1:n
k[1] = h * f( t[j], u[j] )
for i = 1:s-1
k[i+1] = h * f( t[j] + h * c[i], u[j] + sum( A[i,ℓ] * k[ℓ] for ℓ ∈ 1:i ) )
end
u[j+1] = u[j] + sum( b[ℓ] * k[ℓ] for ℓ ∈ 1:s )
end
return u
end RK (generic function with 1 method)Burden, Faires, and Burden (2015) section 5.6
All the methods we have seen so take the form u_{j+1} = u_j + h \phi( t_j, u_j; h). These are known as 1-step methods because u_{j+1} is only given as a function of 1 previously computed value u_j. Multistep methods hope to improve efficiency by using the history of the solution \{u_j, u_{j-1}, \dots, u_{j-m+1}\} for some m \in \mathbb N to compute the next term u_{j+1}. Linear multistep methods produce u_{j+1} as a linear combination of \{ u_{i} \} and \{ f(t_{i}, u_{i}) \}:
Definition 7.1 (Linear multistep method) A linear m-step method is of the form:
\begin{align} u_{j+1} &= \alpha_{m-1} u_{j} + \alpha_{m-2} u_{j-1} + \cdots + \alpha_{0} u_{j-m+1} \nonumber\\ % &+ h \big[ \beta_m f(t_{j+1}, u_{j+1}) + \beta_{m-1} f(t_{j}, u_{j}) + \cdots + \beta_0 f(t_{j-m+1}, u_{j-m+1})\big] \end{align}
for some coefficients \alpha_0,\dots,\alpha_{m-1} and \beta_0, \dots, \beta_{m}.
If \beta_m = 0, the method is explicit (u_{j+1} is given explicitly),
If \beta_m \not= 0, the method is implicit (u_{j+1} solves a non-linear equation that must be solved).
For 1-step methods, we started from u_0 given by the equation.
Now, we need to initiate u_0,\dots,u_{m-1} in order to define u_{m}. This can be done by e.g. using a RK method.
We can store the history f(t_{j-1}, u_{j-1}), \dots, f(t_{j-m+1}, u_{j-m+1}) so that (in the explicit case) to evaluate u_{j+1} we only need to compute one function evaluation f(t_{j}, u_{j}).
Here is a basic implementation which uses RK4 to compute u_1,\dots,u_{m-1}. In the implicit case, we need to solve a nonlinear problem in 1d (See Chapter 2 of Math 5485 for the details). To do this, we use Roots (Exercise: what method does this package use?),
using Rootsfunction linearMultistep( u0, f, T, n, α, β)
m, h, t = length(α), T/n, 0:T/n:T
u = fill(float(u0), n+1)
# initiate u_1,...,u_m-1 using RK4
if (m > 1)
u[1:m] = RK( u0, f, (m-1)*h, m-1,
[1/2 0 0 ; 0 1/2 0 ; 0 0 1],
[1/6 1/3 1/3 1/6], [1/2 1/2 1] )
end
f_hist = [ f(t[m-ℓ],u[m-ℓ]) for ℓ ∈ 1:m-1 ]
for j = m:n
f_hist = [ f(t[j],u[j]), f_hist[1:m-1]... ]
if (length(β) == m)
u[j+1] = sum( α[m-ℓ] * u[j-ℓ] + h * β[m-ℓ] * f_hist[ℓ+1] for ℓ ∈ 0:(m-1) )
else
r = x -> x - h * β[m+1] * f(t[j+1], x) - sum( α[m-ℓ] * u[j-ℓ] + h * β[m-ℓ] * f_hist[ℓ+1] for ℓ ∈ 0:(m-1) )
u[j+1] = find_zero( r, u[j] )
end
end
return u
end linearMultistep (generic function with 1 method)Example 7.1 (Adams-Bashforth (AB)) AB methods are explicit (\beta_{m} = 0) with \alpha_{m-1} = 1 and \alpha_{m-2} = \cdots = 0. They are chosen to have order p=m and denoted \text{(AB}p\text{)}:
\begin{align} \tag{AB1} &m=1, \quad \beta_0 = 1 \\ \tag{AB2} &m=2, \quad \beta_1 = \frac32, \quad \beta_0 = -\frac12 \\ \tag{AB3} &m=3, \quad \beta_2 = \frac{23}{12}, \quad \beta_1 = -\frac{16}{12}, \quad \beta_0 = \frac5{12} \end{align}
Exercise 7.1 Notice that \text{(AB1)} is Euler’s method.
We define the local truncation error in an analogous way to 1-step methods:
Definition 7.2 The local truncation error of a linear multistep method with coefficients \alpha_{m-1}, \dots, \alpha_0 and \beta_{m},\dots,\beta_0 is given by
\begin{align} \tau_{j+1}(h) &= \frac1h \Big[ u(t_{j+1}) - \sum_{i=0}^{m-1} \alpha_{m-1-i} u(t_{j-i}) \Big] \nonumber\\ % &\quad - \sum_{i=0}^m \beta_{m-i} f\big( t_{j+1-i}, u(t_{j+1-i}) \big) \end{align}
We say the multistep method is consistent if \max_j \tau_{j+1}(h) \to 0 as h\to 0.
We say the multistep method has order of accuracy p if \tau_{j+1}(h) = O(h^p) as h\to 0.
Exercise 7.2 Show that \text{(AB2)} has order of accuracy 2.
Example 7.2 (Adams-Moulton (AM)) AM methods are implicit (\beta_{m} \not= 0) with \alpha_{m-1} = 1 and \alpha_{m-2} = \cdots = 0. Since there is another degree of freedom (\beta_m\not=0), m-step implicit methods can have order of accuracy p = m+1:
\begin{align} \tag{AM1} &m=1, \quad \beta_1 = 1, \quad \beta_0 = 0 \\ \tag{AM2} &m=1, \quad \beta_1 = \frac12, \quad \beta_0 = \frac12 \\ \tag{AM3} &m=2, \quad \beta_2 = \frac{5}{12}, \quad \beta_1 = \frac{8}{12}, \quad \beta_0 =- \frac1{12} \end{align}
Here, we use the naming convention (\text{AM}p) where p is the order of accuracy.
Exercise 7.3 Write down the formulas for (\text{AM}1), (\text{AM}2), (\text{AM}3).
(\text{AM}1) is also called backward’s Euler because you can derive it from u'(t_{j+1}) = f(t_{j+1}, u_{j+1}) and replace the derivative with the backwards difference \frac{u_{j+1} - u_j}{h}.
(\text{AM}2) is also called the trapezoid rule or Crank-Nicolson
Example 7.3 Using our function linearMultistep(), we can try and numerically solve
\begin{align} u(0) &= 0 \nonumber\\ u'(t) &= \cos u(t) + \sin t \end{align}
on [0, 5]. Again, we plot the solution using a built-in function with a small tolerance so we have something to compare with.
u0, T, n = 0., 10., 20
f(t,u) = (u == Inf) ? Inf : cos(u) + sin(t)
ivp = ODEProblem((u,p,t)->f(t,u), u0, (0., T))
sol = solve(ivp, Tsit5(), reltol=1e-15, abstol=1e-15)
u1 = linearMultistep(u0, f, T, n, [0 1], [-1/2 3/2]) # AB2
u2 = linearMultistep(u0, f, T, n, [1], [1/2 1/2]) # AM2
plot(sol,
label="Tsit5()",
xlabel=L"t", ylabel=L"u(t)")
plot!(0:T/n:T, u1, m=3, label="AB2")
plot!(0:T/n:T, u2, m=3, label="AM2")
We plot the errors of various methods that we have seen:
nn = [round(Int, 10^ℓ) for ℓ ∈ 1:0.5:4]
errs = zeros(length(nn),8)
for (ℓ,n) ∈ enumerate(nn)
u = zeros(n+1, 7)
# order 1
u[:,1] = RK( u0, f, T, n, [], [1], []) # Euler
u[:,2] = linearMultistep(u0, f, T, n, [1], [0 1]) # Backwards Euler
# order 2
u[:,3] = linearMultistep(u0, f, T, n, [0 1], [-1/2 3/2]) # AB2
u[:,4] = linearMultistep(u0, f, T, n, [1], [1/2 1/2]) # AM2 = trapezoid
# order 4
u[:,5] = linearMultistep(u0, f, T, n, [0 0 0 1], [-9/24 37/24 -59/24 55/24]) # AB4
u[:,6] = linearMultistep(u0, f, T, n, [0 0 1], [1/24 -5/24 19/24 9/24]) # AM4
u[:,7] = RK( u0, f, T, n,
[1/2 0 0 ; 0 1/2 0 ; 0 0 1],
[1/6 1/3 1/3 1/6], [1/2 1/2 1] ) # RK4
mesh = 0:T/n:T
for i ∈ 1:7
errs[ℓ,i] = maximum( @. abs( u[:,i] - sol(mesh) ) )
end
end
plot(nn, errs[:,1:7];
m=3, label=["AB1 = Euler" "AM1 = Backwards Euler" "AB2" "AM2 = Trapezoid" "AB4" "AM4" "RK4"],
xaxis=(:log10, L"n"), yaxis=:log10,
title="Error", legend=:bottomleft)
plot!( nn, (1/2)*errs[end,1]*(nn/nn[end]).^(-1),
label=L"O(n^{-1})", linestyle=:dash, primary=false )
plot!( nn, (1/2)*errs[end,4]*(nn/nn[end]).^(-2),
label=L"O(n^{-2})", linestyle=:dash, primary=false )
plot!( nn, (1/2)*errs[end,7]*(nn/nn[end]).^(-4),
label=L"O(n^{-4})", linestyle=:dash, primary=false)
Here, we can see that all the methods introduced above are convergent (in this particular case). But why?
We saw that, if |\frac{\partial \phi}{\partial u}\big| \leq L, then the 1-step methods (in particular Runge-Kutta methods) are convergent and their rate of convergence is given by the order of accuracy (i.e. the rate of convergence is given by the local truncation error). This is not necessarily true for general linear multistep methods:
Example 7.4 Here, we return to the equation from Example 10.1 but now use the method
\begin{align} u_{j+1} &= (1-\alpha) u_j + \alpha u_{j-1} \nonumber\\ &\,\, + \frac{h}{2} \big( (\alpha + 3) f_j + (\alpha-1) f_{j-1} ) % \tag{$\alpha2$} \end{align}
where f_i := f( t_i, u_i ).
Note: calling this method (\alpha2) is non-standard notation!
Exercise 7.4 Show that this method has order of accuracy p=2. Extra: for what \alpha is this method 3rd order accurate?
What do you expect to happen if we apply this method with different choices of \alpha?
n1, n2 = 20, 100
α = -1.001
u1 = linearMultistep(u0, f, T, n1, [α, 1-α], [(α-1)/2, (α+3)/2])
u2 = linearMultistep(u0, f, T, n2, [α, 1-α], [(α-1)/2, (α+3)/2])
P = plot(sol,
label="Tsit5()",
xlabel=L"t", ylabel=L"u(t)")
plot!(0:T/n1:T, u1, m=3, label="α = $α, n = $n1")
plot!(0:T/n2:T, u2, m=3, label="α = $α, n = $n2", ylims=(0,3))
nn = [round(Int, 10^ℓ) for ℓ ∈ 1:0.5:5]
errs = zeros(length(nn))
for (ℓ,n) ∈ enumerate(nn)
u = linearMultistep(u0, f, T, n, [α, 1-α], [(α-1)/2, (α+3)/2])
mesh = 0:T/n:T
errs[ℓ] = maximum( @. abs( u - sol(mesh) ) )
end
P′ = plot(nn, errs;
m=3, xaxis=(:log10, L"n"), yaxis=(:log10, "Error"),
legend=false)
plot( P, P′, layout=(1,2), size=(1000,500))
Aim of today: understand what is happening here!
In order to understand what is happening here, we need to define the following characteristic polynomials:
Definition 7.3 We define the characteristic polynomials
\begin{align} \rho(z) &= z^m - \alpha_{m-1} z^{m-1} - \alpha_{m-2} z^{m-2} - \cdots - \alpha_0 \nonumber\\ \sigma(z) &= \beta_m z^m + \beta_{m-1} z^{m-1} + \cdots + \beta_0. \end{align}
These are also called generating polynomials.
Exercise 7.5 Write down the characteristic polynomials associated to the method (\alpha2). What are the roots of \rho in this case.
We now introduce the following (very simple) IVP:
\begin{align} u(0) &= u_0 \not= 0 \nonumber\\ u'(t) &= 0 \tag{test}. \end{align}
Notice that u(t) = u_0 is the exact solution to this equation. When approximating (\text{test}) with our linear multistep method, we obtain the difference equation
\begin{align} u_{j+1} = \alpha_{m-1} u_j + \alpha_{m-2} u_{j-1} + \dots + \alpha_0 u_{j+m-1}. \end{align}
Suppose that \zeta_0 is a root of \rho. Then, u_{j} := (\zeta_0)^j solves (6). In fact, if all roots, \zeta_0,\dots,\zeta_{m-1}, of \rho are simple, then the general solution to (6) is given by
\begin{align} u_j := c_0 (\zeta_0)^j + \cdots + c_{m-1} (\zeta_{m-1})^j \end{align}
for some constants c_0, \cdots, c_{m-1}. By consistency of the multistep method, we must have \alpha_0 + \dots + \alpha_m = 1 and so u_j := u_0 must be a solution of (6). As a result, the general solution to (6) is of the form
\begin{align} u_j := u_0 + c_1 (\zeta_1)^j + \cdots + c_{m-1} (\zeta_{m-1})^j \end{align}
where 1,\zeta_1,\dots,\zeta_{m-1} are the roots of \rho and c_1,\dots,c_{m-1} are some constants. The constants c_1,\dots,c_{m-1} are in general non-zero due to round-off errors. This error grows exponentially unless |\zeta_j| \leq 1 for each of the roots \zeta_j of \rho.
Definition 7.4 We say a linear multistep method satisfies the root condition if all roots \zeta of the characteristic polynomial \rho satisfy |\zeta|\leq 1 and, if |\zeta| = 1, then \zeta is a simple root.
Exercise 7.6 For what \alpha does (\alpha2) satisfy the root condition?
Theorem 7.1 A linear multistep method is zero-stable if and only if it satisfies the root condition.
A consistent linear multistep method is convergent if and only if it is zero-stable.
Next: Chapter 2: Iterative Methods in Matrix Algebra