ZWKa-W NATIONAL ADVISORY COMMITTEE FOR AERONAUTICS TECHNICAL MEMORANDUM 1403 ON THE INSTABILITY OF METHODS FOR THE INTEGRATION OF ORDINARY DIFFERENTIAL EQUATIONS By Heinz Rutishauser Translation of "Uber die Instabilitat von Methoden zur Integration gewbhnlicher Differentialgleichungen," ZAMP, Kurze Mitteilungen, vol. Ill, 1952 UNIVERSITY OF FLORIDA DOCUMENTS DEPARTMENT I MARSTON SCIENCE LIBRARY ■ BOX 11 7011 GAINESVILLE, FL 32611-7011 USA Washington April 1956 NATIONAL ADVISORY COMMITTEE FOR AERONAUTICS TECHNICAL MEMORANDUM 1^03 ON THE INSTABILITY OF METHODS FOR THE INTEGRATION OF ORDINARY DIFFERENTIAL EQUATIONS 1 By Heinz Rutishauser^ In spite of the remarkable publication of J. Todd (ref . l) the essential points of which are related below, the author has since observed several times methods for the nummerical 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 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 Uber die Instabilitat von Methoden zur Integration gewohnlicher Dif f erentialgleichungen , " ZAMP, Kurze Mitteilungen, vol. Ill, 1952, pp. 65-7^. ^Institut fur angewandte Mathematik der ETH, Zurich. TM 1U03 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 estimates. DIFFERENTIAL EQUATIONS OF THE FIRST ORDER y' = f(x,y) Variation equation T) 1 = M TJ (a) Integration by Means of Simpson's Rule- 3 y k+ i - y k -i + 1 (y'k+i + k *\ + y'k-i) (1) This relationship, together with the differential equation, yields two equations for the unknown quantities J v ,-, 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 (l) is satisfied. The difference variation equation pertaining to (l) is evidently The phenomenon was observed on this example, and correctly inter- preted, also by Mr. G. Dahlquist, Stockholm (Lecture at the GaMM- convention 1951 &t Freiburg im Breisgau) . Compare also: Z. angew. Math. Mech., 31, 239, 1951- ^Where h signifies the length of the integration interval, and y. stands as an abbreviation for y(kh). TM li+03 1*1(1 - 1 f y , k+1 ) - f f y , k \ -(i ♦ | vjvi - ° If one assumes f to be constant, and chooses the expression t], = A for the solution of this equation, one obtains for A a quadratic equation with the solutions h 2 2 hf„ k \ = 1 + hf + — fy + . . . «, e y One recognizes easily that, of the two fundamental solutions T), , = A]_ and ^2 k = ^2 °^ 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 A 2 ~ (-1) e J i represents for f y <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 1 (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 f y < 0; there is nothing to be apprehensive of for f y = which is very important particularly in regard to the ordi- nary Simpson's integration rule (fy = 0). (b) Integration According to Runge-Kutta and Similar Methods Since these methods calculate y, ■, from y^ according to a prescribed rule and without use of the preceding values Y\_\> ^-2' • * •> TM 1U03 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. 3). (c) Integration According to Adams We consider a four-point formula y k+ i = y k + ^ fer'n+i + 19 ^'k - 5y' k _i + y\_ 2 ) + n? . . . (2) This again yields, together with the differential equation, two equa- tions for the unknown quantities y. ■, and y' k+ -i which are mostly solved by iteration. The difference variation equation pertaining to (2) becomes ^ - f f y,fcn)%n - ( x + lr f y, k ) nk + S ^k-i^-i - ^ c v_^>_ 5 = ^ 1 y,k-2 T) k-2 If one again considers f„ as constant, the expression T], = A yields an equation of the third order for A. A solution A, of this equation hfy v lies very close to e , therefore ^l k = ^1 corresponds to the ' k k solution of the differential variation equation whereas Ag and A* are extraneous solutions . However, the equation for A is reduced to A"* - A^ = when h tends toward so that, for a sufficiently small h, one will have at any rate small A2 and A*, namely ~ ±\|-hfy/2 1 +. The extraneous solu- k k tions T]2 k = ^2 and - Ix k = 'N 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 11+03 corresponding five- or six-point formulas are involved, there are no objections, but one has to be careful when y k +l is not calculated from y k and the derivatives as in (2) but perhaps from yv_]_ or yk-3 an( ^ ^e derivatives, as for instance in yk+i = y k -3 + ff ( 7y 'i-+i + 52y' k + I2y' k _! + 32y' k _ 2 + ly\. 3 \ + h T In fact the pertinent difference variation equation has a solution Ifc. A k With X=-(l-^fy+. • >j so that the method is unstable for f v < 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) y k+ i = 2 yk - y-k-i + ^(y'W + 10 A + A-i) y k+ i = y, k-i + !( y,, k + i + ^\ + y "k-i) (5) ^Compare reference 2, p. 80. TM 1U03 They yield, together with the differential equation, three equations for the unknown quantities yv + ]_> y 'k+i> an A3 = 1 \ = -fl-|f y . + . . ) where a-j_ and a,o are the solutions of the equation a 2 - af v > - f y TM 1U05 7 Evidently A]_ k and Ap k are the regular solutions of the difference variation equation; they correspond to two fundamental solutions of the differential variation equation; A* and A^, in contrast, are extraneous. As long as f 1 ^ 0, there is nothing to fear, in partic- ular, the method may be strongly recommended for a y'-free equation, ed oscillations) % k = A 4 k . (-l)V(W3)f y . but for fyi < (damped oscillations) increases, and A*, too, may still become dangerous because A* ^ = 1 also becomes finally very large, compared to a function converging to- ward zero. The author completely calculated the example y" + y 1 + 1.25y = 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 = O.l): X y y' k.8 k.9 5-0 5.1 5.2 -0.0905699 - .0847792 - .0787152 - .0722891 - .0656175 0.0551227 .0581+842 .06261+10 .0656575 .0676070 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 A^ becomes pronounced for the y' -values and also for the differences of the y-values : 8 TM 1403 X y y" 17.0 -0.00019574 0.00005017 17.1 - .00019061 .00005253 17.2 - .00018366 .00008620 17.3 - .0001752^ .00008239 17.4 - .00016548 .00011235 17.5 - .00015475 .00010258 The considerably weaker oscillation of the y-values follows also from the equations (6): for A = \ one obtains from the first of these equations — ~ — D 8 h '— = - — f,_|. here therefore p ~ — *- 4 6 y ' 600 The further course of numerical integration does not require any comment: X y y' 22.8 -0.00000815 0.00005320 22.9 - .00000864 - .00006078 23.0 - .00000862 .00005968 23.1 - .00000887 -.00006247 23-2 - .00000861 .00006601 23-3 - .00000868 - .00006486 29.5 - .00000140 - .00053037 29.6 .00000041 .00054889 29.7 - .00000144 - .00056693 29.8 .00000050 .00058682 29-9 - .00000148 .00060603 30.0 .00000060 - .00062735 The author is well aware that the assumption of a constant f y and f v i in the above considerations greatly restricts the generality. 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- gration, for instance when method (a) is being used, a special case arises since the occurring oscillations alternately increase and are damped again. TM 1U03 GENERAL CONSIDERATIONS Consideration of the examples suggests the conjecture that insta- bility may occur precisely in the case of integration methods which form y^+x by integration of y 1 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 contained in the general formula 'k+l H Oj •'k+j Z_ lj J k+j -m -m (0) (N) h > a, . y, , . -m k+l " Z_ a lj y k+j + -m -m h £. 4JV\ +J ♦ . . . + ** t « -m -m N-n+1 -m nj k+j (n-1) (N) ^j -m k+j ; (7) n-l = v^ a (n-l) (n-1) + h ^ a (n-l) (n) k+l Z_ n-l,,fk+i /=_ n.i J k+j 10 TM 1^03 There are, in addition, differential equations ■n k+i^k+i' • • •> y k +i = ° n+1 ( x k+i> • k+i' * in+1)' ? w \ x k+l' y k+l' ' •' y k+l = (N) r k+l )• (8) which form, together with the equations (7), N+1 equations for the (N) N + 1 unknowns Yv + i> ^ k+l> * " •> V ■<• Generally, N = n; however, W. E. Milne (ref . 3) 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 F , . . . to F^j by differentiating the initial equation F n . If one, furthermore, brings everything in the equations (7) to one side, the variation equations for the entire system (7) and (8) read evidently I with a. , = -1] a (i) h H-i>) = H=i j=-m HJ k+j (i = 0, 1, . . ., n - l) ^± >) . T], =0 3y(n) X (i = n, n + 1, . . ., N) (n) i If one uses for this the expression T) . = p V, one obtains in sub- stituting a system of N+1 simultaneous equations for the p which TM 1U03 11 can be satisfied only when the determinant of the system disappears |i=i\j=-m ^° / = (l = 0, 1, . . ., n - 1) X bf oy ( ^ (i = n, n + 1, . . ., N) If one defines, in addition, with the coefficients a . appearing in the formulas (7) the functions A iu (A) = 5 a (i) AJ +m wherein A. = for |i < i, the characteristic determinant reads as follows D(A) = AOO hAd All h 2 Ao2 hA 12 An-l^n-l h N AQN h N-n+l A , „ n A n-1,N *n,y * n ,y' F n+l,y F n+l,y' F n,y (n) F N,y' N,y \y W 12 TM 11+03 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 y (l)= 1 y(i+l) = y(i+2) = y(i+3) = . . . = y (N) = J J .3 J j must be compatible. Hence follows, however, by substitution into the ith of the equations (7):V~ a^ 1 ) = 0, therefore A.,.(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 = where D(A) is, except for one factor, reduced to AqcAu " ' ' ^n-1 n-1 "that A = 1 is an n-fold zero of D(A), because of A i;L (l) = 0. All other zeros of D(A) 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 = 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(A) =7"a (i KJ +m -m ^ possess, aside from the trivial simple zero A = 1 only zeros with |a| < 1. Necessary condition : None of the functions A^(A) has a zero outside the unit circle. PP^ 537-5^2. NACA • LanRley Fluid, Va. 3 c- §2 x « CO CO £ S co oj CU .C 01 X os 5 i-ga » co ^ -c ■» ■« S 3 < CO OS Z N c c . (i 3 a :s. THE LAL thod -C u a c3 -. H cu 3«ZS CO M bi> X 9 O H _ g^K g .2 ^ c dS < o u*-2 u —• 'J for ETH YDI abili r Dlff 1956. rom Z aj S OS ao Commltt ITY OF ORDINA er die In n gewohnliche user. April 03. Trans. 1 « S < z M 14C Advi INST ATIO ONS. egratlo Rutisha TM 14 NACA T National ON THE INTEGR EQUAT] zur Int Heinz (NACA 3«- tinique, (9.2 Heinz 3 as (U ■C . a™. ° S-3 ■ CO ~ Hj« > — S > rch T matic auser TM 1 hrift angewandte und Physlk, p.65-74 esea [a the utish ACA eitsc X 2 OS Z N a st a "J co a. c S jz x 55 u 2 CU 00 -^ u 6 3 01 « -£ < 10 'S — ^ X S X Z S2 CU gl> w OJ _^~ T "° co t, c "f co f Oh i - &TH m ESI rt 3 Q, -5' w °° £ cj H SB CO O g fc '- CO «> Q < O - X o H co X cu 2 X 6c so EX s <5 c X 55 Z S X o W > x'~ 3 = 5 c CU BO c I ,§« — talfl CU w O) "5b SJ — X ^ - co jr co S 3 >, c: _ .r?5 co to co x ^ j: ■ 2 5 3_- tX ■— — >* CO < z 3- 3 " £ p Oh 0)1- — c OJ X -° o : B z 2» H Z < O tt i ^ ~ _ H X. •a; co Z ZZ O S w j= x a u 0) — Q. C C <• CS = | H h/> 0) - 00 C*5 C D O O rt rr w -C ^ rt £ *-; fc* -^ ^ ou a* 3 ^ Z N X _ M 3 _ n O — 35 "! 5 3 3 00 U D 51 d O 5 &.sr "3 " * -2 1 1 IS » g 3S s s? = « ° = £■§ oZ^ rt - O -S O] : cu-2 c S __ -D S - « i J o .2 CO w J= CO J C w r 5 o - CO X 3 O J CO = 2 gS 5 , rt ^T C t- CO X S, 1 " 5 ^ 3 - Ssssp o> o £ c _ o co c = c X 2 2- J - ; CO J3 — . cu c « j* > o u C 3) ix'- m 2 Mp — rt c ^ 1o ■= ~ n c -o £ cu .2 c o £ X co c H - « x a . w 00 fE ioO g fe tH CO 3 Q u X cu X cu «s X b. SO ax 8 J ^ m « 2 < ,, c X <: co z zz o cu c j "o 3 id o j= m c g- m fct ■" cu |2 X b. : i3 Js oj tsi JHH *J 0) Q x^ •SQ H^ Li CO •-- X « SShbx B d-a-3^ CCO CJ3 1-1 CD CO ■" co x JO > rent usu tlon well if 01 CO 3 . 00 :». JS X ^ o co •e s « a Uh *is- ST 2 3 el . e ■8* S| » J. So c it S a, ■s g be u ® m o C ■" iH ■— . m S o •- H :►> a < S S U .Q C < z ■z-s < < z < Z .2 = Q Es — o u ^ w fas HI a ■c J3 M E "3 c £ a 91 S s s >. a X) c . rt a> to ao -^ — i o < < 2 s o < Z .2 => .2 S — o x: a> 5-5 «» J. So! « S & 10 i! o C -" ■"* T3 O -S >>"B. 5 « .O C . rt Q> 00 a) *j — . < o < z I C" £ C\) cu 3 a- rs m oi I cu aa cr> s 11 £ Ed c3> co" in „ -3" :3 s > H x: u « cu ■ a; X CJ a £ 1 Rutishauser NACA TM 1 Zeitschrift la" cu w s c id 3 ~ «aa w cd O u co « Q 0) eg ^ os s-a> a- a u - .2 -G co c a a > tc\3 d> m « ^ a, 3 Q — >H a « < d Z g a •3 o cu X XI -J z o r/i H Z to 5- en c 91 U ~4 ■C r cd ** j= cd cu •J] 01 5 CA itsc CU S 3 < CU « X Z N crj in 3S > Cu ^ I* < be T3 55 cd O t. CO < o u X 3 H w a 6c ax > 2 S Z S S3 = X o Cd > ^ , ra cu xi cd o co d o fie ™ 2 u S°~ cd a 5-0 C ■- C - «i « 5 <£> Q m L, 2 -cSS-S-o _ 3. 3 c <; cd : o cd ■«• 5; — .e -1 * < Z ' CU X 3 ^ Z ™ 2 ciS « _ =d £ 3 « c «; ~ T3 ^ a " * 5 CO w 5 =*> 25 c bkfl T3 en B c " ~ 2 4-1 cu -. -- -C ■J CD BO 3* 9J as p C e ri rt 3J S - > cu cu 3 CU 3 O !S 3J * - •a D a « « ci oi- ede "3 ^3 X •a cu cu — "o. 3 SP x Bu« C O o '3 «i a r- 5 o 3 ^ » ^ — m — a> cd S > u c 'Si a -8 - cd 3 § co M - " c ^J .c .3 cd •o £ J! " 5 H as e 3 c- cd in as « ?> t -2^ ' h y *! co 5 !t3 33 -^ 14/ ^3 3 S S cd H x: ■" <: u 3uS 3 < 0) « Z N ™ x: ii 3 q co cd 3 a. ri _■ aa 3 Eh iag cd O [d t. CO W 0) Q b, <06h ^ as q ^ w >* SJSf5 ■t; fo z i°l a> -h *<^, in cu cu S U) X3 ^ C35 cd o — ™ x: co w cu in *j B «> rt « ■ 2 5 N cu •- ^ «• Sa ac ■O C < cd : x; ^ u PS H ' So *> ■ . ' CO CO C 3 O O Ud" cd 2 «s - — ^ tush ii K s s <5 a g -3 5 cd a. ^ hSS o H idO §■* (1 CO c g- E : i5 It, a a ** °a « St 2 ■7 "-' x: ^5 S m U t. a "a a °S ! g 0-3 MS ^-* CQ Z C 3 O S . o rid" 9 oi S •« - * hz g -2 g : < O 603 H - M§S go -2 w> ^3 ™ 0) CO 1-0 Q co 3 « C3" 1 cu , cd .SP x 3 cd .a .a 3- »a a sr: _ cd o *5 — K cd cd a ^ 3 t< ni C3 o> -w a .* cu fcl C TJ JX : 3»C« g x: a § bo a u P £0 CO °-cd 3 C G c — .2 .2 o a o x; — t: £ ° ** CO - "" Xi fd o is CO ^ — 1 cu Cd J > ^ O bfi XI CO T3 CU x: 2 M 3 xl .S id ■o £ J! cd 3 H 2 => ~ o u *•"• S> c J; o 2 & XZ CD •2 o %5 c IS 3 «> "S £ be [0 " m o 1— t s o ■£ H >. o. " 9- < 5 * U .Q C . < z rt cd co -£ JO "> 0} w — i < < z c S sS fc. »- 0) c Jj o fc4 £ " rt 2 bo £ 0) H 5 . c to s ■o o 2 c £ S a -a c S 3 co to w — 1 < < z < z .2 a p 5; ~ o " rt 3 bo H s . c to 8 6 *5J c & 3 ^ 13 S be to * s XI c . rt co 00 *_» _* < < z .3 3 s z a> B O u G) ■ & -C QJ H C a to s ■a JO j* 0) e "a! c is QJ id s & CO u O c •O <*.< >i 1 —j a. s rt q c . 3 CU CO CO -e « o bo < < z 3 1262 08106 643 2 UNIVERSITY OF FLORIDA DOCUMENTS DEPARTMENT 120 MARSTON SCIENCE LIBRARY P.O. BOX 117011 GAINESVILLE, FL 32611-7011 USA