L I B R.AFLY OF THE "~'- : UNIVERSITY , OT ILLINOIS . 510. 84 \£6r , no. 524- 33d cop. ,2 !N. The person charging this material is re- sponsible for its return on or before the Latest Date stamped below. Theft, mutilation, and underlining of books are reasons for disciplinary action and may result in dismissal from the University. University of Illinois Library m 1 3 m ■ K .. Jwmwflff 2 a 5 ll&fEBl MAY 2 9 RpD OCT 2 \ WO NOV 1 4 1970 NOV m ■ 7 i OEC 19MEB SEP 3 1977 - NOV 1 2 £74 JJ^V2 JUL 1 19 tS FEB - 1 1971 JAN 2 <> S£Pl "W 2 5 flat) MAR 8 1978 MAR 8 RECD m a « MWR ?f UK! APR 2 » MAR i 4 1978 1REFI 2004 A3T7 L161— O-1096 Digitized by the Internet Archive in 2013 http://archive.org/details/studyofstability327vohr Report No. 327 ri^L^Oii coo-ii+69-0118 STUDY OF STABILITY FOR NON-LINEAR ORDINARY DIFFERENTIAL EQUATIONS by Raj Vohra June 1969 JUM i ' ttil DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN URBANA, ILLINOIS Report No. 327 STUDY OF STABILITY FOR NON-LINEAR ORDINARY DIFFERENTIAL EQUATIONS* by Raj Vohra June 1969 Department of Computer Science University of Illinois Urbana, Illinois 6l801 This paper was supported in part by Grant U. S. AEC AT(ll-l)lU69 and submitted in partial fulfillment of the requirements for the degree of Master of Science in Computer Science, June 1969 . iii ACKNOWLEDGEMENT The investigation described herein was undertaken as a Master's thesis project under the direction of Professor C. W. Gear of the Department of Computer Science of the University of Illinois. The author wishes to express his appreciation to his advisor, Professor Gear, for his guidance and comments during the investigation. IV TABLE OF CONTENTS Page 1. INTRODUCTION 1 2. DESCRIPTION OF METHODS USED AND RESULTS U 3. SUMMARY 3*+ LIST OF REFERENCES 35 1 . INTRODUCTION Predictor-corrector methods furnish attractive algorithms for the solution of ordinary differential equations because of the relatively small number of derivative evaluations required per step. For example, fourth degree predictor-corrector methods require two derivative evaluations per integration step while the corresponding Runge-Kutta fourth degree algorithm requires four derivative evalu- ations per integration step. The difficulty of missing starting values is typical for all multistep methods. One step methods such as Runge-Kutta methods or Taylor's expansion are frequently used to find missing starting values. Other starting methods which are more within the spirit of difference methods are also known (Collatz, "The Numerical Treatment of Differential Equations," 3rd Ed. Springer, Berlin (i960), p. 8l). One of the key factors to be considered in the selection of a particular predictor-corrector pair is the stability of the numerical algorithm. This is particularly crucial when the differential equations being solved correspond to a system with a forcing function whose time duration or period is relatively long compared to the transient time constants of the system. Instability is also introduced by the fact that a first order equation is solved with the help of higher order difference equations. The principal root of the difference equation approximates the true solution of the equation. Any other root is extraneous. Aside from these parasitic roots, presence of nonlinearity in a differential equation can cause instability, while establishing the region of stability for a particular method this is generally ignored. In this work we will try to study the effect of nonlinearity in an ordinary differential equation. Adams -Bashforth predictor and Adams-Moulton corrector of orders three through six will he studied and regions of their absolute and relative stability established for different values of parameters. This study is not exhaustive because only a few values of the parameters will be considered. We will introduce the following terminology to describe the purpose of this study in more detail. y' = f(x, y) a < x < b y(a) = y g(x) = f y (x, y(x)) g n = g(x n ) hg n = h and g, (x n ) = <*g^ We will ignore higher derivatives of g(x). Define a function such that 0(o, h) = [0 if stable j • 1 if unstable We will establish i. The region of numerical stability for Adams -Bashforth predictor and Adams -Moult on corrector of orders three through six for a ^ and compare it with the corresponding region of stability for a = 0. ii. It was conjectured that region of stability for a ^ for a particular predictor-corrector pair lies inside the region of stability for a = 0. This was not so in some examples . 3 iii. Another conjecture we wished to test was whether 0(0, h(l-J2h)) = for j = 0, 1,..., k for a k step method. =s> 0(a, h) = This was also false. In the following pages these three properties will be studied in some detail. 2. DESCRIPTION OP METHODS USED AND RESULTS All the methods we use for solving differential equations have various possible sources of error. Circumstances make it likely that the true solution differs to some extent from our computed solution, and we try to make the discrepancy as small as required for any particular application. When we try to estimate the errors, a number of factors must be taken into account. They are i. Inherent instability in the problem. ii. Instability introduced by our particular method of solution, iii. Complete or partial neglect of the truncation error, iv. The presence of some form of singularity, v. Nonlinearity in a differential equation. The problem of inherent instability is independent of the method of solution. We will not consider this problem in this work. We will mainly be concerned with the effect of nonlinearity on the stability of predictor-corrector methods. To get trustworthy results stability with a large margin of safety under all circumstances is essential. In lengthy automatic computations the number of elementary steps may be as large as 10 or more, and disturbances in unstable methods typically grow exponentially with the number of steps. Predictor-corrector methods may cause more catastrophic strong unstability through the introduction of spurious solutions . A result of this kind is hardly surprising since the higher order difference equation must introduce spurious solutions having no relation to that of the first order differential equation. The important feature is 5 that the solution is not improved by a reduction in the size of the interval. The parasitic solutions exist. Another important factor which can cause numerical instability is nonlinearity in an ordinary differential equation. We will attempt to show this in this work. Suppose that the initial value problem (2.1) y ' = f( x , y) a < x < b yCa) = y Q where f is continuous and satisfies Lipschitz conditions with respect to y, is to be solved by a predictor-corrector method of the type PECE, i.e. calculate a predicted value of y, evaluate f(x, y) at the predicted value, correct y, and evaluate f(x, y) at the corrected value. The additional starting values of y which are needed to start the solutions are assumed to have been calculated by some other method. The predicted value y P of y at x = x , is calculated r "'n+k J n+k from the formula (2 - 2) Ck + «k-l *n+k-l + "-- + Vn " h[b k-l Wl + "- + Vn ] where f stands for f(x ,,, y .).y=0, 1,... k-1. The function n+y n+j ' J n+j J f is then evaluated at x ,, , y ,, to give n+k' J n+k & f P = f(x v P ) n+k IV n+k' y n+k ; The corrected value y , of y is then calculated from the formula. J n+k J (2 ' 3) ^n + k + C k-1 y n+ k-l + --- +C & = h[d k f n + k + \-l W-l ♦— ♦Vnl- Finally, the function f is again evaluated at x , y to give f n+k = f(x n+k' y n +k } where predictor-corrector pair is of suitable order. Henrici defines the error in each predicted value by e ~ = y 5 - y(x ) and error in each corrected value by e = y - y(x ) m m m m m m where m = 0, 1, 2,... and he shows that e~ , satisfies UM ** '- -jio [U j " hb J w ^j 1 " Vi hP+1 *** ll *J * 0(h* +2 ) and that e n+k satisfies k-1 (2 ' 5) e n + k = ~l [C j " hd j Wj ] e n + j ] + hd k ^n+k e n + k " C q+ 1 h4+1 y " ^ Vj « k Vj + J W - h6 Vj Vj W 1 £ n +J - c ; + i \ \ + ^ + » 2 w hP+2 * (p+1)(5) - C p+1 h* +1 y ( 1 +1) (U + 0(h^ 2 ) In conventional methods to find the region of stability g! is taken to he zero. While solving a differential equation absolute values of g' and g are not small and thus may have a marked affect on the stability of solution. In the literature one can find root locus plots of the characteristic equations of the difference equation (2.8) giving the magnitudes of roots plotted against E * hg is constant. Krough gives these plots for Adams- Bashforth and Adams -Moulin type methods. In this work we will take g to be constant. In the analysis that is given here, third and fourth order terms in the righthand side of (2.15) are disregarded. The characteristic equation of the simplified difference equation is k— 1 (2.16) Y k + Z [[C + hg n (ad - d ) - hV (db ) 2 «I t AA _ V * M ~dl . - h ^ ljd j - k ajV ] Y 1 = _ 2 To study the region of stability we define hg = h and g'/g^ = ot in (2.16) and assume that g/ i- 0. Equation (2.l6) thus becomes k k_1 2 (2.17) Y + £ [[C, + h (a.d, - d.) - h" (d,bj i-Q J J -K- J K J - h 2 a(jd - k i)] Y" 3 ] = Both a and h are taken to be complex valued. Of interest are those values of h and a in (2.17) which will make the roots y = Y-> ip = 1> 2,..., k of this polynomial satisfy the condition of stability. If |y. | £l, i = 1, 2,..., k and roots of unit magnitude are simple, the method is 10 said to be absolutely stable. The region of absolute stability is defined to be the set of points in h plane (for particular value of a) for which the method is absolutely stable. One of the roots ' y ' of the characteristic equation, the P principal root, is approximately e . Any other root we label Y e "to indicate that it is an extraneous root. The extraneous roots result from approximating a first order differential equation with a difference equation of higher order. If each of y satisfies |y_| < e , and if the y with the magnitude of e are simple, then the method is said to be relatively stable. The region of relative stability is defined to be the set of points in the h plane for which the method is relatively stable. Inside this region introduced [7] errors grow (or decrease) at the same rate as the solution. Ralston has defined a method to be relatively stable provided all extraneous roots satisfy \y <_ \ y and those y with magnitude of y are simple. In this region useful solutions can be obtained and are computationally more convenient. To keep computations at a reasonable level only eight values of a are considered for each pair of predictor-corrector methods. (2.18) a = 0-5e ij7r/U j = 1, 2,..., 8 Lehmer technique is used to test whether all roots of the characteristic equation lie inside the region |z| <, 1. The method is as follows . 11 In characteristic equations substitute y = l/z. Let the resulting polynomial be denoted by (2.19) g(z) = a Q + a 1 z + ...+8^ where the coefficients a. are complex. Let (2.20) g~(z) = a + a + z +...+a„ z B n n The linear combination T(g) = T(g(z)) = a Q g(z) - a n g (z) has no terms in z so that T(g) is of lower degree than g. The constant term of T(g) is of particular interest because it is real. In fact T(g(0)) = a Q a - a n a n = I & r I |2 |2 1 If this constant term is different from zero T(g) is a fit polynomial to apply this transformation. Continuing as far as possible we obtain the finite sequence T(g), T 2 (g),... -Ag) where ^g(O)) * If for same i > 0, T (g(0)) < 0, then g has at least one root inside the unit circle. This implies that equation (2.17) has at least one root outside the unit circle. In this way we are able to determine those points for which there is no root of g(z) inside the unit circle. This provides us the region of absolute stability. This algorithm also provides a convenient method to determine the region of relative stability. If the principal root = y then we make the transformation 12 z = z/y in equation (2.19) and follow the above procedure, Methods Investigated Adams -Bashforth predictor and Adams-Moulton corrector of orders three through six are considered. TABLE 1. Coefficients of Adams -Bashforth Predictor Ck + a k-l y iwk-l + --- + Vn = h [\-l f n + k-l + --- + Vn' + Vl hP+lyP+1 a, = -1, all other a's are zero. Fig. 5-12 Fig. 13-20 Fig. 21-28 Fig. 29-36 p 3 h 5 6 b o 5/12 -9/2U 251/720 -U25/1UU0 b l -16/12 3T/2U -127^/720 2627/lUUo b 2 23/12 -59/24 2616/720 -6798/lUUo b 3 55/24 -277^/720 9U82/1UU0 \ 1901/720 -7673/lUUo b 5 U227/1UU0 13 TABLE 2. Coefficients of Adams -Moult on Corrector y + C y +...+C y = h[d, f 5 , + d, , f , , +...+d n f ] °n+k k-1 y n+k-l 0°n k n+k Tc-1 n+k-J. On ♦ C + _ h* +1 y4 +1 (x n ) + 0(h* +2 ) . , = -1, all other C's are zero Fig. 5-12 Fig. 13-20 Fig. 21-28 Fig. 29-36 q. 3 1+ 5 6 d o d l -1/12 1/24 -19/720 27/1440 d 2 2/3 -5/24 106/720 -173/lUUo d 3 5/12 19/24 -264/720 482/1440 \ 3/8 646/720 -798/1440 d 5 251/720 1427/1440 d 6 475/1440 XX NOTE: x x gives boundary for the region of absolute stability • ■ gives boundary for the region of absolute and relative stability oo oo gives boundary for the region of relative stability lU h-PLANE \» / / FIG.I /' a = / ^"""" '/ k = 3 / y^ * • i > / x / I • s / I / / / s / t ' s / t i s / i [.*' • i / i / 1 ..J I , . -3-2-10 I 2 THIRD ORDER ADAMS- BASHFORTH PREDICTOR WITH ADAMS- MOULTOR CORRECTOR + 2i J L FIG. 2 a=0 k = 3 -3-2-10 I 2 FOURTH ORDER ADAMS -BASHFORTH PREDICTOR WITH AOAMS- MOULTOR CORRECTOR h-PLANE -•2i -L. FIG. 3 a = k = 5 15 -3-2-10 I 2 FIFTH ORDER ADAMS- BASHFORTH PREDICTOR WITH ADAMS - MOULTON CORRECTOR -2i FIG. 4 a=0 k=6 •■i -3-2-10 I 2 SIXTH ORDER ADAMS - BASHFORTH PREDICTOR WITH ADAMS- MOULTON CORRECTOR 16 -2 FIG. 6 k=3 THIRD OROER ADAMS -BASHFORTH PREDICTOR WITH ADAMS- MOULTON CORRECTOR 2i t THIRO ORDER ADAMS -BASHFORTH PREDICTOR WITH ADAMS- MOULTON CORRECTOR FIG. 9 a = 0.5t i5 7% k = 3 FIG. 10 a=0.5e'3^ k= 3 THIRD ORDER ADAMS -BASHFORTH PREDICTOR WITH ADAMS -MOULTON CORRECTOR -2 19 FIG. 1 1 k = 3 2i t THIRD OROER ADAMS -BASHFORTH PREDICTOR WITH ADAMS -MOULTON CORRECTOR h-PLANE 2i 20 FIG. 13 k=4 2i t -2 -I FIG. 14 a=0.5«'3% k = 4 FOURTH ORDER ADAMS -BASHFORTH PREDICTOR WITH ADAMS -MOULTON CORRECTOR h-PLANE 2'it 21 FIG. 15 a = 0.5« i3 ^ k = 4 2i t -2 -I i ■• X* : X*. FIG. 16 a = 0.5« i7r k = 4 FOURTH ORDER ADAMS -BASHFORTH PREDICTOR WITH ADAMS -MOULTON CORRECTOR FIG. 17 a = 0.5t i5 % k = 4 -2 -I 2i t i •• xocx:x FIG. 18 a = 0.5e'37^ k = 4 FOURTH ORDER ADAMS -BASHFORTH PREDICTOR WITH ADAMS - MOULTON CORRECTOR 23 FIG. 19 a = 0.5t 17 ^ k = 4 -2 2i t -2 FIG. 20 a=0.5«'27T k=4 FOURTH ORDER ADAMS -BASHFORTH PREDICTOR WITH ADAMS - MOULTON CORRECTOR h-PLANE 2J T 2k FIG. 21 a = 0.5«i"% k= 5 -2 -I 2ir -2 FIG. 22 a=0.5«i "^ FIFTH ORDER ADAMS -BASHFORTH PREDICTOR WITH ADAMS - MOULTON CORRECTOR h-PLANE 2i T 25 FIG. 23 k=5 2i t FIG. 24 a=0.5t l7r k = 5 FIFTH ORDER ADAMS -BASHFORTH PREDICTOR WITH ADAMS -MOULTON CORRECTOR h-PLANE 2i T 26 I ■■ FIG. 25 k = 5 X : &#£4 < « 1 *70vsfl -2 -1 — *E$S K : . : **¥ *W* #! *!*:* ft" 2i T -2 i -• FIG. 26 a = 0.5t'37J> k=5 FIFTH ORDER ADAMS -BASHFORTH PREDICTOR WITH ADAMS-MOULTON CORRECTOR -2 -I 27 FI6.27 k= 5 2i t FIG. 28 a = 0.5t' 27T k= 5 -2 -I FIFTH OROER AOAMS-BASHFORTH PREDICTOR WITH ADAMS - MOULTON CORRECTOR -2 -I h-PLANE 2i T :*:*:■* \ >p •*:*:; FIG. 29 k = 6 28 2i T -2 -I FIG. 30 k = 6 SIXTH ORDER ADAMS -BASHFORTH PREDICTOR WITH ADAMS -MOULTON CORRECTOR 29 FIG. 31 a = 0.5e*3^ k= 6 -2 -I 2i T FIG. 32 a=0.5e' ,7r k = 6 -2 -I SIXTH ORDER ADAMS -BASHFORTH PREDICTOR WITH ADAMS- MOULTON CORRECTOR -2 -I 30 FIG. 33 a*0.5t 15 ^ k»6 -2 -I 2i T i ■■ x* ■*.* ! > x*x FIG. 34 a=0.5« l3 ^£ k = 6 SIXTH ORDER ADAMS -BASHFORTH PREDICTOR WITH ADAMS - MOULTON CORRECTOR h-PLANE 31 2'it i " FIG. 35 a = 0.5fl77^ k = 6 -2 -I 2i T ***** FIG. 36 a = 0.5«i27r k = 6 -2 -I Ktt:* SIXTH ORDER ADAMS -BASHFORTH PREDICTOR WITH ADAMS - MOULTON CORRECTOR 32 The curves of Figure 1-k result from solving the characteristic equation for a = 0. With the symmetry in mind, only the parts of the curves in the upper half plane are given. Portions to the right of solid curves are the regions of absolute stability and portions to the right of dashed curves are the regions of relative stability. For a $ stability was calculated at an interval of «1 along x-axis and y-axis. The figures 5-36 result from solving the characteristic equation for a = o»5 ^ ' j = 1, 2,... 8 for different predictor-corrector pairs. Since in this case there is no symmetry, complete figures are given. When we compare the region of stability for a = o and for a f it is quite obvious that even a small value of a has a marked effect on the region of stability • For any predictor-corrector pair it was conjectured that regions of absolute (relative) stability for a # will always lie inside the region for a = 0. But this is not true. Compare Figure 22 with Figure 3. In Figure 3 the point + . 8i is not stable where as in Figure 22 point + . 8i is stable. But in general, the region of stability shrinks when a 4- 0. Another conjecture which is disproved is the following. We have defined a function 6(a, h) s if stable 1 if unstable 6(0, 5(1 - jah)) =0 j = 0, 1,... k => 9(a, h) = 33 Absolute Stability . In Figure lU o a ..'5i. The point E = -.1 + .3i is unstable or 8(a, E) = 1 but 6(0, 5(1 - Jah)) =0 j = 0, 1,... h from Figure 2. k = h for both of these figures. Relative Stability . Consider again Figure lU, a = .5i. The point E = +.1 - . i is unstable but 6(0, E(l - jotE)) = J = 0, 1,... k from Figure 2. Another interesting fact to note from the study of these figures is as the order of method is increased the region of stability narrows down considerably. 3k 3. SUMMARY For any specified predictor-corrector pair the a's, b's, C's, and d's in (2.17) are known. We have seen that stability depends upon the roots of characteristic equations. We have replaced « . = g(x +j) = £L+ j^gl. We neglect the higher order terms. It is not clear how the neglected terms will effect the situation. When we compare the region of stability for a ^ with a = it is evident that the region of stability shrinks in most cases. Suppose that in the course of the solution, estimates of g and g* for each n are known and it is found that neither is small in absolute value. Then it is very likely that the variation in g and g* can have a marked effect on the numerical stability of the solution. In higher order methods the region of stability is quite smaller than the region of stability for lower order methods. Through the use of Figures U-36, it is possible to predict stability characteristics of predictor- corrector methods as a function of h and a. 35 LIST OF REFERENCES [l] Chase, P. E. "Stability Properties of Predictor-Corrector Methods, JACM, 9, 1962. [2] Gear, C. W. "Hybrid Methods for Initial Value Problems in Ordinary Differential Equations" SIAM Journal on Numerical Analysis, 13, #2, 1965. [3] Hamming, R. "Stable Predictor-Corrector Methods for Ordinary Differential Equations ," JACM, 6, 1959- [k] Henri ci, P. "Discrete Variable Methods in Ordinary Differential Equations," John Wiley, New York, 1962, p. 260. [5] Krogh, F. T. "Predictor-Corrector Methods of Higher Order with Improved Stability Characteristics, JACM, 1965. [6] Lambert, R. J. "An Analysis of the Numerical Stability of Predictor-Corrector Solutions of Nonlinear Ordinary Differential Equations ," SIAM Journal on Numerical Analysis, 4, #6, 1967. [7] Ralston, A. "Relative Stability in the Numerical Solution of Ordinary Differential Equations ," SIAM Review, 7, 1965, p. 114-125. Form AEC-427 (6/68) AECM 3201 U.S. ATOMIC ENERGY COMMISSION U r^'o^o S,TY " TYPE CONTRACTOR'S RECOMMENDATION FOR DISPOSITION OF SCIENTIFC AND TECHNICAL DOCUMENT ( See Instructions on Reverse Side ) 1. AEC REPORT NO. 000-1469-0118 2. TITLE Study of Stability for Non-Linear Ordinary Differential Equations 3. TYPE OF DOCUMENT (Check one): a. Scientific and technical report □ b. Conference paper not to be published in a journal: Title of conference Date of conference _^ Exact location of conference Sponsoring organization □ c - Other (Specify) RECOMMENDED ANNOUNCEMENT AND DISTRIBUTION (Check one): j3 a. AEC-s normal announcement and distribution procedures may be followed n c Mat n Vailable ° n,V Withi " AEC 3nd t0 AEC COntraCt ° rS and ° th6r US - G ° V — l -"*- "* ** contractor, LJ c. Make no announcement or distrubution. 5. REASON FOR RECOMMENDED RESTRICTIONS: 6. SUBMITTED BY: Organization Signature NAME AND POSITION (Please print or type) C W. Gear, Professor of Computer Science and Applied Mathematics Department of Computer Science University of Illinois Urbana, Illinois 6l801 ^hJ2^k> C i<^'' Date June 2, 1969 7. AEC CONTRACT ADMINISTRATOR'S COMMENTS IF FOR AEC USE ONLY RECOMMENDATION: ANY. ON ABOVE ANNOUNCEMENT AND DISTRIBUTION PATENT CLEARANCE: Q a. AEC patent clearance has been granted by responsible AEC patent group. U b. Report has been sent to responsible AEC patent group for clearance. U c. Patent clearance not required.