The person charging this material is re- sponsible for its return to the library from which it was withdrawn 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 AT URBANA-CHAMPAIGN I DEC 16 iEW L161 — O-1096 Digitized by the Internet Archive in 2013 http://archive.org/details/optimalstifflyst402jain REPORT NO. 1+02 ^f / n c^aa\ OPTIMAL STIFFLY STABLE METHODS FOR ORDINARY DIFFERENTIAL EQUATIONS by M. K. Jain V. K. Sri vast ava COO-1U69-016U June 1970 IHE LIBRARY OF THE JUL 3 1 1 UNIVERSITY OF ILLINOIS REPORT NO. H02 OPTIMAL STIFFLY STABLE METHODS FOR ORDINARY DIFFERENTIAL EQUATIONS* *y M. K. Jain V. K. Srivastava June 1970 Department of Computer Science University of Illinois Urbana, Illinois 6l801 * Supported in part by grant U. S. AEC AT(ll-l)lU69 and by the National Science Foundation US NSF GJ 812 . ERRATA TO REPORT NO. 1+02 COO-lU69-Ol6U OPTIMAL STIFFLY STABLE METHODS FOR ORDINARY DIFFERENTIAL EQUATIONS by M. K. Jain V. K. Srivastava June 1970 p. 8, 1. -2 171 13 9 L , 300 , + 726 Y n-k + 726 y n-5 " 726 y n -6 h 726 y n-l ACKNOWLEDGMENT We wish to thank Professor C. W. Gear for his helpful suggestions. Thanks are also extended to Miss Barbara Hurdle for typing of the final manuscript. iii TABLE OF CONTENTS Page^ ABSTRACT iv LIST OF FIGURES. . . v 1. INTRODUCTION 1 2. OPTIMAL STIFFLY STABLE METHODS 2 3. CONCLUSIONS 10 LIST OF REFERENCES 11 IV ABSTRACT The purpose of this report is to discuss one parameter methods for the numerical integration of stiff differential equations, i.e., systems with widely separated time constants. The methods which are optimal in terms of least amount of solution history data to "be saved at each step of the corrector iteration have been known for order as high as six. We find that in one parameter methods the stability and accuracy are controlled by the same variable y. Suitable values of y lead to methods which have either optimal accuracy or optimal stability. The methods of order six have been presented. The section of the locus p(£)/o(0 for methods of order six is shown in Figure 7. LIST OF FIGURES FIGURE Pape 1. STABILITY AND ACCURACY REGIONS 2 2. THIRD ORDER FORMULA 12 3. LOCUS OF p(£)/a(5) FOR THIRD ORDER FORMULA 13 k. FOURTH ORDER FORMULA lU 5. FIFTH ORDER FORMULA 15 6. SIXTH ORDER FORMULA 16 7. SECTION OF LOCUS OF p(€)/a(0 FOR p < 6 IT 1. INTRODUCTION Many numerical problems in applied mathematics entail the solution of systems of ordinary differential equations with a property given by the following definition. Definition. A system of ordinary differential equations y' = f(t y y) 3 y(a) = n is said to be stiff if the eigen values of the matrix (-^-) has at every -point t negative real parts, and differ greatly in magnitude. Equations of this form frequently arise in physical equilibrium problems such as chemical or nuclear reactions where there are many components with a large range of time constants. The standard numerical integration methods are unstable if the step size used is much greater than the smallest time constant. Dahlquist [l] has introduced the concept of A-stability. The theorem of Dahlquist states that among all linear multistep methods the trapezoidal rule is the most accurate. Gear [2,3] has developed stiffly stable methods based on the necessary requirements of the stiff differential equations. He has obtained stiffly stable methods of order as high as six. Dill [k] has extended this analysis to the methods of order seven and eight. Jain and Srivastava [5] examined a class of methods and obtained stiffly stable methods of order eleven. The purpose of this report is to obtain optimal stiffly stable methods for the stiff differential equation (1) |r= f(t, y), y(a) = n 2. OPTIMAL STIFFLY STABLE METHODS The requirements for the stiffly stable methods are shown in Figure 1. STABLE ACCURATE S T A B L E •0 hX PLANE Fi gure 1 STABILITY AND ACCURACY REGIONS The existence of stiffly stable methods depends on the parameters D, 0, a and on the definition of accuracy. The usual definition of accuracy is that of the order. Dahlquist has shown that if D is zero then the order of the method cannot exceed two. Here we have obtained stiffly stable methods with extended region of stability of order as high as six for suitable parameters D, 0, and a. Third Order Formula : Multistep formulas may be obtained by writing a linear relation between values of the function and derivative at specified points and determining the coefficients in the relation by expansions in Taylor series about an arbitrary point. For instance, we may write ( 2 ) y j.n = a -, y + a o y -, + a -> y « + a i y o + h $o y'^n + T n+l In 2 "n-1 3 n-2 k n-3 n+l n The fourth order formula, that is, the one whose truncation error is of order h , has been obtained by Gear (3) = U8 36 +16 _I + h 12 . y n+l " 25 y n 25 y n-l 25 y n-2 25 y n-3 25 y n+l 12 ,5 (5) "125 h y n If instead of requiring that T , the truncation error, should be of fifth order, we specify a fourth order truncation error we find that there is a one-parameter family of methods. This may be written (h) . ,U8 _ _26_ ) , 36 _ _5T v ,16 ^2 > y n+l " V 25 " 300 Y;y n ^25 300 Y; y n -l ^25 " 300 Y; y n-2 - (— - — v)y + h(— + — y) v' + T v 25 300 r/,y n-3 25 300 w °n+l n or alternatively 1+8 36 16 3 .12 . v = — v - — v + — v - — v + h — v °n+l 25 "n 25 °n-l 25 °n-2 25 ^n-3 25 n+l [26y^ - 5Ty^ n + ^ _ - lly _ - h6y « . ] + T where 300 'n " J n-1 °n-2 J n-3 "'n+l n YhV U) /: m _ i n { n 5 (S) c 12b/y^V125 This later way of writing emphasizes the nature of the integration formula. It is seen to be composed of the fourth order formula with the minus of a third order restrictive condition. Note also that y = — gives the third order stiffly stable method 18 9 2 L , 6 , 3 A (U) y n + l = iT y n " IT y n-l + IT y n-2 + h IT y n + l " 22 h y n The discretization error of the formula {h) is defined as the difference between the value y calculated from (k) and the exact solution y(x ). n n Define e = y - y(x ) n n n Then the error e obeys the difference equation (5) (iii _2£ \ (36 _5T n / 16 k2 > n+1 " ^25 " 300 Y;£ n " ^25 " 300 YJ£ n -l + l 25 " 300 Yj£ n -2 - ( 2T-3^ Y)£ n-3 + hA( i + 34 Y ) £ n + l- T n for the differential equation y' = Ay. This linear difference equation with constant coefficients may be solved by setting e = e, , which gives the characteristic equation (6) pU) - hA cr(5) + T =0 n where and PU; * l 25 300 YH + l 25 300 YJ ^ " l 25 300 YH v 25 300 w ff(€) = ( 2? + 300 Yk It is necessary to bound the solution of the inhomogeneous difference equation (5), Henrici [6], It depends on the stability of the corresponding homogeneous equation which is obtained when T is set to zero. The difference equation is stable if and only if all roots of the polynomial equation (T) pU) - hX c(c) * are inside the unit circle or on the unit circle and simple. If stiff differential equations are to be integrated with large values of h, then large values of hX must not make (7) unstable. Letting hX -> °°, the roots of (7) approach those of a(C) = 0. This implies that the polynomial a(£) must not have roots outside the unit circle and those roots £. for l which |£.| = 1 are simple. This condition is already satisfied by a(£) in this case. For stiff stability we want hX values such that (7) has roots inside the unit circle or on the unit circle and simple. This region is bounded by the locus of p(0/a(0 in the hX-plane for £ = e , 0e[O, 2tt]. For different values of y the locus of p(£)/a(£) is shown in Figure 2. The locus boundaries are closed and symmetrical about the real axis, pass through the origin and steadily decrease in area as y increases from to 7, approximately. At y = 7» the formula ceases to be stiffly stable. The intercept of the boundary on the real axis is given by ( 8 ) ™ - k 192_-_lJx hA " 3 2U + y 36 . y = — gives third order formula which is optimal in terms of least amount of solution history data to be saved at each step of the corrector iteration. Here we find that the values of y in the interval < Y < ■=— give third order formulas which are stiffly stable with reduced truncation error coefficient. We shall denote such formulas by P. Similarly for the values of Y in the interval — < y < 7, we obtain third order stiffly stable formulas which have extended region of stability. We shall denote such formulas by Q. The locus of P(0/a(C), £ = e 1 , 0e[O, 2tt ] for y = 6 is shown in Figure 3. The P and Q, type formulas for y = 2 and y = 6 are listed here. The values of parameters C (error constant), D, and max|?| (£. being the roots of pCO = 0) for formulas of order p ^_ 6 are listed in the table. (P) .262 159 _5jt _L_ 1 h T8 ' y n+l " 150 y n " 150 y n-l 150 y n-2 " 150 y n-3 150 y n+l 12 J n (Q) _ 1U _3_ _2_ _1_ +h A. y n+l " 10 y n " 10 y n-l " 10 y n-2 10 y n-3 10 y n+l 1 h (k) One parameter family of formulas together with P and Q type formulas for order six are given here. Fourth Order Formula : One parameter family of fourth order formula is given by (9) . 300_ 300 200_ _J_5_ _12_ y n+l "" 137 y n " 137 y n-l 137 y n-2 " 137 y n-3 137 y n-4 + h §fK + l ' 3288 [TTy n " 21 K-1 + 23 K-2 ~ 122y n _ 3 + 25y n _ u - h I2y i+1 ] - ^ h 5 y [ 5 } The characteristic equation is <»> C 5 - ^ (300 - '%,)+ ♦ ^ (300 - $,H 3 - ^ (200 - f Y H 2 + I5r (75 - W^ - I5T (12 - 1^ " hx iff (60 + W^ " ° The locus of p(£)/a(0 for different values of y is shown in Figure k. pQO Y = ~ gives fourth order formula obtained by Gear, y >_ 25 does not give stiffly stable methods. The formulas for y = 5 and y = 2k are listed below. (P) 6815 6130 3630 1190 163 y n+l " 3288 y n " 3288 y n-l + 3288 y n-2 " 3288 y n -3 3288 Y n-k l u J^OQ 1 x 1,5 (5) + h 3288 y n + l " 2¥ h y n (Q) . 223 _86_ J}L + J±l. 13 y n+l "" 137 y n 137 y n-l " 137 y n-2 137 y n-3 137 V n-k Fifth Order Formula : One parameter family of fifth order formula is (H) 1 l^r, 522 v 1 ,,„ 1753 N y n + l = IW (36 ° " 725 Y)y n " HPf {k5 ° ~ l2T Y)y n-l 1 n no 25U0 N 1 , ooc 1980 s + 1U7 (U0 ° " ^^ Y)y n-2 " IW (225 - 72T Y)y n-3 1 , VQ 810 . _1_ ( 137 x 1U7 " 720^11-1* " lU7 U " 720 Y;y n-5 6 (6) + h lU7 720 wy n+l 720 The locus of p(£)/a(£) for different values of y is given in Figure 5. 7200 . Y - -, o 7 gives fifth order stiffly stable formula and has the virtue of economy of storage. For y > HO, fifth order stiffly stable formula ceases to exist. The stiffly stable formulas for y = 36 and y = 96 becomes (p) . 6678 72U5 5^6o 2520 y n+l " 29^0 y n 29^0 y n-l 29^0 y n-2 29^0 y n-3 630 63 , 1260 , _i_ h 6(6 ) 29U0 y n-l+ " 29^0 y n-5 29^0 y n+l " 20 y n (Q) , U356 32 Uo , 920 + _585_ y n+l " 2205 n ~ 2205 y n-l 2205 Y n-2 2205 Y n-3 5hO 22k L . 1020 , 2 6 (6) 2205 J n-U 2205 ,x n-5 2205 J n+1 15 J n Sixth Order Formula: (12) y = -1_ ( 29 U0 - ^Y)V - -A- (MlO - ^Y )y y n+l 1089 iy 720 YJy n 1089^ 720 ' yy n -l • 1 (U9 oo-il£yW - 1 (3675 -^Y)y 1089 720 T/J n-2 1089 720 "'n-3 + IbV (1T6U - W* )y n-U " 1559 ik9 ° ~ W^n-5 + im (6 ° - # )y n-6 + h IW (U2 ° + W )y n + 1 Y h 7 (7) " 3W h y n U3200 Y = 1 k 7 gives six order stiffly stable formula as obtained by Gear. Y = gives a seventh order formula which is unstable. The locus of p(£)/a(0 as a function of y is shown in Figure 6. The stiffly stable formulas for y = 2^0 and y = 360 are listed here. The formula is stiffly stable in the interval 15 < Y K 620, approximately. (P) _ 8151 10593 9955 6105 y n+l " 3267 y n " 3267 y n -l ' 326>7~ y n-2 " 32oT 7 n-3 2277 U51 33 1320 , 3267 y n-U " 3267 y n-5 3267 y n-6 " 3267 y n+l 21 J n (Q) . 1737 2061 1685 810 y n+l " 726 y n 72 6 y n-l 726 y n -2 ~ 726 y n-3 171 13 9 , , 300_ , 726 y n-U " 726 y n-5 " 726 y n-6 h 726 y n+l 1 ,7 (7) " IT hy n TABLE OF VALUES OF THE PARAMETERS Order and Type of the Method Y max I £ | 5*1 3 P ■ ' ■ 2 -0.1777 8U°lU' O.U3635 -25/156 Q ' 6 -0.0535 86°50' O.I49 -5/12 it P 5 -1.3953 69°20' 0.62757 -137/1500 Q 2h -0.U238 77°26' 0.60328 -137/360 5 P 36 -3.2166 59°12' 0.7^3 -1U7/1260 Q 96 -I.U966 66°3' 0.71381 -1U7/510 6 p 2U0 -7.1812 U8°36» 0.88UUU -363/3080 Q 360 -5.OIU9 50°3 ) +' 0.8^820 -121/700 10 3. CONCLUSIONS The stiffly stable methods have been found very efficient in integrating stiff differential equations because they allow a much larger step size while maintaining stability. Here we have developed one parameter family of methods suited to integrating stiff differential equations. It is found that stability and accuracy are controlled by the same variable y Ih e methods which are optimal in terms of least amount of solution history data to be saved at each step were obtained by Gear for order as high as six. A suitable value of y leads to methods which have optimal accuracy or optimal stability without violating the condition of economy of storages. There is a great advantage to be gained by the use of a formula with optimal stability when solving sets of simultaneous differential equations whose solution approximate to exponentials with real negative arguments. An optimal stability enables the interval of integration to be increased without causing instability. However, in general the value of y will vary with the behavior of the equations we are attempting to solve. Parametric methods for the solution of ordinary differential equations have also been discussed by Hamming [7] and Robertson [8]. Applications of these and similar formulas to predictor-corrector methods will be discussed separately, 11 LIST OF REFERENCES [l] Dahlquist, G. A Special Stability Criterion for Linear Multi-Step Methods . B.I.T. 3 (1963), pp. 22-1+3. [2] Gear, C. W. Numerical Integration of Stiff Ordinary Differential Equations. DCS. Report No. 221 (January 196?). [3] Gear, C. W. The Automatic Integration of Stiff Ordinary Differential Equations , Proceedings IFIP Congress, Edinburgh (August 1968). [h] Dill, C. A Computer Graphic Technique for Finding Numerical Methods for Ordinary Differential Equations , DCS Report No. 295 (January 1969) . [5] Jain, M. K. and Srivastava, V. K. High Order Stiffly Stable Methods for Ordinary Differential Equations , DCS Report No. 39^ (April 1970). [6] Henrici, P. Discrete Variable Methods in Ordinary Differential Equations . John Wiley and Sons, New York (1962). [7] Hamming, R. W. Stable Predictor-Corrector Methods for Ordinary Differential Equations , Journal of the ACM 6 (1959), pp. 37-^7. [8] Robertson, H. H. Some New Formulas for the Numerical Integration of Ordinary Differential Equations , Proceedings IFIP Congress, Paris (June 1959). 12 Y = Figure 2 Third Order Formula STIFFLY STABLE REGION AS A FUNCTION OF y 13 Figure 3 Third Order Formula LOCUS OF p(5)/a(5), i = e 10 , 0e[0, 2tt] 1U ■2. CD 2.00 6.C0 10.00 14.00 18.00 Figure h Fourth Order Formula STIFFLY STABLE REGION AS A FUNCTION OF y 15 o ■6.GC 24. OC 3G.CC Figure 5 Fifth Order Formula STIFFLY STABLE REGION AS A FUNCTION OF y 16 y = 36o \y = \y = 2U0 1+3200 1U7 12. OG 18.00 24.00 30.00 Figure 6 Sixth Order Formula STIFFLY STABLE REGION AS A FUNCTION OF y IT 1 . 00 Fi gure 7 SECTION OF LOCUS OF pU)/o(£) FOR p < 6 Form AEC-427 (6/68) AECM 3201 U.S. ATOMIC ENERGY COMMISSION UNIVERSITY -TYPE CONTRACTOR'S RECOMMENDATION FOR DISPOSITION OF SCIENTIFIC AND TECHNICAL DOCUMENT l Sm Instructions on Reverse Side ) 1. AEC REPORT NO. COO-1U69-016U 2. title OPTIMAL STIFFLY STABLE METHODS FOR ORDINARY DIFFERENTIAL EQUATIONS 3. TYPE OF DOCUMENT (Check one): Kl a. Scientific and technical report I I 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) 4. RECOMMENDED ANNOUNCEMENT AND DISTRIBUTION (Check one): [Xl a. AEC's normal announcement and distribution procedures may be followed. I I b. Make available only within AEC and to AEC contractors and other U.S. Government agencies and their contractors. 3 c. Make no announcement or distribution. 5. REASON FOR RECOMMENDED RESTRICTIONS: 6. SUBMITTED BY: NAME AND POSITION (Please print or type) C. W. Gear, Professor and Principal Investigator Organization Department of Computer Science University of Illinois Urbana, Illinois 6l801 Signature ^h^ tjo^ Date June 1970 FOR AEC USE ONLY 7. AEC CONTRACT ADMINISTRATOR'S COMMENTS, IF ANY, ON ABOVE ANNOUNCEMENT AND DISTRIBUTION RECOMMENDATION: 8. PATENT CLEARANCE: □ a. AEC patent clearance has been granted by responsible AEC patent group. I~~1 b. Report has been sent to responsible AEC patent group for clearance. I 1 c. Patent clearance not required. ^^ \#\