s

bi = 1.

i=1

As usual, we say that (11.70) is a method of order p (≥ 1) with respect to

h if „ (h) = O(hp ) as h ’ 0.

As for convergence, since RK methods are one-step methods, consistency

implies stability and, in turn, convergence. As happens for MS methods,

estimates of „ (h) can be derived; however, these estimates are often too

complicated to be pro¬tably used. We only mention that, as for MS meth-

ods, if an RK scheme has a local truncation error „n (h) = O(hp ), for any

n, then also the convergence order will be equal to p.

The following result establishes a relation between order and number of

stages of explicit RK methods.

Property 11.4 The order of an s-stage explicit RK method cannot be

greater than s. Also, there do not exist s-stage explicit RK methods with

order s ≥ 5.

We refer the reader to [But87] for the proofs of this result and the results we

give below. In particular, for orders ranging between 1 and 10, the minimum

11.8 Runge-Kutta (RK) Methods 511

number of stages smin required to get a method of corresponding order is

shown below

order 1 2 3 4 5 6 7 8

smin 1 2 3 4 6 7 9 11

Notice that 4 is the maximum number of stages for which the order of the

method is not less than the number of stages itself. An example of a fourth-

order RK method is provided by the following explicit 4-stage method

h

un+1 = un + (K1 + 2K2 + 2K3 + K4 )

6

K1 = fn ,

(11.73)

K2 = f (tn + h , un + h K1 ),

2 2

K3 = f (tn + h , un + h K2 ),

2 2

K4 = f (tn+1 , un + hK3 ).

As far as implicit schemes are concerned, the maximum achievable order

using s stages is equal to 2s.

Remark 11.4 (The case of systems of ODEs) An RK method can be

readily extended to systems of ODEs. However, the order of an RK method

in the scalar case does not necessarily coincide with that in the vector

case. In particular, for p ≥ 4, a method having order p in the case of the

autonomous system y = f (y), with f : Rm ’ Rn maintains order p even

when applied to an autonomous scalar equation y = f (y), but the converse

is not true. Regarding this concern, see [Lam91], Section 5.8.

11.8.1 Derivation of an Explicit RK Method

The standard technique for deriving an explicit RK method consists of en-

forcing that the highest number of terms in Taylor™s expansion of the exact

solution yn+1 about tn coincide with those of the approximate solution

un+1 , assuming that we take one step of the RK method starting from the

exact solution yn . We provide an example of this technique in the case of

an explicit 2-stage RK method.

Let us consider a 2-stage explicit RK method and assume to dispose at

the n-th step of the exact solution yn . Then

un+1 = yn + hF (tn , yn , h; f ) = yn + h(b1 K1 + b2 K2 ),

K1 = fn , K2 = f (tn + hc2 , yn + hc2 K1 ),

having assumed that (11.72) is satis¬ed. Expanding K2 in a Taylor series

in a neighborhood of tn and truncating the expansion at the second order,

512 11. Numerical Solution of Ordinary Di¬erential Equations

we get

K2 = fn + hc2 (fn,t + K1 fn,y ) + O(h2 ).

We have denoted by fn,z (for z = t or z = y) the partial derivative of f

with respect to z evaluated at (tn , yn ). Then

un+1 = yn + hfn (b1 + b2 ) + h2 c2 b2 (fn,t + fn fn,y ) + O(h3 ).

If we perform the same expansion on the exact solution, we ¬nd

h2 h2

yn + O(h3 ) = yn + hfn + (fn,t + fn fn,y ) + O(h3 ).

yn+1 = yn + hyn +

2 2

Forcing the coe¬cients in the two expansions above to agree, up to higher-

order terms, we obtain that the coe¬cients of the RK method must satisfy

1

b1 + b2 = 1, c2 b2 = 2 .

Thus, there are in¬nitely many 2-stage explicit RK methods with second-

order accuracy. Two examples are the Heun method (11.10) and the modi-

¬ed Euler method (11.91). Of course, with similar (and cumbersome) com-

putations in the case of higher-stage methods, and accounting for a higher

number of terms in Taylor™s expansion, one can generate higher-order RK

methods. For instance, retaining all the terms up to the ¬fth one, we get

scheme (11.73).

11.8.2 Stepsize Adaptivity for RK Methods

Since RK schemes are one-step methods, they are well-suited to adapting

the stepsize h, provided that an e¬cient estimator of the local error is

available. Usually, a tool of this kind is an a posteriori error estimator,

since the a priori local error estimates are too complicated to be used in

practice. The error estimator can be constructed in two ways:

- using the same RK method, but with two di¬erent stepsizes (typically 2h

and h);

- using two RK methods of di¬erent order, but with the same number s of

stages.

In the ¬rst case, if an RK method of order p is being used, one pretends

that, starting from an exact datum un = yn (which would not be available

if n ≥ 1), the local error at tn+1 is less than a ¬xed tolerance. The following

relation holds

yn+1 ’ un+1 = ¦(yn )hp+1 + O(hp+2 ), (11.74)

where ¦ is an unknown function evaluated at yn . (Notice that, in this

special case, yn+1 ’ un+1 = h„n+1 (h)).

11.8 Runge-Kutta (RK) Methods 513

Carrying out the same computation with a stepsize of 2h, starting from

tn’1 , and denoting by un+1 the computed solution, yields

yn+1 ’ un+1 = ¦(yn’1 )(2h)p+1 + O(hp+2 ) = ¦(yn )(2h)p+1 + O(hp+2 )

(11.75)

having expanded also yn’1 with respect to tn . Subtracting (11.74) from

(11.75), we get

(2p+1 ’ 1)hp+1 ¦(yn ) = un+1 ’ un+1 + O(hp+2 ),

from which

un+1 ’ un+1

yn+1 ’ un+1 = E.

(2p+1 ’ 1)

If |E| is less than the ¬xed tolerance µ, the scheme moves to the next time

step, otherwise the estimate is repeated with a halved stepsize. In general,

the stepsize is doubled whenever |E| is less than µ/2p+1 .

This approach yields a considerable increase in the computational e¬ort,

due to the s ’ 1 extra functional evaluations needed to generate the value

un+1 . Moreover, if one needs to half the stepsize, the value un must also

be computed again.

An alternative that does not require extra functional evaluations consists

of using simultaneously two di¬erent RK methods with s stages, of order

p and p + 1, respectively, which share the same set of values Ki . These

methods are synthetically represented by the modi¬ed Butcher array

c A

2

bT

(11.76)

bT 2

2

ET

where the method of order p is identi¬ed by the coe¬cients c, A and b,

while that of order p + 1 is identi¬ed by c, A and b, and where E = b ’ b.

Taking the di¬erence between the approximate solutions at tn+1 pro-

duced by the two methods provides an estimate of the local truncation

error for the scheme of lower order. On the other hand, since the coe¬-

s

cients Ki coincide, this di¬erence is given by h i=1 Ei Ki and thus it does

not require extra functional evaluations.

Notice that, if the solution un+1 computed by the scheme of order p is

used to initialize the scheme at time step n + 2, the method will have order

p, as a whole. If, conversely, the solution computed by the scheme of order

p + 1 is employed, the resulting scheme would still have order p + 1 (exactly

as happens with predictor-corrector methods).

The Runge-Kutta Fehlberg method of fourth-order is one of the most

popular schemes of the form (11.76) and consists of a fourth-order RK

514 11. Numerical Solution of Ordinary Di¬erential Equations

scheme coupled with a ¬fth-order RK method (for this reason, it is known

as the RK45 method). The modi¬ed Butcher array for this method is shown

below

0 0 0 0 0 0 0

1 1

0 0 0 0 0

4 4

3 3 9

0 0 0 0

8 32 32

’ 7200

12 1932 7296

0 0 0

13 2197 2197 2197

’8 ’ 4104

439 3680 845

1 0 0

216 513

’ 27 ’ 3544 ’ 11

1 8 1859

2 0

2 2565 4104 40

’1

25 1408 2197

0 0

216 2565 4104 5

’ 50

16 6656 28561 9 2

0

135 12825 56430 55

’ 4275 ’ 75240

1 128 2197 1 2

0

360 50 55

This method tends to underestimate the error. As such, its use is not com-

pletely reliable when the stepsize h is large.

Remark 11.5 MATLAB provides a package tool funfun, which, besides

the two classical Runge-Kutta Fehlberg methods, RK23 (second-order and

third-order pair) and RK45 (fourth-order and ¬fth-order pair), also imple-

ments other methods suitable for solving sti¬ problems. These methods are

derived from BDF methods (see [SR97]) and are included in the MATLAB

program ode15s.

11.8.3 Implicit RK Methods

Implicit RK methods can be derived from the integral formulation of the

Cauchy problem (11.2). In fact, if a quadrature formula with s nodes in

(tn , tn+1 ) is employed to approximate the integral of f (which we assume,

for simplicity, to depend only on t), we get

tn+1

s

f („ ) d„ h bj f (tn + cj h)

j=1

tn

having denoted by bj the weights and by tn + cj h the quadrature nodes. It

can be proved (see [But64]) that for any RK formula (11.70)-(11.71), there

exists a correspondence between the coe¬cients bj , cj of the formula and

the weights and nodes of a Gauss quadrature rule.

In particular, the coe¬cients c1 , . . . , cs are the roots of the Legendre

polynomial Ls in the variable x = 2c ’ 1, so that x ∈ [’1, 1]. Once the

11.8 Runge-Kutta (RK) Methods 515

s coe¬cients cj have been found, we can construct RK methods of order

2s, by determining the coe¬cients aij and bj as being the solutions of the

linear systems

s

ck’1 aij = (1/k)ck , k = 1, 2, . . . , s,

i

j

j=1

s

ck’1 bj = 1/k, k = 1, 2, . . . , s.

j

j=1

The following families can be derived:

1. Gauss-Legendre RK methods, if Gauss-Legendre quadrature nodes are

used. These methods, for a ¬xed number of stages s, attain the maximum

possible order 2s. Remarkable examples are the one-stage method (implicit

midpoint method) of order 2

1 1

2 2

1 1

un+1 = un + hf tn + 2 h, 2 (un + un+1 ) ,

1

and the 2-stage method of order 4, described by the following Butcher array

√ √

3’ 3 3’2 3

1

6 4 12

√ √

3+ 3 3+2 3 1

6 12 4

1 1

2 2

2. Gauss-Radau methods, which are characterized by the fact that the

quadrature nodes include one of the two endpoints of the interval (tn , tn+1 ).

The maximum order that can be achieved by these methods is 2s ’ 1, when

s stages are used. Elementary examples correspond to the following Butcher

arrays

’ 12

1 5 1

3 12

3 1

1

0 1 11 4 4

, ,

1 1

3 1

4 4

and have order 1, 1 and 3, respectively. The Butcher array in the middle

represents the backward Euler method.

3. Gauss-Lobatto methods, where both the endpoints tn and tn+1 are quadra-

ture nodes. The maximum order that can be achieved using s stages is

2s ’ 2. We recall the methods of the family corresponding to the following

516 11. Numerical Solution of Ordinary Di¬erential Equations

Butcher arrays

’1

1 1

0 6 3 6

1

0 0 0 0 0

’ 12

2 1 1 5 1

1 1 1

1 1 0 2 6 12

2 2 2

, , 1 2 1

1 6 3 6

1 1 1 1

2 2 2 2 1 2 1

6 3 6

which have order 2, 2 and 3, respectively. The ¬rst array represents the

Crank-Nicolson method.

As for semi-implicit RK methods, we limit ourselves to mentioning the

case of DIRK methods (diagonally implicit RK), which, for s = 3, are

represented by the following Butcher array

1+µ 1+µ

0 0

2 2

’µ 1+µ

1

0

2 2 2

1 + µ ’1 ’ 2µ

1’µ 1+µ

2 2

1’

1 1 1

6µ2 3µ2 6µ2

The parameter µ represents one of the three roots of 3µ3 ’ 3µ ’ 1 = 0 (i.e.,

√ √ √

(2/ 3) cos(10o ), ’(2/ 3) cos(50o ), ’(2/ 3) cos(70o )). The maximum or-

der that has been determined in the literature for these methods is 4.

11.8.4 Regions of Absolute Stability for RK Methods

Applying an s-stage RK method to the model problem (11.24) yields

s s

Ki = un + h» aij Kj , un+1 = un + h» bi Ki , (11.77)

i=1 i=1

that is, a ¬rst-order di¬erence equation. If K and 1 are the vectors of com-

ponents (K1 , . . . , Ks )T and (1, . . . , 1)T , respectively, then (11.77) becomes

un+1 = un + h»bT K,

K = un 1 + h»AK,

from which, K = (I ’ h»A)’1 1un and thus

un+1 = 1 + h»bT (I ’ h»A)’1 1 un = R(h»)un

where R(h») is the so-called stability function.

The RK method is absolutely stable, i.e., the sequence {un } satis¬es

(11.25), i¬ |R(h»)| < 1. Its region of absolute stability is given by

A = {z = h» ∈ C such that |R(h»)| < 1} .

11.9 Systems of ODEs 517

If the method is explicit, A is strictly lower triangular and the function R

can be written in the following form (see [DV84])

det(I ’ h»A + h»1bT )

R(h») = .

det(I ’ h»A)

Thus since det(I’h»A) = 1, R(h») is a polynomial function in the variable

h», |R(h»)| can never be less than 1 for all values of h». Consequently, A

can never be unbounded for an explicit RK method.

In the special case of an explicit RK of order s = 1, . . . , 4, one gets (see

[Lam91])

s

1

(h»)k .

R(h») =

k!

k=0

The corresponding regions of absolute stability are drawn in Figure 11.10.

Notice that, unlike MS methods, the regions of absolute stability of RK

methods increase in size as the order grows.

4

3.5

3 s=4

s=3

2.5

2

s=2

1.5

1