Lecture 5
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.
This lecture is mostly based on Burden, Faires, and Burden (2015) section 5.4 & 5.5
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 define 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, then we have the error bound
\begin{align} |u(t_j) - u_j| \leq \frac{Ch^p}{L} ( e^{Lt_j} - 1). \end{align}
Then, we looked at particular examples of \phi:
Taylor methods: use the expansion u(t_{j+1}) = u(t_j) + h u'(t_j) + \frac{h^2}{2} u''(t_j) + \dots and use the fact that u'(t) = f(t, u(t)) to evaluate the higher derivatives:
\begin{align} \frac{u(t_{j+1}) - u(t_j)}{h} % &= f(t_j, u(t_j)) \nonumber\\ % &+ \frac{h}{2} \left[ \frac{\partial f}{\partial t}(t_j, u(t_j)) + \frac{\partial f}{\partial u}(t_j, u(t_j)) f(t_j, u(t_j))\right] + \dots \end{align}
\begin{align} \phi = f + \frac{h}{2}\Big( \frac{\partial f}{\partial t} + \frac{\partial f}{\partial u} f \Big). \end{align}
The downside of these methods is that they require higher order derivatives of f to be known exactly.
Last week, we looked for an order 2 method with \phi = f(t + \alpha, u + \beta). That is, we expanded
\begin{align} \kappa_{j+1}(h) := \frac{u(t_{j+1}) - u(t_j)}{h} % - f\big( t_j + \alpha, u(t_j) + \beta \big) \nonumber \end{align}
in \alpha, \beta and found that, under the choice \alpha := \frac{h}{2} and \beta := \frac{h}{2} f(t_j, u_j), the local truncation error is \kappa_{j+1}(h) = O(h^2). This method is known as the midpoint method:
\begin{align} u_{j+1} = u_j + h f\big( t_j + \tfrac{h}{2}, u_j + \tfrac{h}{2} f(t_j,u_j) \big). \tag{midpoint} \end{align}
We can also write this 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}
Definition 5.1 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 5.1 Write down the Butcher tableau for (i) Euler’s method, (ii) Midpoint method (also known as improved Euler’s method), and the following 2-stage, order 2 methods:
(iii): Modified Euler \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): Heun’s method \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}
Example 5.1 Show that all 2 stage, order 2 Runge-Kutta methods have a Butcher tableau of the form
\begin{array} {c|cccc} 0\\ \tfrac{1}{2b} & \tfrac{1}{2b} \\ \hline & 1-b & b \end{array}for some b \not= 0.
Here, we implement a general RK method:
function RK( u0, f, T, n, A, b, c )
h = T/n
t = 0:h: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)Example 5.2 Here, we solve
\begin{align} u(0) &= u_0 \nonumber\\ u'(t) &= \cos u(t) + \sin t \end{align}
with different initial conditions using the OrdinaryDiffEq package:
f(t,u) = cos(u) + sin(t)
u0, T = 0., 10.
ivp = ODEProblem((u,p,t)->f(t,u), u0, (0., T))
sol = solve(ivp, Tsit5());
P = plot(sol,
title=L"u'=f(t,u)", legend=false, label="u0 = $u0",
xlabel=L"t", ylabel=L"u(t)")
for u0 = vcat(0:.1:2π, [π/2,π,3π/2,2π])
ivp = ODEProblem((u,p,t)->f(t,u), u0, (0.,T))
sol = solve(ivp, Tsit5());
plot!(sol,
label="u0 = $u0",
xlabel=L"t", ylabel=L"u(t)")
end
display(P)
We next apply (i) Euler, (ii) midpoint, (iii) modified Euler, and (iv) Heun’s method and plot the errors as a function of n:
u0 = 1.
ivp = ODEProblem((u,p,t)->f(t,u), u0, (0.,T))
sol = solve(ivp, Tsit5(), reltol=1e-15, abstol=1e-15);
nn = [round(Int, 10^ℓ) for ℓ ∈ 0:0.5:4]
errs = zeros(length(nn),5)
i = 0
for (ℓ,n) ∈ enumerate(nn)
u1 = RK( u0, f, T, n, [], [1], [])
u2 = RK( u0, f, T, n, [1/2], [0, 1], [1/2])
u3 = RK( u0, f, T, n, [1], [1/2, 1/2], [1])
u4 = RK( u0, f, T, n, [2/3], [1/4, 3/4], [2/3])
mesh = 0:T/n:T
errs[ℓ,1] = maximum( @. abs( u1 - sol(mesh) ) )
errs[ℓ,2] = maximum( @. abs( u2 - sol(mesh) ) )
errs[ℓ,3] = maximum( @. abs( u3 - sol(mesh) ) )
errs[ℓ,4] = maximum( @. abs( u4 - sol(mesh) ) )
end
plot(nn, errs[:,1:4];
m=3, label=["Euler" "Midpoint" "Modified Euler" "Heun"],
xaxis=(:log10, L"n"), yaxis=:log10,
title="Error")
plot!( nn, (1/2)*errs[end,1]*(nn/nn[end]).^(-1), label=L"O(n^{-1})", linestyle=:dash )
plot!( nn, (1/2)*errs[end,4]*(nn/nn[end]).^(-2), label=L"O(n^{-2})", linestyle=:dash )
Section 5.4 of Burden, Faires, and Burden (2015).
The most common RK method is 4-stage, order 4 method given by the followng Butcher tableau:
\begin{align}\tag{RK4} \begin{array} {c|cccc} 0\\ \frac12 & \frac12 \\ \frac12 & 0 & \frac12 \\ 1 & 0 & 0 & 1 \\ \hline & \frac16 & \frac13 & \frac13 & \frac16 \end{array} \end{align}
Example 5.3 We apply (\text{RK4}) to Example 5.2:
A = [1/2 0 0 ; 0 1/2 0 ; 0 0 1]
b = [1/6 1/3 1/3 1/6]
c = [1/2 1/2 1]
for (ℓ,n) ∈ enumerate(nn)
u5 = RK( u0, f, T, n, A, b, c)
mesh = 0:T/n:T
errs[ℓ,5] = maximum( @. abs( u5 - sol(mesh) ) )
end
plot(nn, errs;
m=4, label=["Euler" "Midpoint" "Modified Euler" "Heun" "RK4"],
xaxis=(:log10, L"n"), yaxis=:log10, ylims=(1e-11,1e2),
title="Error")
plot!( nn, (1/2)*errs[end,1]*(nn/nn[end]).^(-1), label=L"O(n^{-1})", linestyle=:dash )
plot!( nn, (1/2)*errs[end,4]*(nn/nn[end]).^(-2), label=L"O(n^{-2})", linestyle=:dash )
plot!( nn, (1/100)*errs[end,5]*(nn/nn[end]).^(-4), label=L"O(n^{-4})", linestyle=:dash )
We want to compare the error made with the computational cost involved (rather than n).
Euler:
\begin{align} u_{j+1} = u_j + h f( t_j, u_j) \end{align}
General 2 stage RK:
\begin{align} u_{j+1} = u_j &+ b_1 h f(t_j, u_j) \nonumber\\ &+ b_2 h f\big( t_j + c_1 h, u_j + a_{11} h f(t_j, u_j) \big) \end{align}
Exercise 5.2 How many function evaluations are needed to implement these methods? What about a general s-stage RK method?
Example 5.4 We return to Example 5.2 and plot the errors against the function evaluations:
A general s-stage RK method requires s function evaluations.
(For 5 \leq s \leq 7 the maximal order of accuracy is s-1 and, for 8 \leq s \leq 9 the maximal order of accuracy is s-2 and, for s \geq 10, the maximal order of accuracy is s-3.)
Here, we look at example 6.1.8 from Driscoll and Braun (2018):
Example 5.5 A model for the velocity v(t) of a skydiver at time t is given by
\begin{align} v(0) &= 0 \nonumber\\ \dot{v}(t) &= -g + \frac{k(t)}{m} v(t)^2 \end{align}
where g \approx 9.81 \text{m/s}^2 is the acceleration due to gravity, m is the mass of the skydiver, and k(t) models the effects of air-resistance. A simple model for k is given as a step function:
\begin{align} k(t) &= \begin{cases} k_1 & \text{if }t \leq t_1 \nonumber\\ k_2 & \text{if }t > t_1 \end{cases} \end{align}
where k_2 > k_1 and t_1 is the time at which the parachute is deployed.
What do you expect the solution to look like?
Here, we implement the RK4 method and the scheme included in OrdinaryDiffEq:
t1 = 2π
k(t) = t<=t1 ? 1 : 50
f(t,v) = -9.81 + (1/80) * k(t) * v^2
v0, T, n = 0., 10., 100
ivp = ODEProblem((v,p,t)->f(t,v), v0, (0., T))
sol = solve(ivp, Tsit5())
N = length(sol.t)
v = RK( v0, f, T, n, A, b, c)
plot(sol,
title="Velocity of a skydiver",
label="Tsit5 (n=$N)", xlabel=L"t", ylabel=L"v(t)",
ylims=(-28,0))
plot!(0:T/n:T, v, label="RK4 (n=$n)")
plot(sol.t, sol.u, m=5,
title="Velocity of a skydiver",
label="Tsit5 (n=$N)", xlabel=L"t", ylabel=L"v(t)",
xlims=(2π-.5,2π+.5), ylims=(-28,-1))
plot!(0:T/n:T, v, m=3, label="RK4 (n=$n)")
Exercise 5.3 Complete exercise 6.1.8 of Driscoll and Braun (2018).
Section 5.5 of Burden, Faires, and Burden (2015)
Here, we implement a adaptive Runge-Kutta algorithm:
function ARK( u0, f, T; A=[1/2 0 0 ; 0 3/4 0 ; 2/9 1/3 4/9], b=[2/9 1/3 4/9 0; 7/24 1/4 1/3 1/8], c=[1/2 3/4 1], tol=1e-5, p=2)
t = [0.]
u = [float(u0)]
s, = size(A)
s = s+1
k = fill(float(u0), s)
i = 1
h = 1
k = fill(float(u0), s)
# Time stepping.
while t[i] < T
if h < 1e-15
@warn "Step size is too small, t=$(t[i])"
break
end
k[1] = h * f( t[i], u[i] )
for j = 1:s-1
k[j+1] = h * f( t[i] + h * c[j], u[i] + sum( A[j,ℓ] * k[ℓ] for ℓ ∈ 1:j ) )
end
next_u_1 = u[i] + sum( b[1,ℓ] * k[ℓ] for ℓ ∈ 1:(s-1) )
next_u_2 = u[i] + sum( b[2,ℓ] * k[ℓ] for ℓ ∈ 1:s )
E = abs(next_u_1 - next_u_2)
if E < tol
push!(t, t[i] + h)
push!(u, next_u_2)
i += 1
end
# step size
q = 0.85 * (tol/E)^(1/(p+1))
q = min(q, 2)
h = min(q*h, T-t[i])
end
return t, u
endARK (generic function with 1 method)Example 5.6 We use this method to numerically solve Example 5.5:
t1 = 2π
k(t) = t<=t1 ? 1 : 50
f(t,v) = -9.81 + (1/80) * k(t) * v^2
v0, T, n = 0., 10., 78
ivp = ODEProblem((v,p,t)->f(t,v), v0, (0., T))
sol = solve(ivp, Tsit5())
N = length(sol.t)
t, v = ARK( v0, f, T; tol=1e-2 )
n = length(t)
plot(sol,
title="Velocity of a skydiver",
label="Tsit5 (n=$N)", xlabel=L"t", ylabel=L"v(t)",
ylims=(-28,0))
plot!(t, v, m=3, label="ARK (n=$n)")