51. Numerical solution of ordinary differential equations: High order Runge-Kutta methods#
André Massing
Date: March 17, 2021
If you want to have a nicer theme for your jupyter notebook, download the cascade stylesheet file tma4125.css and execute the next cell:
import matplotlib.pyplot as plt
import numpy as np
from IPython.core.display import HTML
from numpy import pi
from numpy.linalg import norm, solve
def css_styling():
try:
with open("tma4125.css") as f:
styles = f.read()
return HTML(styles)
except FileNotFoundError:
pass # Do nothing
# Comment out next line and execute this cell to restore the default notebook style
css_styling()
As always, we start by import some important Python modules.
%matplotlib inline
# Use a funny plotting style
plt.xkcd()
newparams = {
"figure.figsize": (6.0, 6.0),
"axes.grid": True,
"lines.markersize": 8,
"lines.linewidth": 2,
"font.size": 14,
}
plt.rcParams.update(newparams)
52. Runge-Kutta Methods#
In the previous lectures we introduced Euler’s method and Heun’s method as particular instances of the One Step Methods, and we presented the general error theory for one step method.
In this Lecture, we introduce a large family of the one step methods which go under the name Runge-Kutta methods (RKM). We will see that Euler’s method and Heun’s method are instance of RKMs.
52.1. Derivation of Runge-Kutta Methods#
For a given time interval \(I_i = [t_i, t_{i+1}]\) we want to compute \(y_{i+1}\) assuming that \(y_i\) is given. Start from the exact expression
The idea is now to approximate the integral by some quadrature rule \(\mathrm{Q}[\cdot](\{\xi_j\}_{j=1}^s,\{b_j\}_{j=1}^s)\) defined on \(I_i\). Then we get
Now we can define \(\{c_j\}_{j=1}^s\) such that \(\xi_j = t_{i} + c_j \tau\) for \(j=1,\ldots,s\)
52.2. Exercise 1: A first condition on \(b_j\)#
Question: What value do you expect for \(\sum_{j=1}^s b_{j}\)?
Choice A: \(\sum_{j=1}^s b_{j} = \tau\)
Choice B: \(\sum_{j=1}^s b_{j} = 0\)
Choice C: \(\sum_{j=1}^s b_{j} = 1\)
In contrast to pure numerical integration, we don’t know the values of \(y(\xi_j)\). Again, we could use the same idea to approximate
But then again we get a closure problem if we choose new quadrature points. The idea is now to not introduce even more new quadrature points but to use same \(y(\xi_j)\) to avoid the closure problem. Note that this leads to an approximation of the integrals \(\int_{t_i}^{t_i+c_j \tau}\) with possible nodes outside of \([t_i, t_i + c_j \tau ]\).
This leads us to
where we set \( c_j \tilde{a}_{jl} = a_{jl}\).
52.3. Exercise 2: A first condition on \(a_{jl}\)#
Question: What value do you expect for \(\sum_{l=1}^s a_{jl}\)?
Choice A: \(\sum_{l=1}^s a_{jl} = \tfrac{1}{c_j}\)
Choice B: \(\sum_{l=1}^s a_{jl} = c_j \)
Choice C: \(\sum_{l=1}^s a_{jl} = 1 \)
Choice D: \(\sum_{l=1}^s a_{jl} = \tau \)
52.4. Definition 1: Runge-Kutta methods#
Given \(b_j\), \(c_j\), and \(a_{jl}\) for \(j,l = 1,\ldots s\), the Runge-Kutta method is defined by the recipe
Runge-Kutta schemes are often specified in the form of a Butcher table:
If \(a_{ij} = 0\) for \(j \geqslant i\) the Runge-Kutta method is called explicit. (Why?)
Note that in the final step, all the function evaluation we need to perform have already been performed when computing the stages \(Y_j\).
Therefore one often rewrite the scheme by introducing stage derivatives
so the resulting scheme will be (swapping index \(l\) and \(j\))
52.5. Exercise 3: Butcher table for the explicit Euler#
Write down the Butcher table for the explicit Euler.
Solution. Define \(k_1 = f(t_i, y_i) = f(t_i + 0\cdot \tau, y_i + \tau \cdot 0 \cdot k_1)\). Then the explicit Euler step \(y_{i+1} = y_i + \tau k_1 = y_i + \tau \cdot 1 \cdot k_1\), and thus the Butcher table is given by
52.6. Exercise 4: The improved explicit Euler method#
We formally derive the explicit midpoint rule or improved explicit Euler method. Applying the midpoint rule to our integral representatio yields
Since we cannot determine the value \(y(t_i + \tfrac{1}{2}\tau)\) from this system, we approximate it using a half Euler step
leading to the scheme
a) Is this a one-step function? Can you define the increment function \(\Phi\)?
Solution. Yes it is, and it’s increment function is given by
b) Can you rewrite this as a Runge-Kutta method? If so, determine the Butcher table of it.
Solution. Define \(k_1\) and \(k_2\) as follows,
Then \begin{equation} y_{i+1} = y_i + \tau k_2 \label{ode:eq:improved_euler_k2} \tag{19} \end{equation}
Thus the Butcher table is given by
52.7. Implementation of explicit Runge-Kutta methods#
Below you will find the implementation a general
solver class
Explicit_Runge_Kutta
which at its initialization takes
in a Butcher table and has __call__ function
def __call__(self, y0, f, t0, T, n):
and can be used like this
# Define Butcher table
a = np.array([[0, 0, 0],
[1.0/3.0, 0, 0],
[0, 2.0/3.0, 0]])
b = np.array([1.0/4.0, 0, 3.0/4.0])
c = np.array([0, 1.0/3.0, 2.0/3.0])
# Create solver
rk3 = Explicit_Runge_Kutta(a, b, c)
# Solve problem (applies __call__ function)
ts, ys = rk3(y0, t0, T, f, Nmax)
The complete implementation is given here:
class Explicit_Runge_Kutta:
def __init__(self, a, b, c):
self.a = a
self.b = b
self.c = c
def __call__(self, y0, t0, T, f, Nmax):
# Extract Butcher table
a, b, c = self.a, self.b, self.c
# Stages
s = len(b)
ks = [np.zeros_like(y0, dtype=np.double) for s in range(s)]
# Start time-stepping
ys = [y0]
ts = [t0]
dt = (T - t0) / Nmax
while ts[-1] < T:
t, y = ts[-1], ys[-1]
# Compute stages derivatives k_j
for j in range(s):
t_j = t + c[j] * dt
dY_j = np.zeros_like(y, dtype=np.double)
for l in range(j):
dY_j += dt * a[j, l] * ks[l]
ks[j] = f(t_j, y + dY_j)
# Compute next time-step
dy = np.zeros_like(y, dtype=np.double)
for j in range(s):
dy += dt * b[j] * ks[j]
ys.append(y + dy)
ts.append(t + dt)
return (np.array(ts), np.array(ys))
52.8. Example 1: Implementation and testing of the improved Euler method#
We implement the improved explicit Euler from above and plot the analytical and the numerical solution. Finally, we determine the convergence order.
def compute_eoc(y0, t0, T, f, Nmax_list, solver, y_ex):
errs = []
for Nmax in Nmax_list:
ts, ys = solver(y0, t0, T, f, Nmax)
ys_ex = y_ex(ts)
errs.append(np.abs(ys - ys_ex).max())
print(f"For Nmax = {Nmax:3}, max ||y(t_i) - y_i||= {errs[-1]:.3e}")
errs = np.array(errs)
Nmax_list = np.array(Nmax_list)
dts = (T - t0) / Nmax_list
eocs = np.log(errs[1:] / errs[:-1]) / np.log(dts[1:] / dts[:-1])
return errs, eocs
# Define Butcher table for improved Euler
a = np.array([[0, 0], [0.5, 0]])
b = np.array([0, 1])
c = np.array([0, 0.5])
# Define rk2
rk2 = Explicit_Runge_Kutta(a, b, c)
t0, T = 0, 1
y0 = 1
lam = 1
# rhs of IVP
f = lambda t, y: lam * y
# Exact solution to compare against
y_ex = lambda t: y0 * np.exp(lam * (t - t0))
# EOC test
Nmax_list = [4, 8, 16, 32, 64, 128]
errs, eocs = compute_eoc(y0, t0, T, f, Nmax_list, rk2, y_ex)
print(errs)
print(eocs)
For Nmax = 4, max ||y(t_i) - y_i||= 2.343e-02
For Nmax = 8, max ||y(t_i) - y_i||= 6.441e-03
For Nmax = 16, max ||y(t_i) - y_i||= 1.688e-03
For Nmax = 32, max ||y(t_i) - y_i||= 4.322e-04
For Nmax = 64, max ||y(t_i) - y_i||= 1.093e-04
For Nmax = 128, max ||y(t_i) - y_i||= 2.749e-05
[2.34261385e-02 6.44058991e-03 1.68830598e-03 4.32154479e-04
1.09316895e-04 2.74901378e-05]
[1.86285442 1.93161644 1.96595738 1.98303072 1.99153035]
52.9. Exercise 5#
If you have time, determine the experimental order of convergence EOC for the following methods:
Heun’s method with 3 stages $\( \begin{array}{c|c c c} 0 & 0 & 0 & 0 \\ \tfrac{1}{3} & \tfrac{1}{3} & 0 & 0 \\ \tfrac{2}{3} & 0 & \tfrac{2}{3} & 0 \\ \hline & \tfrac{1}{4} & 0 & \tfrac{3}{4} \end{array} \)$
# Insert code here
The classical Runge-Kutta method with 4 stages $\( \begin{array}{c|c c c c} 0 & 0 & 0 & 0 & 0 \\ \tfrac{1}{2} & \tfrac{1}{2} & 0 & 0 & 0 \\ \tfrac{1}{2} & 0 & \tfrac{1}{2} & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ \hline & \tfrac{1}{6} & \tfrac{1}{3}& \tfrac{1}{3} & \tfrac{1}{6} \end{array} \)$
# Insert code here
52.10. Order conditions for Runge-Kutta Methods#
The convergence theorem for one-step methods gave us some necessary conditions to guarantee that a method is convergent order of \(p\):
“consistency order \(p\)” + “Increment function satisfies a Lipschitz condition” \(\Rightarrow\) “convergence order \(p\).
To put it in slightly different words:
“local truncation error behaves like \(\mathcal{O}(\tau^{p+1})\)” + “Increment function satisfies a Lipschitz condition” \(\Rightarrow\) “global truncation error behaves like \(\mathcal{O}(\tau^{p})\)”
It turns out that for \(f\) is at least \(C^1\) with respect to all its arguments then the increment function \(\Phi\) associated with any Runge-Kutta methods satisfies a Lipschitz condition.
52.11. Theorem 1: Order conditions for Runge-Kutta methods#
A Runge - Kutta method has consistency order $p$ if and only if all the conditions up to and including $p$ in the table below are satisfied.where sums are taken over all the indices from 1 to \(s\).
Proof. Without proof.