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,\dots ,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}=\frac{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,\dots 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⩾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+\frac{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 $\mathrm{\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}{cccc}0& 0& 0& 0\\ \frac{1}{3}& \frac{1}{3}& 0& 0\\ \frac{2}{3}& 0& \frac{2}{3}& 0\\ & \frac{1}{4}& 0& \frac{3}{4}\end{array}$$

# Insert code here

The classical Runge-Kutta method with 4 stages $$\begin{array}{ccccc}0& 0& 0& 0& 0\\ \frac{1}{2}& \frac{1}{2}& 0& 0& 0\\ \frac{1}{2}& 0& \frac{1}{2}& 0& 0\\ 1& 0& 0& 1& 0\\ & \frac{1}{6}& \frac{1}{3}& \frac{1}{3}& \frac{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 $\mathrm{\Phi }$ associated with any Runge-Kutta methods satisfies a Lipschitz condition.

52.11. Theorem 1: Order conditions for Runge-Kutta methods#

where sums are taken over all the indices from 1 to $s$.

Proof. Without proof.

52.12. Theorem 2: Convergence theorem for Runge-Kutta methods#