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

d y d t = f ( t , y ) , y ( 0 ) = y 0 .

Calculate K 1 = f ( n h , y n ) then K 2 = f ( ( n + 1 2 ) h , y n + 1 2 h K 1 ) then K 3 = f ( ( n + 1 2 ) h , y n + 1 2 h K 2 ) then K 4 = f ( ( n + 1 ) h , y n + h K 3 ) finally y n + 1 = y n + h 6 K 1 + 2 K 2 + 2 K 3 + K 4

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

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 ( 1 2 h , y 0 + 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 ( 1 2 h , y 0 + 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

y 1 = y 0 + h 6 K 1 + 2 K 2 + 2 K 3 + K 4 = 3 + 0.1 6 0.989992 + 2 × 0.981797 + 2 × 0.981875 0.971390 = 2.901855 to 6 decimal places.
Task!

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

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 ( 1 2 h , y 0 + 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 ( 1 2 h , y 0 + 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

y 1 = y 0 + h 6 K 1 + 2 K 2 + 2 K 3 + K 4 = 0.7 + 0.1 6 0.210000 + 2 × 0.205690 + 2 × 0.205780 + 0.201345 = 0.720571 to 6 d.p.
Exercises
  1. Assuming the notation established earlier, write down the linear multistep scheme corresponding to the choices k = 2 , α 0 = 0 , α 1 = 1 , α 2 = 1 , β 0 = 1 12 , β 1 = 2 3 , β 2 = 5 12 .
  2. A numerical scheme has been used to approximate the solution of

    d y d t = t 2 y 2 y ( 0 ) = 2

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

    y ( 0.3 ) 1.471433 , y ( 0.32 ) 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 ) .

  3. 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.

  4. Find the order of the 2-step linear multistep scheme

    y n + 2 + 2 y n + 1 3 y n = h 10 f n + 2 + 16 f n + 1 + 17 f n

    (Would you recommend using this method?)

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

    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 ) .

  1. y n + 2 y n + 1 = h 12 5 f n + 2 + 8 f n + 1 f n
  2. Evidently the value h = 0.02 will serve our purposes and we seek y 17 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 × 5 f 16 4.6 f 15 = 1.425279

    And we conclude that y ( 0.34 ) 1.425279 , to 6 decimal places.

  3. 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.
  4. 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 = α j = 0 and c 1 = j α j β j = 0

    from which we conclude that the method is consistent.

    We also find that c 2 = 1 2 j 2 α j j β j = 0 c 3 = 1 6 j 3 α j 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).

  5. 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 ( 1 2 h , y 0 + 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 ( 1 2 h , y 0 + 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

    y 1 = y 0 + h 6 K 1 + 2 K 2 + 2 K 3 + K 4 = 2 + 0.4 6 0.250000 + 2 × 0.237954 + 2 × 0.238514 + 0.227753 = 2.095379