Stability of the fourth order Runge-Kutta method for the solution
of systems of differential equations
By Abbas I. Abdel Karim*
This paper tackles the problem of the region of stability of the fourth order Runge-Kutta method
for the solution of systems of differential equations. Such a region can be characterized by means
of linear transformation but can not be given in a closed form. In the present case, the region
has been determined using the electronic digital computer Z22.
A system of m ordinary differential equations can be The corresponding exact solution is:
Downloaded from https://academic.oup.com/comjnl/article/9/3/308/406366 by guest on 25 December 2021
written in vector form as follows:
Y' = dY/dx = F(x; Y), (1)
Y(xm + h)= Y{xn) X (4)
with the initial value where
Y(a) = Y0(a<x<b) 7) = F{xn + 8jh; Y(xn) + - 1)).
where Subtracting (3) from (4) we have
~f\ (x;yuy2, . . .,yj' Yn+l- Y(xn+h)
h(x;yuy2, . . .,ym) = Yn- Y(xn) + -A YJ(KJ - AC/)) + R (5)
Y = , F(x; Y) =
where R is the difference between the truncation and
the round-off error.
fm(x; yu y2, . . .,ym) Let the error en of the solution at the point xn be
denned by
Let:
Y(x) be the exact solution of (1), £n = Yn - Y{xn). (6)
Y* the exact solution of the difference equation Hence we obtain
at x = xn,
Yn the practical (computed) solution at x — xn. K{ - AC1) = F(xn; Yn) - F(xn; Yn - en). (7)
In the fourth order Runge-Kutta method (1901) the Assume that
vector Y*+1 is then given by: (a) F(xn; Y) is a continuous function of Y contained
in the closed interval whose end-points are Y(xn) and
(2) Yn.
i(xn; Y)
where (b) The matrix ( ) for i j = 1,2,. .., m
K; = F(xn + 8jh; Y: + SjhKjL,),
exists for Y in the open interval whose end-points are
1 0 y(jcn) and Yh. Therefore, by the mean-value theorem
2 we get
((Yj),
2 i
1 1 A, - AO) = ( | ^ K 1 , 7 = 1 , 2 , . . . , in. (8)
and
id hh =
= the
the step
step of integration.
of integration. ,~df. \
The computed (practical) solution at x = xn+l is as The arguments of the matrix ( ^ J are(x n ; Yn-Qen),
follows:
where 0 is a diagonal matrix = diag (#,) and O<0,<1,
Y Y
(3)
i— 1,2, . . ., m.
n+ 1 = n YjKj
Similarly we get
where
Kj = F(xn + 8jh; Yn
and Rn is the round-off error. where / is the unit matrix,
* Helupolice, El Sheek Mohamed Refaat St. 5, Cairo, Egypt, U.A.R.
308
Differential equations
K" = T-lJ"T. (16)
From (15) it follows that
(17)
i+T^ En- The corresponding homogeneous system of (17) has
the form
Substituting (8), (9), (10) and (11) in (5) we get the
difference-equation for the error (18)
Case I. If the eigenvalues pt of the matrix / are all
Downloaded from https://academic.oup.com/comjnl/article/9/3/308/406366 by guest on 25 December 2021
distinct, then the matrix K has these eigenvalues in its
principal diagonal and zeros elsewhere. In this case the
system (18) has the form
sn. (19)
The Jacobian matrices ( =T-0 = / are written in (12)
individually to denote that the argument of the matrices 0
are in general different (see Carr, 1958). v=0
In studying the local growth of the error en it is
The general solution of this system is clearly of the
reasonable to assume that R and the matrices (.r-'- ) in form
a small interval (xn, xn+x) are both constants, and they
may be written in the form R and J, respectively (see
Hamming, 1962). In practice they vary slowly from where the ct are constants. Here e/>n denotes the /th
step to step, otherwise the step size of integration is too element of the vector en.
large.
Under these assumptions we get from (12) the non- Case II. If the eigenvalues p, of the matrix / are
homogeneous system of linear difference-equations with not all distinct, then the matrix / may be transformed to
constant real coefficients in the form
0
(13) Kl2
\
To solve the system (13) we introduce a new vector Kle.
variable e in place of the variable e by means of a
non-singular linear transformation
\
e = r - 1 e , where det T # 0, (14) K,
so that the system (13) is transformed into 0 \J
4 W where each submatrix Ku of order j is of the form
(15)
'P. 1 0-
l
where T~ R = A. Pi 1
P; •
In particular the transformation T may be chosen so Ku = •
as to put the matrix T~XJT into the classical canonical
form (see Zurmuhl, 1961). . 1
With .0 Pi.
= T~XJT with pi in each position in the leading diagonal, unity
we get in each position in the diagonal immediately above it,
309
Differential equations
and with zeros elsewhere. In this case of non-distinct With the assumption that
eigenvalues, we have for every submatrix of order y a
system of equations in the form (25)
we
(21)
where •-('-!.£*) (26)
In conclusion it can be said that the general solution
P* of the non-homogeneous difference equation (13) can be
Pn 1
written as follows whether the eigenvalues pf are distinct
or not:
Downloaded from https://academic.oup.com/comjnl/article/9/3/308/406366 by guest on 25 December 2021
*
0
It can be proved that
,=o
(see Schmeidler, 1949).
Thus the system (21) has the form
Av v - v , . . . u
v= 0 l!v-o
h
- o
v 0
4! • • • ° (22)
We can write (22) in the form Where c,- are constants for distinct eigenvalues or also
polynomials in n for repeated eigenvalues.
fi/.n+i — 2 -7; 2 - £,• + ,-, „, (i = 1, 2, . . ., y ;
y=0 J • v=0 Definition
0 < j < 4). (23) The numerical integration method (3) is said to be
stable according to Wilf's criterion (1959), when
It can be proved that the general solution of the system
(23) can be written in the form 2 —f- (for all i = 1, 2, . . ., m) are inside the unit
v=0 VI
circle in the complex plane, i.e. when
C
Ht-D.n = t-i
< 1, i = l , 2 , . . ., m. (28)
where cY_,-, are polynomials in n of the fth degree;
(1 = 0, 1 . . ., y - 1). When this definition of stability is satisfied, it is clear
For the particular solution of the non-homogeneous from (27) that the error en will be dominated by the
difference equation (17) we set
en+, = "en — e = constant vector ¥= zero vector. It is clear that det ( S -V jv _ / ) = 0, when S = 1
Hence \v=0l > ! J v>=0 v\
at least for one /, which is the case where one of hpi lies on the
boundary of the region R (see Fig. 1). Furthermore, the particular
= R. solution of the non-homogeneous difference equation (13) is a
polynomial in n.
310
Differential equations
hv~ V - l
2 - 7 ) . R when n is large.
v=l v\ /
From the condition (28) it can be easily proved that the
fourth order Runge-Kutta method for the solution of
systems of ordinary differential equations is stable, if
and only if
(1) for real Pi: i{-^[4V(2349) - 172] -
- -^[4V(2349) + 172] - 4} < hPl < 0
i.e. -2-785 . . . < hPi < 0,
(2) for pure imaginary p, : 0 < \hp,\ <
Downloaded from https://academic.oup.com/comjnl/article/9/3/308/406366 by guest on 25 December 2021
(3) for complex pt : hp{ lie inside the region /?,f (see
Fig. 1).
Special case
In the case of an individual differential equation
y' =f(x; y), y(a) = y0, a < x < b, we have our final
result as follows: The fourth order Runge-Kutta method
is stable, if -2-785 . . . < h ^ < 0.
7>y
This is a well known result which is given by Liniger
(1957).
• The exact values in (1) and (2) can be obtained by solving the
inequality (28).
t The region R is estimated by means of the electronic digital
Fig. 1 computer Z22.
References
CARR, J. W., Ill (1958). "Error bounds for the Runge-Kutta single-step integration process", /. Assoc. Comp. Mach. Vol. 5,
pp. 39-44.
HAMMING, R. W. (1962). Numerical Methods for Scientists and Engineers, Part II, Chap. 14, New York: McGraw-Hill Book Co.
KUTTA, W. (1901). "Beitrag zur naherungsweisen Integration totaler Differentialgleichungen", Zeitschr. fur Math. u. Physik,
Vol. 46, pp. 435^53.
LINIGER, W. (1957). "Zur Stabilitat der numerischen Integrationsmethoden fur Differentialgleichungen." Doctoral thesis U.
of Lausanne, Zurich.
SCHMEIDLER, W. (1949). Vortrdge iiber Determinanten und Matrizen, Berlin: Akademie-Verlag Berlin.
WILF, S. (1959). "A stability criterion for numerical integration", /. Assoc. Comp. Mach. Vol. 6, Nr. 3, pp. 363-365.
ZURMUHL, R. (1961). Matrizen, Berlin: Springer-Verlag.
311
