Runge–Kutta methods - Fundamentals of Numerical Computation

We come now to one of the major and most-used types of methods for initial-value problems: Runge–Kutta (RK) methods.^[1] They are one-step methods in the sense of (6.2.7), though they are not often written in that form. RK methods boost the accuracy past first order by evaluating the ODE function $f(t,u)$ more than once per time step.

6.4.1A second-order method¶

Consider a series expansion of the exact solution to $u'=f(t,u)$ ,

\hat{u}(t_{i+1}) = \hat{u}(t_i) + h \hat{u}'(t_i) + \frac{1}{2}h^2 \hat{u}''(t_i) + O(h^3) .

(6.4.1)

If we replace $\hat{u}'$ by $f$ and keep only the first two terms on the right-hand side, we would obtain the Euler method. To get more accuracy we will need to compute or estimate the third term as well. Note that

\hat{u}'' = f' = \frac{d f}{d t} = \frac{\partial f}{\partial t} + \frac{\partial f}{\partial u} \frac{d u}{d t} = f_t + f_u f,

(6.4.2)

where we have applied the multidimensional chain rule to the derivative, because both of the arguments of $f$ depend on $t$ . Using this expression in (6.4.1), we obtain

\begin{split} \hat{u}(t_{i+1}) &= \hat{u}(t_i) \\ & + h\left[f\bigl(t_i,\hat{u}(t_i)\bigr) + \frac{h}{2}f_t\bigl(t_i,\hat{u}(t_i)\bigr) + \frac{h}{2}f\bigl(t_i,\hat{u}(t_i)\bigr)\,f_u\bigl(t_i,\hat{u}(t_i)\bigr)\right] \\ &+ O(h^3). \end{split}

(6.4.3)

We have no desire to calculate and then code those partial derivatives of $f$ directly; an approximate approximation is called for. Observe that

\begin{split} f\bigl(t_i+\alpha,\hat{u}(t_i)+\beta\bigr) & = f\bigl(t_i,\hat{u}(t_i)\bigr) \\ & + \alpha f_t\bigl(t_i,\hat{u}(t_i)\bigr) + \beta f_u\bigl(t_i,\hat{u}(t_i)\bigr) \\ & + O\bigl(\alpha^2 + |\alpha\beta| + \beta^2\bigr). \end{split}

(6.4.4)

Matching this expression to the term in brackets in (6.4.3), it seems natural to select $\alpha = h/2$ and $\beta = \frac{1}{2}h f\bigl(t_i,\hat{u}(t_i)\bigr)$ . Doing so, we find

\hat{u}(t_{i+1}) = \hat{u}(t_i) + h\left[f\bigl(t_i+\alpha,\hat{u}(t_i)+\beta\bigr)\right] + O(h\alpha^2 + h|\alpha \beta| + h\beta^2 + h^3).

(6.4.5)

Truncation of the series here results in a new one-step method.

Thanks to the definitions above of α and β, the omitted terms are of size

O(h\alpha^2 + h|\alpha \beta| + h\beta^2 + h^3) = O(h^3).

(6.4.7)

Therefore, $h\tau_{i+1}=O(h^3)$ , and the order of accuracy of improved Euler is two.

6.4.2Implementation¶

Runge–Kutta methods are called multistage methods. We can see why if we interpret (6.4.6) from the inside out. In the first stage, the method takes an Euler half-step to time $t_i+\frac{1}{2}h$ :

\begin{split} k_1 &= h f(t_i,u_i), \\ v &= u_i + \tfrac{1}{2}k_1. \end{split}

(6.4.8)

The second stage employs an Euler-style strategy over the whole time step, but using the value from the first stage to get the slope:

\begin{split} k_2 &= h f\left(t_i+\tfrac{1}{2}h,v\right),\\ {u}_{i+1} &= u_i + k_2. \end{split}

(6.4.9)

Our implementation of IE2 is shown in Function 6.4.1.

Algorithm 6.4.1 (ie2)

Julia

MATLAB

Python

Improved Euler method for an IVP

ie2.jl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
"""
    ie2(ivp, n)

Apply the Improved Euler method to solve the given IVP using `n`
time steps. Returns a vector of times and a vector of solution
values.
"""
function ie2(ivp, n)
    # Time discretization.
    a, b = ivp.tspan
    h = (b - a) / n
    t = [a + i * h for i in 0:n]

    # Initialize output.
    u = fill(float(ivp.u0), n+1)

    # Time stepping.
    for i in 1:n
        uhalf = u[i] + h / 2 * ivp.f(u[i], ivp.p, t[i])
        u[i+1] = u[i] + h * ivp.f(uhalf, ivp.p, t[i] + h / 2)
    end
    return t, u
end

Improved Euler method for an IVP

ie2.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
function [t, u] = ie2(du_dt, tspan, u0, n)
% IE2    Improved Euler method for an IVP.
% Input:
%   du_dt   defines f in u'(t) = f(t, u) 
%   tspan   endpoints of time interval (2-vector)
%   u0      initial value (m-vector)
%   n       number of time steps (integer)
% Output:
%   t       selected nodes (vector, length n+1)
%   u       solution values (array, n+1 by m)

    % Define the time discretization.
    a = tspan(1);  b = tspan(2);
    h = (b - a) / n;
    t = a + (0:n)' * h;

    % Initialize solution array.
    u = zeros(length(u0), n+1);
    u(:, 1) = u0(:);

    % Time stepping.
    for i = 1:n 
        uhalf = u(:, i) + h/2 * du_dt(t(i), u(:, i));
        u(:, i+1) = u(:, i) + h * du_dt(t(i) + h/2, uhalf);
    end

    u = u.';    % conform to MATLAB output convention
end

Improved Euler method for an IVP

ie2.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
def ie2(du_dt, tspan, u0, n):
    """
    ie2(du_dt, tspan, u0, n)

    Apply the Improved Euler method to solve the vector-valued IVP u'=du_dt(u,p,t) over the
    interval tspan with u(tspan[1])=u0, using n subintervals/steps. Returns a vector
    of times and a vector of solution values/vectors.
    """
    # Time discretization.
    a, b = tspan
    h = (b - a) / n
    t = np.linspace(a, b, n + 1)

    # Initialize output.
    u = np.tile(np.array(u0), (n+1, 1))

    # Time stepping.
    for i in range(n):
        uhalf = u[i] + h / 2 * du_dt(t[i], u[i])
        u[i+1] = u[i] + h * du_dt(t[i] + h / 2, uhalf)

    return t, u.T

6.4.3More Runge–Kutta methods¶

While it is interesting to interpret IE2 as a pair of Euler-like steps, the Taylor series derivation is the only way to see that it will be more accurate than Euler, and it is also the path to deriving other methods. Moving to higher orders of accuracy requires introducing additional stages, each having free parameters so that more terms in the series may be matched. The amount of algebra grows rapidly in size and complexity, though there is a sophisticated theory for keeping track of it. We do not give the derivation details.

A generic $s$ -stage RK method takes the form

\begin{split} k_1 &= h f(t_i,u_i),\\ k_2 &= h f(t_i+c_1h,u_i + a_{11}k_1),\\ k_3 &= h f(t_i+c_2h, u_i + a_{21}k_1 + a_{22}k_2),\\ &\vdots\\ k_s &= h f(t_i + c_{s-1}h, u_i + a_{s-1,1}k_1 + \cdots + a_{s-1,s-1}k_{s-1}),\\ \mathbf{u}_{i+1} &= u_i + b_1k_1 + \cdots + b_s k_s. \end{split}

(6.4.10)

This recipe is completely determined by the number of stages $s$ and the constants $a_{ij}$ , $b_j$ , and $c_i$ . Often an RK method is presented as just a table of these numbers, as in

\begin{array}{r|ccccc} 0 & & & & & \\ c_1 & a_{11} & & &\\ c_2 & a_{21} & a_{22} & & &\\ \vdots & \vdots & & \ddots & &\\ c_{s-1} & a_{s-1,1} & \cdots & & a_{s-1,s-1}&\\[1mm] \hline \rule{0pt}{2.25ex} & b_1 & b_2 & \cdots & b_{s-1} & b_s \end{array}

(6.4.11)

For example, IE2 is given by

\begin{array}{r|cc} \rule{0pt}{2.75ex}0 & & \\ \rule{0pt}{2.75ex}\frac{1}{2} & \frac{1}{2} &\\[1mm] \hline \rule{0pt}{2.75ex}& 0 & 1 \end{array}

(6.4.12)

Here are two more two-stage, second-order methods, modified Euler and Heun’s method, respectively:

\begin{array}{r|cc} \rule{0pt}{2.75ex}0 & & \\ \rule{0pt}{2.75ex}1 & 1 &\\[1mm] \hline \rule{0pt}{2.75ex}& \frac{1}{2} & \frac{1}{2} \end{array} \qquad \qquad \begin{array}{r|cc} \rule{0pt}{2.75ex} 0 & & \\ \rule{0pt}{2.75ex} \frac{2}{3} & \frac{2}{3} &\\[1mm] \hline \rule{0pt}{2.75ex} & \frac{1}{4} & \frac{3}{4} \end{array}

(6.4.13)

The most commonly used RK method, and perhaps the most popular IVP method of all, is the fourth-order one given by

\begin{array}{r|cccc} \rule{0pt}{2.75ex}0 & & & & \\ \rule{0pt}{2.75ex}\frac{1}{2} & \frac{1}{2} & & &\\ \rule{0pt}{2.75ex}\frac{1}{2} & 0 & \frac{1}{2} & &\\ \rule{0pt}{2.75ex}1 & 0 & 0 & 1\\[1mm] \hline \rule{0pt}{2.75ex}& \frac{1}{6} & \frac{1}{3} & \frac{1}{3} & \frac{1}{6} \end{array}

(6.4.14)

This formula is often referred to as the fourth-order RK method, even though there are many others, and we refer to it as RK4. Written out, the recipe is as follows.

Our implementation is given in Function 6.4.2.

Algorithm 6.4.2 (rk4)

Julia

MATLAB

Python

Fourth-order Runge-Kutta for an IVP

rk4.jl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
"""
    rk4(ivp, n)

Apply the common Runge-Kutta 4th order method to solve the given
IVP using `n` time steps. Returns a vector of times and a vector of
solution values.
"""
function rk4(ivp, n)

    # Time discretization.
    a, b = ivp.tspan
    h = (b - a) / n
    t = [a + i * h for i in 0:n]

    # Initialize output.
    u = fill(float(ivp.u0), n+1)

    # Time stepping.
    for i in 1:n
        k₁ = h * ivp.f(u[i], ivp.p, t[i])
        k₂ = h * ivp.f(u[i] + k₁ / 2, ivp.p, t[i] + h / 2)
        k₃ = h * ivp.f(u[i] + k₂ / 2, ivp.p, t[i] + h / 2)
        k₄ = h * ivp.f(u[i] + k₃, ivp.p, t[i] + h)
        u[i+1] = u[i] + (k₁ + 2(k₂ + k₃) + k₄) / 6
    end
    return t, u
end

Fourth-order Runge-Kutta for an IVP

rk4.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
function [t, u] = rk4(du_dt, tspan, u0, n)
% RK4    Fourth-order Runge-Kutta for an IVP.
% Input:
%   du_dt   defines f in u'(t) = f(t, u) 
%   tspan   endpoints of time interval (2-vector)
%   u0      initial value (m-vector)
%   n       number of time steps (integer)
% Output:
%   t       selected nodes (vector, length n+1)
%   u       solution values (array, n+1 by m)

    % Define the time discretization.
    a = tspan(1);  b = tspan(2);
    h = (b - a) / n;
    t = a + (0:n)' * h;

    % Initialize solution array. 
    u = zeros(length(u0), n+1);
    u(:, 1) = u0(:);

    % Time stepping.
    for i = 1:n
      k1 = h * du_dt( t(i),       u(:, i)       );
      k2 = h * du_dt( t(i) + h/2, u(:, i) + k1/2);
      k3 = h * du_dt( t(i) + h/2, u(:, i) + k2/2);
      k4 = h * du_dt( t(i) + h,   u(:, i) + k3  );
      u(:, i+1) = u(:, i) + (k1 + 2*(k2 + k3) + k4) / 6;
    end

    u = u.';    % conform to MATLAB output convention
end

Fourth-order Runge-Kutta for an IVP

rk4.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
def rk4(du_dt, tspan, u0, n):
    """
    rk4(du_dt, tspan, u0, n)

    Apply "the" Runge-Kutta 4th order method to solve the vector-valued IVP u'=du_dt(u,p,t)
    over the interval tspan with u(tspan[1])=u0, using n subintervals/steps.
    Return a vector of times and a vector of solution values/vectors.
    """
    # Time discretization.
    a, b = tspan
    h = (b - a) / n
    t = np.linspace(a, b, n + 1)

    # Initialize output.
    u = np.tile(np.array(u0), (n+1, 1))

    # Time stepping.
    for i in range(n):
        k1 = h * du_dt(t[i], u[i])
        k2 = h * du_dt(t[i] + h / 2, u[i] + k1 / 2)
        k3 = h * du_dt(t[i] + h / 2, u[i] + k2 / 2)
        k4 = h * du_dt(t[i] + h, u[i] + k3)
        u[i+1] = u[i] + (k1 + 2 * (k2 + k3) + k4) / 6

    return t, u.T

Example 6.4.1 (Convergence of Runge–Kutta methods)

Julia

MATLAB

Python

Example 6.4.1

We solve the IVP

u'=\sin[(u+t)^2], \quad t \in [0,4], \quad u(0)=-1.

using OrdinaryDiffEq
f(u, p, t) = sin((t + u)^2)
tspan = (0.0, 4.0)
u₀ = -1.0
ivp = ODEProblem(f, u₀, tspan)

ODEProblem with uType Float64 and tType Float64. In-place: false
Non-trivial mass matrix: false
timespan: (0.0, 4.0)
u0: -1.0

We use a DifferentialEquations solver to construct an accurate approximation to the exact solution.

u_ref = solve(ivp, Tsit5(), reltol=1e-14, abstol=1e-14);

Now we perform a convergence study of our two Runge–Kutta implementations.

n = [round(Int, 2 * 10^k) for k in 0:0.5:3]
err = zeros(length(n), 2)
for (k, n) in enumerate(n)
    t, u = FNC.ie2(ivp, n)
    err[k, 1] = norm(u_ref.(t) - u, Inf)
    t, u = FNC.rk4(ivp, n)
    err[k, 2] = norm(u_ref.(t) - u, Inf)
end
@pt :header=["n", "IE2 error", "RK4 error"] [n err]

The amount of computational work at each time step is assumed to be proportional to the number of stages. Let’s compare on an apples-to-apples basis by using the number of $f$ -evaluations on the horizontal axis.

using Plots
plot([2n 4n], err;
    m=3, label=["IE2" "RK4"], legend=:bottomleft,
    xaxis=(:log10, "f-evaluations"),  yaxis=(:log10, "inf-norm error"),
    title="Convergence of RK methods")

plot!(2n, 0.1 * err[end,1] * (n / n[end]) .^ (-2), l=:dash, label=L"O(n^{-2})")
plot!(4n, 0.1 * err[end,2] * (n / n[end]) .^ (-4), l=:dash, label=L"O(n^{-4})")

The fourth-order variant is more efficient in this problem over a wide range of accuracy.

Example 6.4.1

We solve the IVP

u'=\sin[(u+t)^2], \quad t \in [0,4], \quad u(0)=-1.

f = @(t, u) sin((t + u)^2);
u0 = -1;
a = 0;  b = 4;

We use a built-in solver to construct an accurate approximation to the exact solution.

opt = odeset('RelTol', 1e-13, 'AbsTol', 1e-13);
sol_ref = ode45(f, [a, b], u0, opt);
u_ref = @(t) deval(sol_ref, t)';

Now we perform a convergence study of our two Runge–Kutta implementations.

n = round(2 * 10.^(0:0.5:3)');
err = zeros(length(n), 2);
for i = 1:length(n)
    [t, u] = ie2(f, [a, b], u0, n(i));
    err(i, 1) = norm(u_ref(t) - u, Inf);
    [t, u] = rk4(f, [a, b], u0, n(i));
    err(i, 2) = norm(u_ref(t) - u, Inf);
end

table(n, err(:, 1), err(:, 2), variableNames=["n", "IE2 error", "RK4 error"])

clf, loglog([2*n 4*n], err, '-o')
hold on
loglog(2*n, 1e-5 * (n / n(end)) .^ (-2), '--')
loglog(4*n, 1e-10 * (n / n(end)) .^ (-4), '--')
xlabel("f-evaluations");  ylabel("inf-norm error")
title("Convergence of RK methods")
legend("IE2", "RK4", "O(n^{-2})", "O(n^{-4})", "location", "southwest");

The fourth-order variant is more efficient in this problem over a wide range of accuracy.

Example 6.4.1

We solve the IVP

u'=\sin[(u+t)^2], \quad t \in [0,4], \quad u(0)=-1.

We use the solve_ivp function to create an accurate reference solution.

from scipy.integrate import solve_ivp
du_dt = lambda t, u: sin((t + u)**2)
tspan = (0.0, 4.0)
u0 = -1.0
sol = solve_ivp(du_dt, tspan, [u0], dense_output=True, atol=1e-13, rtol=1e-13)
u_ref = sol.sol

Now we perform a convergence study of our two Runge–Kutta implementations.

n = array([int(2 * 10**k) for k in linspace(0, 3, 7)])
err = {"IE2" : [], "RK4" : []}
results = PrettyTable(["n", "IE2 error", "RK4 error"])
for i in range(len(n)):
    t, u = FNC.ie2(du_dt, tspan, u0, n[i])
    err["IE2"].append( abs(u_ref(4)[0] - u[0][-1]) )
    t, u = FNC.rk4(du_dt, tspan, u0, n[i])
    err["RK4"].append( abs(u_ref(4)[0] - u[0][-1]) )
    results.add_row([n[i], err["IE2"][-1], err["RK4"][-1]])

print(results)

+------+------------------------+------------------------+
|  n   |       IE2 error        |       RK4 error        |
+------+------------------------+------------------------+
|  2   |   1.7690264118810441   |   0.8206513302232612   |
|  6   |   0.5126838225133257   |   0.7919245473433536   |
|  20  |  0.002966971266360252  | 4.0177650143746746e-05 |
|  63  | 0.00021416270501584123 | 3.581705267929891e-07  |
| 200  | 1.951309601522233e-05  | 3.326041442264227e-09  |
| 632  | 1.9058382192405077e-06 | 3.2657876403163755e-11 |
| 2000 | 1.8883901087285437e-07 | 2.6711965972481266e-13 |
+------+------------------------+------------------------+

loglog(2 * n, err["IE2"], "-o", label="IE2")
loglog(4 * n, err["RK4"], "-o", label="RK4")
plot(2 * n, 0.5 * err["IE2"][-1] * (n / n[-1])**(-2), "--", label="2nd order")
plot(4 * n, 0.5 * err["RK4"][-1] * (n / n[-1])**(-4), "--", label="4th order")

xlabel("f-evaluations"),  ylabel("inf-norm error")
legend()
title("Convergence of RK methods");

The fourth-order variant is more efficient in this problem over a wide range of accuracy.

6.4.4Efficiency¶

As with rootfinding and integration, the usual point of view is that evaluations of $f$ are the only significant computations and are therefore to be minimized in number. One of the most important characteristics of a multistage method is that each stage requires an evaluation of $f$ ; that is, a single time step of an $s$ -stage method requires $s$ evaluations of $f$ .

The error decreases geometrically as $s$ is incremented, so trading a stage for an increase in order is a good deal. But $s=5$ , 6, or 7 gives a maximal order of accuracy of $s-1$ ; this decreases to $s-2$ for $s=8$ and $s=9$ , etc. Fourth order is considered adequate and the sweet spot for many applications.

6.4.5Exercises¶

Exercise 6.4.9

⌨ For each IVP, compute the solution using Function 6.4.2. (i) Plot the solution for $n=300$ . (ii) For $n=100,200,300,\ldots,1000$ , compute the error at the final time and make a log-log convergence plot, including a reference line for fourth-order convergence.

(a) $u''+ 9u = 9t, \: 0< t< 2\pi, \: u(0) = 1,\: u'(0) = 1; \: \hat{u}(t) = t+\cos (3t)$

(b) $u''+ 9u = \sin(2t), \: 0< t< 2\pi, \: u(0) = 2,\: u'(0) = 1$ ;

$\quad \hat{u}(t) = (1/5) \sin(3t) + 2 \cos (3t)+ (1/5) \sin (2t)$

(c) $u''- 9u = 9t, \: 0< t< 1, \: u(0) = 2,\: u'(0) = -1; \: \hat{u}(t) = e^{3t} + e^{-3t}-t$

(d) $u''+ 4u'+ 4u = t, \: 0< t< 4, \: u(0) = 1,\: u'(0) = 3/4; \: \hat{u}(t) = (3t+5/4)e^{-2t} + (t-1)/4$

(e) $x^2 y'' +5xy' + 4y = 0,\: 1<x<e^2, \: y(1) = 1, \: y'(1) = -1, \: \hat{y}(x) = x^{-2}( 1 + \ln x)$

(f) $2 x^2 y'' +3xy' - y = 0,\: 1<x<16, \: y(1) = 4, \: y'(1) = -1, \: \hat{y}(x) = 2(x^{1/2} + x^{-1})$

(g) $x^2 y'' -xy' + 2y = 0,\: 1<x<e^{\pi}, \: y(1) = 3, \: y'(1) = 4$ ;

$\quad \hat{y}(x) = x \left[ 3 \cos \left( \ln x \right)+\sin \left( \ln x \right) \right]$

(h) $x^2 y'' + 3xy' + 4y = 0,\: e^{\pi/12} < x < e^{\pi}, \: y(e^{\pi/12}) = 0, \: y'(e^{\pi/12}) = -6$ ;

$\quad \hat{y}(x) = x^{-1} \left[ 3 \cos \left( 3 \ln x \right)+\sin \left( 3 \ln x \right) \right]$

Footnotes¶

Americans tend to pronounce these German names as “run-ghuh kut-tah.”
↩