An example of a Runge-Kutta method

2 An example of a Runge-Kutta method

A full discussion of the so-called Runge-Kutta methods is not required here, but we do need to touch on them to resolve a remaining issue in the implementation of linear multistep schemes.

The problem with linear multistep methods is that a zero-stable, 1-step method can never be better than second order (you need not worry about why this is true, it was proved in the latter half of the last century by a man called Dahlquist). We have seen methods of higher order than 2, but they were all at least 2-step methods. And the problem with 2-step methods is that we need 2 starting values to implement them and we are only ever given 1 starting value: the initial condition $y (0)$ .

One way out of this “Catch 22" is to use a Runge-Kutta method to generate the extra starting value(s) we need. Runge-Kutta methods are not linear multistep methods and do not suffer from the problem mentioned above. There is no such thing as a free lunch, of course, and Runge-Kutta methods are generally more expensive in effort to implement than linear multistep methods because of the number of evaluations of $f$ required at each time step.

The following Key Point gives a statement of what is, perhaps, the most popular Runge-Kutta method (sometimes called “RK4").

Key Point 14

Runge Kutta method (RK4)

Consider the usual initial value problem

$\frac{d y}{d t} = f (t, y), y (0) = y_{0} .$

\begin{array}{rcl} Calculate K_{1} & = & f (n h, y_{n}) \\ then K_{2} & = & f ((n + \frac{1}{2}) h, y_{n} + \frac{1}{2} h K_{1}) \\ then K_{3} & = & f ((n + \frac{1}{2}) h, y_{n} + \frac{1}{2} h K_{2}) \\ then K_{4} & = & f ((n + 1) h, y_{n} + h K_{3}) \\ finally y_{n + 1} & = & y_{n} + \frac{h}{6} (K_{1} + 2 K_{2} + 2 K_{3} + K_{4}) \end{array}

Notice that each calculation is explicit, all of the right-hand sides in the formulae in the Key Point above involve known quantities.

Example 11

Suppose that $y = y (t)$ is the solution to the initial value problem

$\frac{d y}{d t} = cos (y) y (0) = 3$

Carry out one time step of the Runge-Kutta method RK4 with a step size of $h = 0.1$ so as to obtain an approximation to $y (0.1)$ .

Solution

The iteration must be carried out in four stages. We start by calculating

$K_{1} = f (0, y_{0}) = f (0, 3) = - 0.989992$

a value we now use in finding

$K_{2} = f (\frac{1}{2} h, y_{0} + \frac{1}{2} h K_{1}) = f (0.05, 2.950500) = - 0.981797$

This value $K_{2}$ is now used in our evaluation of

$K_{3} = f (\frac{1}{2} h, y_{0} + \frac{1}{2} h K_{2}) = f (0.05, 2.950910) = - 0.981875$

which, in turn, is used in

$K_{4} = f (h, y_{0} + h K_{3}) = f (0.1, 2.901812) = - 0.971390$

All four of these values are then used to complete the iteration

\begin{array}{rcl} y_{1} & = & y_{0} + \frac{h}{6} (K_{1} + 2 K_{2} + 2 K_{3} + K_{4}) \\ = & 3 + \frac{0.1}{6} (- 0.989992 + 2 \times - 0.981797 + 2 \times - 0.981875 - 0.971390) \\ = & 2.901855 to 6 decimal places. \end{array}

Task!

Suppose that $y = y (t)$ is the solution to the initial value problem

$\frac{d y}{d t} = y (1 - y) y (0) = 0.7$

Carry out one time step of the Runge-Kutta method RK4 with a step size of $h = 0.1$ so as to obtain an approximation to $y (0.1)$ .

The time step must be carried out in four stages. We start by calculating

$K_{1} = f (0, y_{0}) = f (0, 0.7) = 0.210000$

a value we now use in finding

$K_{2} = f (\frac{1}{2} h, y_{0} + \frac{1}{2} h K_{1}) = f (0.05, 0.710500) = 0.205690$

This value $K_{2}$ is now used in our evaluation of

$K_{3} = f (\frac{1}{2} h, y_{0} + \frac{1}{2} h K_{2}) = f (0.05, 0.710284) = 0.205780$

which, in turn, is used in

$K_{4} = f (h, y_{0} + h K_{3}) = f (0.1, 0.720578) = 0.201345$

All four of these values are then used to complete the time step

\begin{array}{rcl} y_{1} & = & y_{0} + \frac{h}{6} (K_{1} + 2 K_{2} + 2 K_{3} + K_{4}) \\ = & 0.7 + \frac{0.1}{6} (0.210000 + 2 \times 0.205690 + 2 \times 0.205780 + 0.201345) \\ = & 0.720571 to 6 d.p. \end{array}

Exercises

Assuming the notation established earlier, write down the linear multistep scheme corresponding to the choices $k = 2$ , $α_{0} = 0$ , $α_{1} = - 1$ , $α_{2} = 1$ , $β_{0} = \frac{- 1}{12}$ , $β_{1} = \frac{2}{3}$ , $β_{2} = \frac{5}{12}$ .
A numerical scheme has been used to approximate the solution of
$\frac{d y}{d t} = t^{2} - y^{2} y (0) = 2$

and has given the following estimates, to 6 decimal places,

$y (0.3) \approx 1.471433, y (0.32) \approx 1.447892$

Now use the 2-step, explicit linear multistep scheme

$y_{n + 2} - 1.6 y_{n + 1} + 0.6 y_{n} = h (5 f_{n + 1} - 4.6 f_{n})$

to approximate $y (0.34)$ .
Find the roots of the first characteristic polynomial for the linear multistep scheme
$5 y_{n + 2} + 3 y_{n + 1} - 2 y_{n} = h (f_{n + 2} + 2 f_{n + 1} + f_{n})$

and hence determine whether or not the scheme is zero stable.
Find the order of the 2-step linear multistep scheme
$y_{n + 2} + 2 y_{n + 1} - 3 y_{n} = \frac{h}{10} (f_{n + 2} + 16 f_{n + 1} + 17 f_{n})$

(Would you recommend using this method?)
Suppose that $y = y (t)$ is the solution to the initial value problem
$\frac{d y}{d t} = 1 ∕ y^{2} y (0) = 2$

Carry out one time step of the Runge-Kutta method RK4 with a step size of $h = 0.4$ so as to obtain an approximation to $y (0.4)$ .

$y_{n + 2} - y_{n + 1} = \frac{h}{12} (5 f_{n + 2} + 8 f_{n + 1} - f_{n})$
Evidently the value $h = 0.02$ will serve our purposes and we seek $y_{17} \approx y (0.34)$ . The values we will need to use in our implementation of the 2-step scheme are $y_{16} = 1.447892$ , $y_{15} = 1.471433$ and $f_{16} = f (0.32, y_{16}) = - 1.993991 f_{15} = f (0.3, y_{15}) = - 2.075116$ since $f (t, y) = t^{2} - y^{2}$ . It follows that
$y_{17} = 1.6 y_{16} - 0.6 y_{15} + 0.02 \times (5 f_{16} - 4.6 f_{15}) = 1.425279$

And we conclude that $y (0.34) \approx 1.425279$ , to 6 decimal places.
The first characteristic polynomial is $ρ (z) = α_{2} z^{2} + α_{1} z + α_{0} = 5 z^{2} + 3 z - 2$ and the roots of $ρ (z) = 0$ can be found from the quadratic formula. In this case the roots are real and distinct and are equal to $0.4 and - 1.$ In the case of roots that are distinct zero-stability requires that the absolute values have magnitude less than or equal to 1 . Consequently we conclude that the method is zero stable.
In the established notation we have $α_{2} = 1$ , $α_{1} = 2$ and $α_{0} = - 3$ . The $b e t a$ terms similarly come from the coefficients on the right hand side (remembering the denominator of 10).
Now $c_{0} = \sum α_{j} = 0 and c_{1} = \sum j α_{j} - β_{j} = 0$

from which we conclude that the method is consistent.

We also find that $c_{2} = \sum \frac{1}{2} j^{2} α_{j} - j β_{j} = 0 c_{3} = \sum \frac{1}{6} j^{3} α_{j} - \frac{1}{2} j^{2} β_{j} = - 0.533333$

so that the method is of order 2 . This method is not to be recommended however (check the zero stability).
Each time step must be carried out in four stages. We start by calculating
$K_{1} = f (0, y_{0}) = f (0, 2) = 0.250000$

a value we now use in finding $K_{2} = f (\frac{1}{2} h, y_{0} + \frac{1}{2} h K_{1}) = f (0.2, 2.050000) = 0.237954$

This value $K_{2}$ is now used in our evaluation of

$K_{3} = f (\frac{1}{2} h, y_{0} + \frac{1}{2} h K_{2}) = f (0.2, 2.047591) = 0.238514$

which, in turn, is used in $K_{4} = f (h, y_{0} + h K_{3}) = f (0.4, 2.095406) = 0.227753$

All four of these values are then used to complete the time step
$\begin{array}{rcl} y_{1} & = & y_{0} + \frac{h}{6} (K_{1} + 2 K_{2} + 2 K_{3} + K_{4}) \\ = & 2 + \frac{0.4}{6} (0.250000 + 2 \times 0.237954 + 2 \times 0.238514 + 0.227753) = 2.095379 \end{array}$