<< . .

. 74
( : 95)



. . >>

(see [Lam91]) that this happens i¬
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

<< . .

. 74
( : 95)



. . >>