This item is only available as the following downloads:
k(-4 f M -) q3
NATIONAL ADVISORY COMMITTEE FOR AERONAUTICS
TECHNICAL MEMORANDUM 1405
ON THE INSTABILITY OF METHODS FOR THE INTEGRATION
OF ORDINARY DIFFERENTIAL EQUATIONS1
By Heinz Rutishauser2
In spite of the remarkable publication of J. Todd (ref. 1) the
essential points of which are related below, the author has since
observed several times methods for the numerical integration of differ-
ential equations which, although subject only to a temptingly small
truncation error, nevertheless involve the great danger of numerical
instability. I want to state beforehand that this danger hardly exists
for the well-tested methods of Runge-Kutta and Adams (extrapolation
methods) if they are applied correctly.
It is a natural characteristic that a differential equation to be
solved numerically is approximated by a difference equation, and that
the latter is then solved. In order not to be forced to select an all
too small interval, one prefers difference equations which approximate
the differential equation as closely as possible but in compensation
are of higher order than the original differential equation. Precisely
in this, however, there lies a danger because the difference equation
thereby has a greater diversity of possible solutions, and it may well
happen that the numerical integration yields precisely one of the
extraneous solutions which only at the beginning is in any way related
to the desired solution of the differential equation. In the paper of
J. Todd mentioned before several examples of this type have been
Consideration of the pertinent variation equation is particularly
informative. It is very well possible that the differential variation
equation is stable, that is, that it contains only converging solutions,
whereas the difference variation equation is unstable since it possesses,
due to the increased diversity of solutions, aside from the converging
solutions, also solutions which increase exponentially. A deviation
from the correct solution, once it exists, small as it may be and such
1"iber die Instabilitat von Methoden zur Integration gewohnlicher
Differentialgleichungen." ZAMP, Kurze Mitteilungen, vol. III, 1952,
2Institut fur angewandte Mathematik der ETH, Zurich.
deviations are unavoidable, because of the rounding-off errors there-
fore increases rapidly and may finally falsify considerably the solution
obtained. Yet we want to emphasize this once more this instability
is caused only by an inappropriate integration method.
In the following discussion, several customary methods are examined
from this viewpoint, and simple criteria for the stability of such meth-
ods are indicated. For the rest, this report does not deal with error
DIFFERENTIAL EQUATIONS OF THE FIRST ORDER
y' = f(x,y)
(a) Integration by Means of Simpson's Rule5
yk+i = Yk-i -3 k+i + k k-) (14
This relationship, together with the differential equation, yields
two equations for the unknown quantities yk+1 and y' which are
solved mostly by iteration. If the differential equation is linear or
quadratic in y, the iteration can be avoided. We assume, however, that
one passes over, in any case, to the next integration interval only when
the relationship (1) is satisfied.
The difference variation equation pertaining to (1) is evidently
5The phenomenon was observed on this example, and correctly inter-
preted, also by Mr. G. Dahlquist, Stockholm (Lecture at the GaMM-
convention 1951 at Freiburg im Breisgau). Compare also: Z. angew.
Math. Mech., 51, 259, 1951.
4Where h signifies the length of the integration interval, and
Yk stands as an abbreviation for y(kh).
flk+l( yk+) y k 3 1 + f ,k- k-1 =
If one assumes f to be constant, and chooses the expression nk = ,k
for the solution of this equation, one obtains for A a quadratic
equation with the solutions
A = 1 + hf + f2 + e y
'2 -1 y T8 f +18 h -e
One recognizes easily that, of the two fundamental solutions 11,k = 1k
and 12,k 2k of the difference variation equation, the first one
approximates the solution of the differential variation equation whereas
the second one is brought in by the numerical method.
A2 ~ (-1) e
represents for fy <0, thus precisely when the differential equation is
stable, an oscillation which is slowly exponentially increasing. This
has the effect that a small disturbance of the numerical solution -
caused by a rounding-off or a truncation error is intensified in the
further course of the integration and finally gets completely out of
hand. In Collatz' (ref. 2) book, the phenomenon is denoted as "roughening
phenomenon"; means for elimination of this inconvenience are given. On
the other hand, the explanation given there is not complete; the phe-
nomenon occurs only for fy < 0; there is nothing to be apprehensive of
for fy 0 which is very important particularly in regard to the ordi-
nary Simpson's integration rule (f3 0o.
(b) Integration According to Runge-Kutta and Similar Methods
Since these methods calculate yk+l from yk according to a
prescribed rule and without use of the preceding values yk-l' Yk-2' *f
the order remains unchanged in the transition from the differential
equation to the difference equation; thus no foreign solutions are brought
in, and no instability is to be feared.
The same property can be found in a method indicated by W. E. Milne
(c) Integration According to Adams
We consider a four-point formula
k+l = k + rk+l + 1y'k 5'Yk- + y'k-2) + (2)
This again yields, together with the differential equation, two equa-
tions for the unknown quantities yk+l and y'k+1 which are mostly
solved by iteration.
The difference variation equation pertaining to (2) becomes
( fyk+L)%+ L+24 y,k) 'k + 24 fy,k-lk-l y,k-2nk-2 = 0
If one again considers fy as constant, the expression nk = xk yields
an equation of the third order for 7. A solution Al of this equation
lies very close to e therefore 1, k = A1 corresponds to the
solution of the differential variation equation whereas A2k and A k
are extraneous solutions.
However, the equation for A is reduced to A3 A2 = 0 when h
tends toward 0 so that, for a sufficiently small h, one will have at
any rate small A2 and A3, namely \-hfy/24. The extraneous solu-
tions 2,k = 2k and ,k = k thus converge rapidly. For a suffi-
ciently small h, Adams' method is therefore stable.
(d) Variants of (c)
In order to improve the accuracy of Adams' method, one may use also
other expressions for the corrector instead of (2). As long as the
TM 1405 5
corresponding five- or six-point formulas are involved, there are no
objections, but one has to be careful when yk+l is not calculated
from yk and.the derivatives as in (2) but perhaps from yk-, or
Yk-3 and the derivatives, as for instance in
Yk+1 = Yk-5 + 'k + 'k- + 52y'k2 1k + h7 .
In fact the pertinent difference variation equation has a solution
7? X ^k with = l h i + .
so that the method is unstable for fy < 0.
DIFFERENTIAL EQUATIONS OF THE SECOND ORDER
y" = f(x,y,y')
Insofar as these equations are solved by separation into a system
of two equations of the first order, what was said so far is valid.
Particularly in the case of numerical integration of damped oscillations
we must caution against the methods (a) and (d).
However, there exist also methods which solve an equation of the
second order without transformation into a system:
(e) The Method of Central Differences5
The formulas on which this method is based are (especially for
k+1 = 2yk Yk-1 + j"k+l + 10"k + k- 4)
y' y'k + k.l + + Y" k-) (5)
5Compare reference 2, p. 80.
They yield, together with the differential equation, three equations
for the unknown quantities yk+l,, 'k+1' and y"+ The two simul-
taneous difference variation equations pertinent to (4) and (5) are
solved with the expression k = pk, 'k = Q under the assumption
of a constant fy and f r; because of
Ik+1 = fk+i + y'k+
one obtains, with the abbreviations a
hfyi/5, the equations
p [(2 2A + 1) a(X2 + 10 + 1)| -
and b for
q S1 b(A2 + 10A + 1) = 0
q (2 1) b(V2 + 4A + 11- p p a(A2
+ 4A + 1) = 0
These equations can exist simultaneously with (p,q)
the determinant of this equation system for p and
obtains after some calculations
j (0,0) only when
q vanishes; one
(A 1) (A2 1)(A 1) a(A + 1)(A2 + 10A + 1) -
b( 1)(A2 + 4A + 1) = 0
The four solutions of these equations are
Al = 1+ a]h + .
\2 = 1 + a2h +
= = fyi +
where a,1 and a2 are the solutions of
the equation a2 a-fyr fy = 0
Evidently ?1k and x2k are the regular solutions of the difference
variation equation; they correspond to two fundamental solutions of
the differential variation equation; X?3 and A4, in contrast, are
extraneous. As long as f1r I 0, there is nothing to fear, in partic-
uiar, tiie method may be strongly recommended for a y'-free equation,
but for fyl < 0 (damped oscillations)
'4 k (-1)ke-(khf)) ,
TI4 ,k= k (ljhe(kFI.5)f
increases, and ?i, too, may still become dangerous because 3,k 1
also becomes finally very large, compared to a function converging to-
The author completely calculated the example y" + y' + 1.25y = 0
with the initial conditions y. = 0, y' = 1 (exact solution:
e sin x on the sequence-controlled computing machine of the ETH.
There follow a few excerpts from the thus obtained table of func-
tions (we calculated with h = 0.1):
In this region nothing conspicuous is noticeable yet, the y-values
deviate from the true values approximately by one in the last decimal
place, and only formation of the differences for the y'-values reveals
a certain irregularity. For t = 17, however, the influence of Ask
becomes pronounced for the y'-values and also for the differences of
x y y'
4.8 -0.0905699 0.0551227
4.9 .0847792 .0584842
5.0 .0787152 .0626410
5.1 .0722891 .0656575
5.2 .0656175 .0676070
The considerably weaker oscillation of the y-values follows also
from the equations (6): for A = X4 one obtains from the first of
P~ b 8 h f ,,, here therefore p ~ q
q 4 4 6 600
The further course of numerical integration does not require any comment:
The author is well aware that the assumption of a constant fy
in the above considerations greatly restricts the generality.
ever, the results show that what matters is only the sign of these quan-
tities, and this sign is indeed invariable in a great many cases. The
statements are, therefore, qualitatively almost generally valid. Only
when fy changes its sign, from time to time in the course of the inte-
3rstion, for instance when method (a) is being used, a special case
arises since the occurring oscillations alternately increase and are
x y y'
17.0 -0.00019574 0.00005017
17.1 .00019061 .00005253
17.2 .00018566 .00008620
17.5 .00017524 .00008259
17.4 .00016548 .000112355
17 .' .00015475 .00010258
x y y'
22.8 -0.00000815 0.00005520
22.9 .00000864 -.00006078
25.0 .00000862 .00005968
23.1 .00000887 -.00006247
25.2 .00000861 .00006601
23.3 .00000868 -.00006486
29.5 .00000140 -.00055057
29.6 .00000041 .00054889
29.7 .00000144 -.00056695
29.8 .00000050 .00058682
29.9 .00000148 .0oo6060o
50.0 .00000060 -.ooo62755
TM 1403 75
Consideration of the examples suggests the conjecture that insta-
bility may occur precisely in the case of integration methods which
form Yk+1 by integration of y' over several intervals (two intervals
in method (a), one interval in (b) and (c), four intervals in (d), two
intervals in (e)). However, this is not exactly the case and we shall
therefore subject a general integration method to an examination.
Almost all known methods use relationships which are
the general formula
yk+ = i_ 4j Yk+j
a(0) yk+j +
ajj Yvk+j + +
N (0) (N)
h alj Yk+j
y k+1 2
N-1 1 (1) (N)
Y: Ni k+j
n-1 + a(n-l) (n-1)
k+l -m n-l,j k+j
N-n+11 (n-1) (N)
+ .. +
There are, in addition, differential equations
Fn (k+lyk+' k+ ) =
F (+1(xk k (n+l) =0
* y( = 0
which form, together with the equations (7), N + 1 equations for the
N + 1 unknowns yk+', Yk+1' k+1 Generally, N = n; however,
W. E. Milne (ref. 5) uses in the method previously mentioned higher
derivatives than those appearing in the differential equation. There-
fore he differentiates them several times in order to obtain the required
number of relationships. One may thus obtain the equations Fn+1
to FN by differentiating the initial equation Fn.
If one, furthermore, brings everything in the equations (7) to one
side, the variation equations for the entire system (7) and (8) read
with ih -a -
P j>t h ak = O
0- Fi (p
If one uses for this the expression q.
(1 = 0, 1, ., n 1)
(i = n, n + 1, ., N)
= p A, one obtains in sub-
stituting a system of N + 1 simultaneous equations for the p which
can be satisfied only when the determinant of the system disappears
r-)- 7h-ip = 0o
--- P, = 0
(1 = 1, 1, ., n -1)
(i = n, n + 1, ., N)
If one defines, in addition, with the coefficients a
the formulas (7) the functions Pj
A () = a (i)j+m
wherein A. .
0 for p < i, the characteristic determinant reads as
0 0 'An-1,n-
Fn,yI Fn,y(n) 0
. . FN,y(N)
. . .
If the method is to reproduce exactly a polynomial of the ith degree,
which is the solution of a differential equation, together with its
derivatives and this one may require the conditions
(i) (i+) y(i+2) y(i+5) y(N) s o
J J J J j
must be compatible. Hence follows, however, by substitution into the
ith of the equations (7): a(i) = 0, therefore Aii(l) = 0.
The equation D(A) = 0 which is decisive for the stability of the
method must have n solutions in the neighborhood of ?A = 1, corre-
sponding to the n independent solutions of the variation equation.
In fact one finds for h = O0 where D(A) is, except for one factor,
reduced to A00All An-1, n-l, that 'A = 1 is an n-fold zero of
D(o), because of Aii(l) = 0.
All other zeros of D(?) correspond therefore to extraneous solu-
tions of the difference equation; in order to make the method stable,
they must lie, for a sufficiently small h, in the interior or at most
on the periphery of the unit circle. This is certainly the case when
for h = 0 all zeros of D(A) lie in the interior of the unit circle,
and certainly not when individual ones are outside it. Therefore:
Sufficient condition for the stability of the method (7) for
sufficiently small h:
Aii(D) = (i)ij+m
possess, aside from the trivial simple zero A = 1 only zeros with
h < i.
necessary condition: None of the functions Aii(.) has a zero
Outside the unit circle.
0J. Todd considered methods which satisfy not even the necessary
condition. The solution obtained then becomes completely useless
already after a few intervals.
The remaining functions Al(AX) can influence the stability only
when the necessary condition, but not the sufficient condition, is
For the formula (1) there results (one has N = n = 1, m = l)
AO00 = _2 + 1 A01= (2 + 4A + 1) Fl = y' f(x,y)
1 X2 h (N2 + 4A + 1)
D(7A = 5
The fact that A00 has two zeros with |i| = 1 already suggests cau-
tion, but moreover one reads off immediately that D(A) is positive
for fy < O and A = -1, and negative, in contrast, for A = -c.
Thus a zero lies to the left of -1; the method is unstable.
For the formula 5.42 in the book of Collatz mentioned (p. 81),
there is (N = n = 4, m = 1)
F4 = yIV f(x,y,y',y",y')
AOO = A22 = -(A 1)2
All = A5 = -'2 + 1
A12 = 2N
A34 = 1 (X2 + 4X + 1)
All other Ai. occur only with at least h2 in the determinant. Thus
one has, except for terms with h2
h (N2 + 41 + 1)
For A = --m,
D is positive, for ? = -1 E, D has the sign of
-A2 +- 1
+ hA + 1)
fy n< 0 and a sufficiently small
e is negative.
this method is unstable for f "n< 0.
On the other hand it is easy to indicate
stable. One need only shape the formulas (7)
every line begins with
'y ) y(i) + h a(+) l(i+) + h2
+1 k ai+-l,j j+k
methods which are always
in such a manner that
(i = 0,1, ., n-1)
Thereby A (A) = -A m+1 + Am and has therefore only the trivial zero
A = 1 on the periphery of the unit circle.
In the numerical solution of a differential equation as a difference
equation, the latter is usually of higher order and therefore has more
solutions than the original differential equation. It may well be that
some of these "extra" solutions grow faster than any solution of the
given equation; in this case the computational solution has the tend-
ency to follow one of these and has after a certain number of integra-
tion steps nothing to do with the original differential equation.
- + 1
-2 + 1
The author gives some examples and a criterion for stability of
integration methods. This criterion is then applied to some well-known
Translated by Mary L. Mahler
National Advisory Committee
1. Todd, J.: Solution of Differential Equations by Recurrence Relations.
MTAC 4, 1950, pp. 39-44.
2. Collatz, L.: Numerische Behandlung von Differentialgleichungen.
Springer, Berlin, Gottingen und Heidelberg, 1951, section 4.4.V.,
35. Milne, W. E.: A Note on the Numerical Integration of Differential
Equations. J. Res. nat. Bureau Standards 45, 1949, pp. 557-542.
NACA Langley Field. V..
- 4 3
A 0 05 0a s
0 :s., d 0 -
-d 4 0 cx a
cr 0, 00 5. k
o ,am o
CP ad 0,
0 u: C c: 00
_ CCC I 1^C s" 6B 5-
*c0 $4Cd S
vl a c 9 td P
ka V 0)^ t B w
*** fc ^ ? <
: 2 3 -
-- *4~ I- w<
.. w w a> a-
~0 r.4:,g '
Sai '- Sm S>
as 5 .-h^'-
4 0 0' 'e Q
C oS P < :d B 3
E- Z sglso
gz a g bo
S^Bl on -S ^'-
3Zf-Z rtig '
aZOKMB" a a~is
a 04 1= 0 >O
F d-0)- an 4) v
v 0 ^ :^ e u c w Z
00 =O 00 0
j e at s
hca CD ; cu
0 0 w 7c t w
^-c5 uqmc5 cula
" 3 C^ ^ f'
Smbf 3 *Z a SP aai~
a o :Z t-0.o 0
o. = cc r.. E
0if p, -rd
*Sots S.^-pn o
C I- c' 61
0.0 o~o0~weo "
ol 0 0 0 'Ig0
ar"S 0d -- 0
5 ono_ ,.
^g s a S-
4 0 a:' 4.
00'-' 0 *,0
0o N1 0
^ c o5 b^ -3
> 5 P. PI 5
L a 0 -5 o )
S M 1 >
ts o a i .'-
*= o 0 0 S 0 1
5. c -i05. I
em mu -D S a
gti blo y iao c
0 000N- 5.5
_' ffl E a d L. -
a..0 cc 0..r
Sca 0 L ( cL
io o a) J
=1 :I E!2 -SS U
'd 0 e S~
'a =1 a an 1 s
" Q 0 0 0 -a
CD 543 0 .
u v ~S 50
3.2 1 9.5, A
-0 Wa~, bo5
L, .8 o C:0 C
cd 0 = a r H
w iy sl
p -4 -i
0 a g
." g ,4J
91 *^ 3 n
o 2 S t. u
E*Sa0 N -.
0 WS < I g Sw
*e sM -
Ao 0~l 'Lo
80 s. 4
a2 g .14
4:E o =
c 0- &I
E-1 d X "
ut-S9 aa t"
0 w wo,
cc a x -
*ft ^ C ^
co B o 0- Eb
SS-u z a~o
05^, Q l
rg N gi
S 5 5a n
I- c. 4-
a a< 3 g^
4. i Nc a~~u
>4IZa <0 Co g
t C 0 'b- *UI
oiln E4 -a
= A C) 0 L
-~i-i O CB^
^'*"< a a s
~d 10 =I *
a:S cc; ts c3
S &i~ c i "
i; :95 7-i1
_E a)f. -niB B
go 0o S -
a 4.GO .c* 0-
22 "do a) O_
k d 4 2!0 >
-.5 -'Q 5 a, 21 EI
0 e 0 (U 0 C
o e 3- zo e
Ss ^ o'S SLaii "
t*ao s. V wo c 5
.2 rU o w 0
.> LUL 00 U>.
rt ^ li^-i j .Sfi'
42 4g sg>
4.. V 0 0
= = .0 d =0
.4 H 3 .4 H-
vd NM 1 Q* 14
1 f e, IV 1 a0
:s3 M ,I 14 t
2.... g~o42 44
go l u -a
4 it C s 0
044 :3 0 0
; i S S >.a CI
0 irs -c S. 10
Ced 0 9D C 0
-S 0hs s
B5 ^e S -ajgaa
E! go c roic
S'-Ss Ow c-.5 tooa
tr Id e4 ci
So z to k'
.S 4 0o
s 42 ~#4A9 g 44
to i ks1 435
k 0 C0
0- 0S i
C B ..
3 d r. 0 C3 <
o 0 -0 j
U~ .0: 00
ow, W> r- 1 :B
0.44*3'soa^' ZZ S *
w m oB -S 'tSo a
10 4 0 Go 8
55.8-6 r- vti 'p
, o 0 O__ ZE,
gS -5 |*
aP i< a'sc a Si-
- a .
E- 4 0s
a c a
5 7c S
3c am a
- g *0 V
UNIVERSITY OF FLORIDA
3 1262 08106 643 2
xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EMHV8Q4BO_9MEURU INGEST_TIME 2012-03-02T23:19:16Z PACKAGE AA00009205_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC