6.2. Runge-Kutta Methods#

Remark 6.1 (TO DO)

Improve the presentation of examples

References:

import numpy as np
# Shortcuts for some favorite commands:
from numpy import linspace
from matplotlib.pyplot import figure, plot, grid, title, xlabel, ylabel, legend

6.2.1. Introduction#

The original Runge-Kutta method is the fourth order accurate one to be described below, which is still used a lot, though with some modifications.

However, the name is now applied to a variety of methods based on a similar strategy, so first, here are a few simpler methods, all of some value, at least for small, low precision calculations.

6.2.2. Note: the methods described below are#

with comparisons back to Euler’s Method seen in Basic Concepts and Euler’s Method.

6.2.3. Euler’s Method as a Runge-Kutta method#

The simplest of all methods of this general form is Euler’s method. To set up the notation to be used below, rephrase it this way:

To get from \((t, u)\) to an approximation of \((t+h, u(t+h))\), use the approximation

\[\begin{split}\begin{split} K_1 &= h f(t, u) \\ u(t+h) &\approx u + K_1 \end{split}\end{split}\]

6.2.4. Second order Runge-Kutta methods#

We have seen that the global error of Euler’s method is \(O(h)\): it is only first order accurate. This is often insufficient, so it is more common even for small, low precision calculation to use one of several second order methods:

The Explicit Trapezoid Method (a.k.a. the Improved Euler method or Huen’s method)#

One could try to adapt the trapezoid method for integrating \(f(t)\) to solve \(du/dt = f(t)\)

\[ u(t + h) = u(t) + \int_{t}^{t + h} f(s) ds \approx u(t) + h \frac{f(t) + f(t+h)}{2} \]

to solving the ODE \(du/dt = f(t, u)\) but there is a problem that needs to be overcome:

we get

\[ u(t + h) \approx u(t) + h \frac{f(t, u(t)) + h f(t+h, u(t+h))}{2} \]

and inserting the values \(U_i \approx u(t_i)\) and so on gives

\[ U_{i+1} \approx U_i + h \frac{f(t_i, U_i) + f(t_{i+1}, U_{i+1}))}{2} \]

This is known as the Implicit Trapezoid Rule, because the value \(U_{i+1}\) that we seek appears at the right-hand side too: we only have an implicit formula for it.

On one hand, one can in fact use this formula, by solving the equation at each time step for the unknown \(U_{i+1}\); for example, one can use methods seen in earlier sections such as fixed point iteration or the secant method.

We will return to this in a later section; however, for now we get around this more simply by inserting an approximation at right — the only one we know so far, given by Euler’s Method. That is:

  • replace \(u(t+h)\) at right by the tangent line approximation \(u(t+h) \approx u(t) + hf(t, u(t))\), giving

\[ u(t + h) \approx u(t) + h \frac{f(t, u(t)) + f(t+h, u(t) + hf(t, u(t)))}{2} \]

and for the formulas in terms of the \(U_i\), replace \(U_{i+1}\) at right by \(U_{i+1} \approx U_i + h f(t_i, U_i)\), giving

\[ U_{i+1} = U_i + h \frac{f(t_i, U_i) + f(t_{i+1}, U_i + h f(t_i, U_i))}{2} \]

This is the Explicit Trapezoid Rule.

It is convenient to break this down into two stages, one for each evaluation of \(f(t, u)\):

\[\begin{split}\begin{split} K_1 &= h f(t, u) \\ K_2 &= h f(t+h, u + K_1) \\ u(t+h) &\approx u + \frac{1}{2}(K_1 + K_2) \end{split}\end{split}\]

For equal sized time steps, this leads to

Algorithm 6.1 (The Explicit Trapezoid Method)

\[\begin{split}\begin{split} U_0 &= u_0 \\ U_{i+1} &= U_i + \frac{1}{2}(K_1 + K_2), \\ &\quad \text{where} \\ K_1 &= h f(t_i, U_i) \\ K_2 &= h f(t_{i+1}, U_i + K_1) \end{split}\end{split}\]

We will see that, despite the mere first order accuracy of the Euler approximation used in getting \(K_2\), this method is second order accurate; the key is the fact that any error in the approximation used for \(f(t+h, u(t+h))\) gets multiplied by \(h\).

Exercise 1#

A) Verify that for the simple case where \(f(t, u) = f(t)\), this gives the same result as the composite trapezoid rule for integration.

B) Do one step of this method for the canonical example \(du/dt = ku\), \(u(t_0) = u_0\). It will have the form \(U_1 = G U_0\) where the growth factor \(G\) approximates the factor \(g=e^{kh}\) for the exact solution \(u(t_1) = g u(t_0)\) of the ODE.

C) Compare to \(G=1+kh\) seen for Euler’s method.

D) Use the previous result to express \(U_i\) in terms of \(U_0=u_0\), as done for Euler’s method.

def explicitTrapezoid(f, a, b, u_0, n):
    """Use the Explict Trapezoid Method (a.k.a Improved Euler)
    to solve du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0
    """
    h = (b-a)/n
    t = np.linspace(a, b, n+1)  # Note: "n" counts steps, so there are n+1 values for t.
    u = np.empty_like(t)
    u[0] = u_0
    for i in range(n):
        K_1 = f(t[i], u[i])*h
        K_2 = f(t[i]+h, u[i]+K_1)*h
        u[i+1] = u[i] + (K_1 + K_2)/2.
    return (t, u)

As always, this function can now also be imported from numericalMethods, with

from numericalMethods import explicitTrapezoid

Examples#

For all methods in this section, we will solve for versions of Example 2 and 4 in the section Basic Concepts and Euler’s Method.

\[ \frac{du}{dt} = f_1(t, u) = k u \]

with solution

\[ u(t) = u_1(t; a, u_0) = u_0 e^{t-a} \]

and

\[ \frac{du}{dt} = k(\cos(t) - u) - \sin(t) \]

with solution

\[ u(t) = u_2(t; a, u_0, k) = \cos t + c e^{-k(t-a)}, \; c = u_0 - \cos(a) \]

For comparison, the same examples are done below with Euler’s method.

def f2(t, u):
    """The simplest "genuine" ODE, (not just integration)
    The solution is u(t) = u(t; a, u_0) = u_0 exp(k(t-a))
    """
    return k*u
def u2(t, a, u_0, k):
    return u_0 * np.exp(k*(t-a))
    a = 1.
    b = 3.
    u_0 = 2.
    k = 1.5
    n = 40

    (t, U) = explicitTrapezoid(f2, a, b, u_0, n)
    u = u2(t, a, u_0, k)
    h = (b-a)/n
    figure(figsize=[14,5])
    title(f"Solving du/dt = {k}u, u({a})={u_0} by the Explicit Trapezoid Method")
    plot(t, u, "g", label="Exact solution")
    plot(t, U, ".:b", label=f"Solution with h={h:0.4}")
    legend()
    grid(True)

    figure(figsize=[14,5])
    title(f"Error")
    plot(t, u - U, '.:')
    grid(True);
../_images/d4126be2c92d7e4e8d412afc89ae30066d94cbecac4006625aa5f605d47f1fe4.png ../_images/16d6f6a40019be4fc6e7fd309fef683b62a443c22cb5dc10ee2eee24f928c5ac.png
def f4(t, u):
    """A simple more "generic" test case, with f(t, u) depending on both variables.
    The general solution is
        u(t) = u(t; a, u_0) = cos t + C e^(-k t),
        C = (u_0 - cos(a)) exp(k a)
    """
    return  k*(np.cos(t) - u) - np.sin(t)

def u4(t, a, u_0, k):
    return np.cos(t) + (u_0 - np.cos(a)) * np.exp(k*(a-t))
a = 1.
b = a + 4 * np.pi  # Two periods
u_0 = 2.
k = 2.
n = 80

(t, U) = explicitTrapezoid(f4, a, b, u_0, n)
u = u4(t, a, u_0, k)
#h = (b-a)/n
figure(figsize=[14,5])
title(f"Solving du/dt = {k}(cos(t) - u) - sin(t), u({a})={u_0} by the Explicit Trapezoid Method")
plot(t, u, "g", label="Exact solution")
plot(t, U, ".:b", label=f"Solution with h={(b-a)/n:0.4}")
legend()
grid(True)
    
figure(figsize=[14,5])
title(f"Error")
plot(t, U - u, '.:')
grid(True);
../_images/fb63f2ac0cfe93779ade07a930bbb92acf8f36b8ea9f776743bca2ad6e3b7d38.png ../_images/d2256ca707c7bf3e16a975baf466ed018cc336d40260b0609e6f85e285eacf02.png

The Explicit Midpoint Method (a.k.a. Modified Euler)#

If we start with the Midpoint Rule for integration in place of the Trapezoid Rule, we similarly get an approximation

\[ u(t + h) \approx u(t) + h f(t + h/2, u(t + h/2)) \]

This has the slight extra complication that it involves three values of \(u\) including \(u(t+h/2)\) which we are not trying to evaluate. We deal with that by making yet another approximation, using an average of \(u\) values:

\[ u(t+h/2) \approx \frac{u(t) + u(t+h)}{2} \]

leading to

\[ u(t + h) \approx u(t) + h f\left( t + h/2, \frac{u(t) + u(t+h)}{2} \right) \]

and in terms of \(U_i \approx u(t_i)\), the Implicit Midpoint Rule

\[ U_{i+1} = U_i + h f\left( t + h/2, \frac{U_i + U_{i+1}}{2} \right) \]

We will see late that this is a particularly useful method in some situations, such as long-time solutions of ODEs that describe the motion of physical systems with conservation of momentum, angular momentum and kinetic energy.

However, for now we again seek a more straightforward explicit method; using the same tangent line approximation strategy as above gives the Explicit Midpoint Rule

\[\begin{split}\begin{split} K_1 &= h f(t, u) \\ K_2 &= h f(t+h/2, u + K_1/2) \\ u(t+h) &\approx u + K_2 \end{split}\end{split}\]

and thus for equal-sized time steps

Algorithm 6.2 (The Explicit Midpoint Method)

\[\begin{split}\begin{split} U_0 &= u_0 \\ U_{i+1} &= U_i + K_2 \\ &\quad \text{where} \\ K_1 &= h f(t_i, U_i) \\ K_2 &= h f(t_{i}+h/2, U_i + K_1/2) \end{split}\end{split}\]

Exercise 2 (a lot like Exercise 1)#

A) Verify that for the simple case where \(f(t, u) = f(t)\), this give the same result as the composite midpoint rule for integration (same comment as above).

B) Do one step of this method for the canonical example \(du/dt = ku\), \(u(t_0) = u_0\). It will have the form \(U_1 = G U_0\) where the growth factor \(G\) approximates the factor \(g=e^{kh}\) for the exact solution \(u(t_1) = g u(t_0)\) of the ODE.

C) Compare to the growth factors \(G\) seen for previous methods, and to the growth factor \(g\) for the exact solution.

Exercise 3#

A) Apply Richardson extrapolation to one step of Euler’s method, using the values given by step sizes \(h\) and \(h/2\).

B) This should give a second order accurate method, so compare it to the above two methods.

def explicitMidpoint(f, a, b, u_0, n):
    """Use the Explicit Midpoint Method (a.k.a Modified Euler)
    to solve du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0
    """
    h = (b-a)/n
    t = np.linspace(a, b, n+1)  # Note: "n" counts steps, so there are n+1 values for t.
    u = np.empty_like(t)
    u[0] = u_0
    for i in range(n):
        K_1 = f(t[i], u[i])*h
        K_2 = f(t[i]+h/2, u[i]+K_1/2)*h
        u[i+1] = u[i] + K_2
    return (t, u)

Again, available for import with

from numericalMethods import explicitMidpoint

Examples#

    a = 1.
    b = 3.
    u_0 = 2.
    k = 1.5
    n = 40

    (t, U) = explicitMidpoint(f2, a, b, u_0, n)
    u = u2(t, a, u_0, k)
    figure(figsize=[14,5])
    title(f"Solving du/dt = {k}u, u({a})={u_0} by the Explicit Midpoint Method")
    plot(t, u, "g", label="Exact solution")
    plot(t, U, ".:b", label=f"Solution with h={(b-a)/n:0.4}")
    legend()
    grid(True)

    figure(figsize=[14,5])
    title(f"Error")
    plot(t, u - U, '.:')
    grid(True);
../_images/a7b6e4a8a993495399a18f9f1c9b7b16bfbc6a0772e966fa912c78fec65cf8ec.png ../_images/16d6f6a40019be4fc6e7fd309fef683b62a443c22cb5dc10ee2eee24f928c5ac.png

Observation: The errors are very similar to those for the Explicit Trapezoid Method, not “half as much and of opposite sign” as seen with integration; the exercises give a hint as to why this is so.

    a = 1.
    b = a + 4 * np.pi  # Two periods
    u_0 = 2.
    k = 2.
    n = 80

    (t, U) = explicitMidpoint(f4, a, b, u_0, n)
    u = u4(t, a, u_0, k)
    h = (b-a)/n
    figure(figsize=[14,5])
    title(f"Solving du/dt = {k}(cos(t) - u) - sin(t), u({a})={u_0} by the Explicit Midpoint Method")
    plot(t, u, "g", label="Exact solution")
    plot(t, U, ".:b", label=f"Solution with h={(b-a)/n:0.4}")
    legend()
    grid(True)
    
    figure(figsize=[12,5])
    title(f"Error")
    plot(t, u - U, '.:')
    grid(True);
../_images/190bf6b2cba83aef97bd08e6f5a66c040df142c443e5bc2d60f42ac6374d7e4c.png ../_images/fe7065346ef9db0c8c99d09f53c0f925f5d00affd504c3c3a3ccad956ef2a446.png

Observation: This time, the errors are slightly better than for the Explicit Trapezoid Method but still not “half as much “ as seen with integration; this is because this equation has a mix of integration (the “sin” and “cos” parts) and exponential growth (the “Ku” part.)

6.2.5. The “Classical”, Fourth Order Accurate, Runge-Kutta Method#

This is the original Runge-Kutta method:

Algorithm 6.3 (Runge-Kutta)

\[\begin{split}\begin{split} K_1 &= h f(t, u) \\ K_2 &= h f(t + h/2, u + K_1/2) \\ K_3 &= h f(t + h/2, u + K_2/2) \\ K_4 &= h f(t + h, u + K_3) \\ u(t+h) &\approx u + \frac{1}{6}(K_1 + 2 K_2 + 2 K_3 + K_4) \end{split}\end{split}\]

The derivation of this is far more complicated than those above, and is omitted. For now, we will instead assess its accuracy “a postiori”, through the next exercise and some examples.

Exercise 4#

A) Verify that for the simple case where \(f(t, u) = f(t)\), this gives the same result as the composite ximpson’s Rule for integration.

B) Do one step of this method for the canonical example \(du/dt = ku\), \(u(t_0) = u_0\). It will have the form \(U_1 = G U_0\) where the growth factor \(G\) approximates the factor \(g=e^{kh}\) for the exact solution \(u(t_1) = g u(t_0)\) of the ODE.

C) Compare to the growth factors \(G\) seen for previous methods, and to the growth factor \(g\) for the exact solution.

def rungeKutta(f, a, b, u_0, n):
    """Use the (classical) Runge-Kutta Method
    to solve du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0
    """
    h = (b-a)/n
    t = np.linspace(a, b, n+1)  # Note: "n" counts steps, so there are n+1 values for t.
    u = np.empty_like(t)
    u[0] = u_0
    for i in range(n):
        K_1 = f(t[i], u[i])*h
        K_2 = f(t[i]+h/2, u[i]+K_1/2)*h
        K_3 = f(t[i]+h/2, u[i]+K_2/2)*h
        K_4 = f(t[i]+h, u[i]+K_3)*h
        u[i+1] = u[i] + (K_1 + 2*K_2 + 2*K_3 + K_4)/6
    return (t, u)

Yet again, available for import with

from numericalMethods import rungeKutta

Examples#

    a = 1.
    b = 3.
    u_0 = 2.
    k = 1.5
    n = 20

    (t, U) = rungeKutta(f2, a, b, u_0, n)
    u = u2(t, a, u_0, k)
    
    figure(figsize=[14,5])
    title(f"Solving du/dt = {k}u, u({a})={u_0} by the Runge-Kutta Method")
    plot(t, u, "g", label="Exact solution")
    plot(t, U, ".:b", label=f"Solution with h={(b-a)/n:0.4}")
    legend()
    grid(True)

    figure(figsize=[14,5])
    title(f"Error")
    plot(t, u - U, ".:")
    grid(True);
../_images/001bec4998c197a529ab8ad1452efab92e6e8ccb9fc8bbf1c9f07baa023ba01d.png ../_images/2e687ec443b48f142c86afae398003ad87ced2bd323e5ac7aa7897265940f817.png
    a = 1.
    b = a + 4 * np.pi  # Two periods
    u_0 = 2.
    k = 2.
    n = 40

    (t, U) = rungeKutta(f4, a, b, u_0, n)
    u = u4(t, a, u_0, k)
    
    figure(figsize=[14,5])
    title(f"Solving du/dt = {k}(cos(t) - u) - sin(t), u({a})={u_0} by the Runge-Kutta Method")
    plot(t, u, "g", label="Exact solution")
    plot(t, U, ".:b", label=f"Solution with h={(b-a)/n:0.4}")
    legend()
    grid(True)
    
    figure(figsize=[14,5])
    title(f"Error")
    plot(t, u - U, ".:")
    grid(True);
../_images/cb28d9f77384337ad42408906579fe954e589efb20850bea05a8ba7ee7819585.png ../_images/e04c9ba1e741ba37f7261cd3e3499bacfa86baf093fb1ca0aaabe0ae7128c557.png

6.2.6. For comparison: the above examples done with Euler’s Method#

from numericalMethods import eulerMethod
    a = 1.
    b = 3.
    u_0 = 2.
    k = 1.5
    n = 80

    (t, U) = eulerMethod(f2, a, b, u_0, n)
    u = u2(t, a, u_0, k)
    figure(figsize=[14,5])
    title(f"Solving du/dt = {k}u, u({a})={u_0} by Euler's method")
    plot(t, u, "g", label="Exact solution")
    plot(t, U, ".:b", label=f"Euler's answer with h={(b-a)/n:0.4}")
    legend()
    grid(True)

    figure(figsize=[14,5])
    title(f"Error")
    plot(t, u - U, ".:")
    grid(True);
../_images/1c7c0679c08e01698f6529973ce4b4984ccecffd942b5995bb394bef79709c8b.png ../_images/5fd5cee9afc0e299bffdac2c3570ee7ff55634e28f12b9f22da84745bae0f4f7.png
    a = 1.
    b = a + 4 * np.pi  # Two periods
    u_0 = 2.
    k = 2.
    n = 160

    (t, U) = eulerMethod(f4, a, b, u_0, n)
    u = u4(t, a, u_0, k)

    figure(figsize=[12,5])
    title(f"Solving du/dt = {k}(cos(t) - u) - sin(t), u({a})={u_0} by Euler's method")
    plot(t, u, "g", label="Exact solution")
    plot(t, U, ".:b", label=f"Euler's answer with h={(b-a)/n:0.4}")
    legend()
    grid(True)
    
    figure(figsize=[14,5])
    title(f"Error")
    plot(t, u - U, ".:")
    grid(True);
../_images/0ba82b3bab04fab8b566a3a4d25a024446ddeebaed58162f52bd248c0a4764da.png ../_images/e690da7c39f8e3b6dcc9be3a690d026993f54ebb0f7d071e22d2781966051a3c.png