UFDC Home  Search all Groups  World Studies  Federal Depository Libraries of Florida & the Caribbean   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
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 welltested methods of RungeKutta 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 enumerated. 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, pp. 6574. 2Institut fur angewandte Mathematik der ETH, Zurich. TM 1403 deviations are unavoidable, because of the roundingoff 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 estimates. DIFFERENTIAL EQUATIONS OF THE FIRST ORDER y' = f(x,y) Variation equation Cy (a) Integration by Means of Simpson's Rule5 yk+i = Yki 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). TM 1L03 flk+l( yk+) y k 3 1 + f ,k k1 = 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 h2 2hfy 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. In particular 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 roundingoff 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 RungeKutta and Similar Methods Since these methods calculate yk+l from yk according to a prescribed rule and without use of the preceding values ykl' Yk2' *f TM 1405 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 (ref. 5). (c) Integration According to Adams We consider a fourpoint formula k+l = k + rk+l + 1y'k 5'Yk + y'k2) + (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,klkl y,k2nk2 = 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 hfy k 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 k TM 1405 5 corresponding five or sixpoint 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 Yk3 and the derivatives, as for instance in Yk+1 = Yk5 + '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 second order) k+1 = 2yk Yk1 + j"k+l + 10"k + k 4) y' y'k + k.l + + Y" k) (5) 5Compare reference 2, p. 80. TM 1403 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)  for h2fy/l2 and b for q S1 b(A2 + 10A + 1) = 0 4 q (2 1) b(V2 + 4A + 11 p p a(A2 + 4A + 1) = 0  (6) [from (5] 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 afyr fy = 0 [from (4)] TM 1403 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 ward zero. The author completely calculated the example y" + y' + 1.25y = 0 with the initial conditions y. = 0, y' = 1 (exact solution: x e sin x on the sequencecontrolled 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 yvalues 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 the yvalues: 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 TM 1403 The considerably weaker oscillation of the yvalues follows also from the equations (6): for A = X4 one obtains from the first of t;iese equations 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. and How 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 damped again. 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 GENERAL CONSIDERATIONS 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 contained in 0 (0) yk+ = i_ 4j Yk+j m + h m a(0) yk+j + ajj Yvk+j + + N (0) (N) h alj Yk+j m 0 y k+1 2 m (1) alj Y'k+j + h m 2j k+j N1 1 (1) (N) Y: Ni k+j n1 + a(nl) (n1) k+l m nl,j k+j + h m a(n1) (n) nj k+j Nn+11 (n1) (N) h k+j + . + .. + TM 1O05 There are, in addition, differential equations Fn (k+lyk+' k+ ) = F (+1(xk k (n+l) =0 FN(xk+l'Yk+l * * y( = 0 k+1)  which form, together with the equations (7), N + 1 equations for the (N) 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 evidently with ih a  P j>t h ak = O p=3 j=m 0 Fi (p o ay(4) If one uses for this the expression q. j (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 TM 1403 can be satisfied only when the determinant of the system disappears N 1 r) 7hip = 0o 4=1i =L1 1 OF  P, = 0 ay () (1 = 1, 1, ., n 1) (i = n, n + 1, ., N) (i) If one defines, in addition, with the coefficients a the formulas (7) the functions Pj appearing in A () = a (i)j+m j=m wherein A. . follows 0 for p < i, the characteristic determinant reads as hAO1 All "* 0 0 'An1,n Fn,yI Fn,y(n) 0 hNAON hN1AUN hNn+1An1,N 0 Fn+1,y Fn+l,y' . . FN,y(N) FNy FN,y AOO 0 I F 0 Fn,y D(x) = . . . TM 1405 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 An1, nl, that 'A = 1 is an nfold 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: All functions Aii(D) = (i)ij+m m ij 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. TM 1405 The remaining functions Al(AX) can influence the stability only when the necessary condition, but not the sufficient condition, is satisfied. APPLICATIONS 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) Therefore 1 X2 h (N2 + 4A + 1) D(7A = 5 fy 1 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 TM 1405 0 0 0 h (N2 + 41 + 1) 1 For A = m, D is positive, for ? = 1 E, D has the sign of A2 + 1 h (2 + hA + 1) 1 fy n< 0 and a sufficiently small e is negative. Therefore 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 Mm methods which are always in such a manner that (i = 0,1, ., n1) Thereby A (A) = A m+1 + Am and has therefore only the trivial zero A = 1 on the periphery of the unit circle. SUMMARY 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. D(A) = (A 1)2 0 0 0 fy 0  + 1 0 0 f (A 1)2 0 y 0 0 0 2 + 1 f y" Which for TM 1403 The author gives some examples and a criterion for stability of integration methods. This criterion is then applied to some wellknown integration formulas. Translated by Mary L. Mahler National Advisory Committee for Aeronautics REFERENCES 1. Todd, J.: Solution of Differential Equations by Recurrence Relations. MTAC 4, 1950, pp. 3944. 2. Collatz, L.: Numerische Behandlung von Differentialgleichungen. Springer, Berlin, Gottingen und Heidelberg, 1951, section 4.4.V., P. 51. 35. Milne, W. E.: A Note on the Numerical Integration of Differential Equations. J. Res. nat. Bureau Standards 45, 1949, pp. 557542. NACA Langley Field. V.. * ci Z.0 co c  4 3 ,,Bl~i r al3;l i? l" sIS 8d9< gs1 t;S KZNrt )~JN ~42 17 A 0 05 0a s 0 I0 ~~.0 0 :s., d 0  419 d 4 0 cx a  Cc cr 0, 00 5. k o ,am o CP ad 0, 0 u: C c: 00 _ CCC I 1^C s" 6B 5 0 8 *c0 $4Cd S vl a c 9 td P ka V 0)^ t B w *** fc ^ ? < ig~laj^" 5grL'ijBCl el . 0d : 2 3  c a bo  *4~ I w< .. w w a> a ~0 r.4:,g ' M Ej 000 Sai ' Sm S> as 5 .h^' 4 0 0' 'e Q 00 00 = 5:05 *.0 C oS P < :d B 3 E Z sglso gz a g bo 0 0i. SEtZ !G30. oa~i?;.ortas S^Bl on S ^' 3ZfZ rtig ' aZOKMB" a a~is a 04 1= 0 >O 00 F d0) an 4) v v 0 ^ :^ e u c w Z 00 =O 00 0 j e at s r bO ~,MIU, hca CD ; cu 0 0 w 7c t w ^c5 uqmc5 cula " 3 C^ ^ f' Smbf 3 *Z a SP aai~ a o :Z t0.o 0 0 o. = cc r.. E 0if p, rd *Sots S.^pn o C I c' 61 ^ii2imi &+3 0.0 o~o0~weo " ol 0 0 0 'Ig0 id 3OEOz ar"S 0d  0 ; id 'i m 0 .*  5 ono_ ,. ^g s a S QSQ5 W:2 0 Wlg IDpl S""'LS 4 0 a:' 4. 00'' 0 *,0 CS><5Q "'c uI a~ 4' 4 0o N1 0 ^ c o5 b^ 3 > 5 P. PI 5 or mma L a 0 5 o ) wzc ~OWNEcm "3 * S M 1 > ts o a i .' 0D *= o 0 0 S 0 1 5. c i05. I em mu D S a gti blo y iao c 0d 0 000N 5.5 _' ffl E a d L.  a..0 cc 0..r (u r Sca 0 L ( cL io o a) J 0 =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 ^^s~sgls w iy sl o~slS~PiS 1 I Q"i ii I=r " on.  0  0 5& .a .a g i 5 a> S 40 44) 0l 1.. ,0c: Bo "o . 4 o  .0 S e tia Sl oS mS P  p 4 i .44 S S 0 a g ." g ,4J 91 *^ 3 n o 2 S t. u 4242 t E*Sa0 N . 42a aa 0 WS < I g Sw &Mw N .^gl=i Ns *e sM  Ao 0~l 'Lo op4 80 s. 4 a2 g .14 4:E o = c 0 &I E1 d X " UR utS9 aa t" 0 w wo, evi 0) E4r IV: cc a x  ss as *ft ^ C ^ co B o 0 Eb SSu z a~o 05^, Q l 2.S N3*o'< N42 V 142 V s Is., . N rg N gi S 5 5a n o~~d Gi I c. 4 a a< 3 g^ 5wd  r u 4. i Nc a~~u 0 I ,OJ > > 1C o 0 >4IZa <0 Co g OE.0 L)  "'g~aM S t C 0 'b *UI itZ 7 2S^.g> ~4:Z C=rO: 4:zco ~~O4et Q20 oiln E4 a = A C) 0 L g>.got 5^1 ^^ s~ & gSm~n 532.3 ~ii O CB^ ^'*"< a a s Sa gS5!"! ~d 10 =I * SdS6i g'l a:S cc; ts c3 S &i~ c i " i; :95 7i1 _E a)f. niB B gd Cd go 0o S  a 4.GO .c* 0 22 "do a) O_ k d 4 2!0 > 2to .5 'Q 5 a, 21 EI 0.E 0 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> U00 g~js4 4.. V 0 0 = = .0 d =0 (uC ~p~"U .4 H 3 .4 H u 0 1s 0 u 0r a S h ao Un gg 4: (U z~o ' A 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 O w S 0hs s B5 ^e S ajgaa E! go c roic 42 S'Ss Ow c.5 tooa fi"ilS:g? tr Id e4 ci So z to k' 4044 .S 4 0o s 42 ~#4A9 g 44 to i ks1 435 'D 00 k 0 C0 (U w 0 0S i C B .. 0 ba~ o0 t 3 d r. 0 C3 < o 0 0 j U~ .0: 00 kz ow, W> r 1 :B 0.44*3'soa^' ZZ S * w m oB S 'tSo a 10 4 0 Go 8 55.86 r vti 'p , o 0 O__ ZE, gS 5 * ^li~s'sl^ aP i< a'sc a Si 1 t o  a . 0 0 E 4 0s C a .0 cW aa a a c a 5 7c S E & v aa a8 a Sca 3c am a a a a a Sa as SI  g *0 V UNIVERSITY OF FLORIDA 3 1262 08106 643 2 
Full Text 
xml version 1.0 encoding UTF8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd INGEST IEID EMHV8Q4BO_9MEURU INGEST_TIME 20120302T23:19:16Z PACKAGE AA00009205_00001 AGREEMENT_INFO ACCOUNT UF PROJECT UFDC FILES 