fMHBBHMHBJlfl BH m IB MBBlHnHBffiliS 1 If H H Hi ;vV HB9I Bran B I ■ i IH Nn8B hbHI BBSS IB Hi Hi R3H Bnl LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 510. 84- IJl6r cop.2 j^j , j. UIUCDCS-R-TU-672 TrusCZ, coo-2383-0012 y MULTI-DERIVATIVE NUMERICAL METHODS FOR THE SOLUTION OF STIFF ORDINARY DIFFERENTIAL EQUATIONS BY ROY LEONARD BROWN, JR. August, 197^ DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN URBANA, ILLINOIS ^■Hf IJBRARY OF THE FEB 24 1975 UNIVERSITY OF ILLINOIS UIUCDCS-R-7^-672 COO-2383-0012 MULTI -DERIVATIVE NUMERICAL METHODS FOR THE SOLUTION OF STIFF ORDINARY DIFFERENTIAL EQUATIONS BY ROY LEONARD BROWN, JR. August, 197^ DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URB ANA -CHAMPAIGN URBANA, ILLINOIS 6l801 * Supported in part by the Atomic Energy Commission under contracts US AEC AT( 11-1) 2383 and US AEC AT(ll-l)lU69 and submitted in partial fulfillment of the requirements of the Graduate College for the degree of Doctor of Philosophy in Computer Science. Digitized by the Internet Archive in 2013 http://archive.org/details/multiderivativen672brow 11 MULTI -DERIVATIVE NUMERICAL METHODS FOR THE SOLUTION OF STIFF ORDINARY DIFFERENTIAL EQUATIONS Roy Leonard Brown, Jr., Ph.D. Department of Computer Science University of Illinois at Urbana-Champaign , 19TU Current research in the physical sciences and engineering has greatly enhanced the importance of computer programs for the numerical solution of systems of ordinary differential equations. Recent developments have concentrated on the solution of systems characterized as stiff for which the requirement of numerical stability can cause some integration methods to use a smaller stepsize than the control of the local truncation error requires. The stability and convergence of linear multi-step multi-derivative formulas have been investigated since it has been shown that linear multi-step formulas of order greater than two cannot be A-stable, a useful attribute for formulas used to integrate stiff systems. It has been shown that when a stable and consistent linear k-step s-derivative formula implemented in an iterated corrector method is applied to a system satisfying certain continuity requirements, the solution is convergent . A set of 11 multi-derivative multi-step formulas, all A-stable and of orders 1 through 11, have been found and implemented in a variable stepsize, variable order method as a FORTRAN program called T)k . Tests of DU on a set of stiff problems show it compares favorably with another stiff integrator. Ill ACKNOWLEDGMENT I wish to express my gratitude to my advisor, Professor C. W. Gear, for his sure guidance, constructive criticism, and valuable support during the development of this work. D. S. Watanabe, R. D. Skeel, and the other members of the weekly numerical analysis seminar contributed much by listening to several presentations on the related research. The ALTRAN compiler and its software support were made available by the Math and Numerical Services Unit of the Computing Services Office, University of Illinois. A. C. Norman provided the TAYLOR compiler and execution supervisor, as well as valuable suggestions. Computing time was provided by the U. S. Atomic Energy Commission under grants 1U69 and 2383. Special thanks are due Pam Farr for her efficient typing of the draft and final versions of the manuscript. This work is dedicated to the memory of the mathematician C. L. Dodgson (1832-1898). IV TABLE OF CONTENTS Page 1. Introduction: Stiffness and Stability 1 1.1 Stiff Equations 2 1.2 A-Stability k 2. Mult i -Derivative Formulas 9 2.1 A Survey of Higher Derivative Methods 9 2.2 Some Formulas Well -Suited to Stiff Problems 12 2.3 Convergence and Error Estimation for Multi -derivative Formulas 18 3. Implementation 35 3.1 The Iterated Corrector Method 35 3.2 Stepsize and Order Changing Hi 3.3 Program Description hk 3.H Properties of Stepsize and Order Changing Methods ^5 3.5 Higher Derivative Calculation 53 3.6 Numerical Tests 56 k. Alternatives and Conclusions 59 U.l Matrix- Valued Coefficient Methods 59 k.2 Conclusions 66 LIST OF REFERENCES 69 APPENDICES I DU FOR THE INTEGRATION OF STIFF SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS 72 II TEST PROBLEMS FOR COMPARISON OF DU WITH DIFSUB 86 1 . Introduction: Stiffness and. Stability Current research in the physical sciences and engineering has greatly enhanced the importance of computer programs for the numerical solution of systems of ordinary differential equations. The usage of both special purpose and general programs has increased, and, in particular, the availability of methods for solving systems characterized as stiff has led to an increase in attempts to solve such systems. For examples of both stiff problems and current methods used to solve them, see Bjurel et_ al [3]. This work will present a method that should solve a large subclass of all stiff systems; however, the amount and type of information it requires, which includes the second or higher derivatives of a first order system, place it as a special purpose program suitable for certain problems which more general computer packages are unable to solve efficiently. The method will be developed in several steps. First, the underlying formulas on which the method is based will be derived and shown to meet certain previously defined desirable criteria. Second, the program package design will be considered, with justification offered for the design choices made. Third, the results of tests comparing the program with a similar package having a similar program organization will be presented and discussed. The desirable criteria of A-stability, convergence, and stability at infinity are presented in Chapter 1 following a discussion of stiffness in the next section. In Chapter 2, conditions necessary for A-stability and sufficient for stability at infinity are derived and used to search out a set of formulas that satisfy these desirable criteria. A theorem on the con- vergence of these formulas is proved, and an asymptotic error estimate is derived which provides formulas used by the method to calculate stepsize and select order. Some observations are made about a recent conjecture of Genin [2U] concerning the maximum order attainable by mult i- derivative A-stable formulas . Chapter 3 presents some details of the program implementation with notes on program usage. The effect of changing stepsize and order is considered. Test results are presented. Chapter k presents a more easily used matrix-valued coefficient formula with stability properties similar to those of A-stable multi-derivative methods. Concluding observations on the usefulness of both multi- derivative methods and matrix- valued coefficient methods are made. 1.1 Stiff Equations The solution of problems in such areas as control theory, chemical kinetics, circuit theory, and nuclear reactor kinetics requires modeling of physical systems by large systems of ordinary differential equations. If only Lipschitz problems, as defined below, are considered, then a unique solution exists. For a standard proof, see Coppel [9 : P- 12]. Definition 1.1 A Lipschitz problem is the initial value problem y_' = f(y_, t), (i.i) y_(t ) = Zq, such that f.(y_,t) is continuous in t and satisfies the Lipschitz condition II Id.;,,, t) - f(£g S t) || * l || ^ - ^ || for some positive L in the region t <: t < t , y_ < «. A Lipschitz problem is considered stiff if its numerical solution has certain stability properties that make the computation difficult. For example, consider the system iL' = A y_, (1.2) where A is a real n x n matrix with its eigenvalues {X.} satisfying Re (X.) < 0, at least two solutions are not identically zero, and | Re (X.) | o = max IM is large in the sense that it is greater than some threshold value a . Ehle [15] classifies a range of stiff problems with stiffness ratios o varying from mildly stiff for 10 $ a $ 100 to extremely stiff for a > 10 Equations (1.2) are stable since the true solution 6 y_(t) = e (t_t )A y^ (1.3) approaches zero exponentially as t increases. In fact, the non-zero component corresponding to the largest | Re(X.) becomes small enough that its contribution to the truncation error of the numerical solution becomes negligible before any other component. Yet for standard techniques such as the Adams method (see Bashforth and Adams [2] and Moulton [3^]), and explicit Runge-Kutta formulas (see Ralston [37: 199]), absolute stability, to be considered shortly, requires that max I hX. < C (l.U) where h is the steps ize in the independent variable and C is some function dependent on the numerical method. Therefore, a method suitable for stiff equations is one which places no upper limit on | hX. | where X. is an eigenvalue of the locally linearized system v_' = f y_, for Re(X.) < 0, other than restrictions based on controlling the error of the computation. r\ -f* Here f = _j^ > the Jacobian matrix of f_ evaluated at y_ and t. ~ L dv_ Ehle [15] has proposed that a system be called stiff only with respect to a method. Thus , a problem p is step-wise-stiff with respect to method m if the average stepsize required in solving p with m is smaller than that used in solving p with some other method m , other parameters being equal. Thus, finding appropriate stiff methods is equivalent to a search for maximal elements of the set {m. } with respect to step-wise-stiffness for all p in some class of problems P. Further readings in stiff systems can be found in Bjurel et_ al [3], Ehle [15] » and in the works listed in Appendix II of Enright [l6]. 1.2 A-stability Stiff equations were first discussed as a problem in numerical analysis by Curtiss and Hirschfelder [10], who encountered the problem in chemical kinetic research. Dahlquist [11, 12] first provided a general treatment of the stability of multistep formulas, defined below. Definition 1.2 A k-step multistep formula for problem (l.l) has the form = .1 a {0) y . + a U) f(y . ,t .) (1.5) i=0 1 n-i 1 n-i n-i where y approximates y(t ) and t = t_ + nh. n n n Dahlquist [13] subsequently investigated the special properties necessary for stiff multistep methods. Gear [l8: 11 6-132] expanded this work to include his multivalue methods which treat directly the elements of the vector used to store past information. For multistep methods this vector is T 2L. = ( y > • • • »y i »hy ' , . . . ,hy ' ) , 1 4i n 'n-k n n-k ' where T is the transpose operator. In a multivalue method, this vector may be transformed to store differences, Taylor series, etc. Multistep methods are a subset of multivalue methods, so the latter will be discussed here. The main construct in the theoretical treatment is the stability region of the test equation y' = Ay, (1.6) y(t Q ) 4 o. Definition 1.2 The stability region R associated with a multivalue formula is the set R = {hX : the formula applied to (1.6) with constant stepsize h > produces a sequence {y } satisfying y -*■ 0}. When a method is applied to a system of equations y_ f = Ay_, it is absolutely stable for values of hX. I e R where X. are the eigenvalues of the 1 i ■ i Jacobian matrix A. The most desirable property for a multivalue formula suitable to stiff problems is that j hX. | never be restricted by the requirement of absolute stability for Re(X.) < 0. This property is defined as follows . Definition 1.3 A multistep formula is A-stable if its associated stability region contains the open left complex half plane. Dahlquist [13] also proved that the highest order attainable by an A-stable multistep method (1.5) is two i.e. the local truncation error V* (t n» = fa "I ' ^ (t n-i' + *f f ^ * + h '4 lJ f WVi>'Vi> (!•»' = o( h P +1 ) then the formula has order p. If p > 1 , the formula is consistent . Consistency and the strict root condition are sufficient conditions for a fixed step multistep formula to be convergent i.e. y converges to y(t ) as n -*■ °° for stepsize h = (t -t-)/n. *- n f One final requirement precludes methods with the undesirable property that lim n -»■ °° n y n-l -*■ 1 as hA -> -oo where y is the calculated solution at t_ + nh to (1.6). n Definition 1.6 A formula is stable at infinity if there exists a real W < such that sup Re(hA) < W lim n -*■ °° 'n-l < 1, (1.10) where y is the multivalue method solution to the test equation n at t + nh for fixed h. 2. Multi-Derivative Formulas There are many examples in the literature of formulas that require derivatives higher than the first to solve first order systems of ordinary differential equations. Some authors provide justifications for requiring higher derivatives hy demonstrating techniques for calculating them; others do not. Throughout the remainder of this work, a k-step s-derivative multi value method will be defined by k = A ( 4 0> *n-i + ill 4 J> hJ ^ y n-l> (2 ' 2) which can be of order 2s and A-stable but are then not stable at infinity. Order 2s is achieved when aj: = (-1) cc are the coefficients of the jth diagonal Pade approximation to e (see Graves-Morris [25]). Other one-step 10 methods, all based on Pade approximations, are derived in Blue & Gummel [h], Cavendish et al [8], and Varga [^0] . Since these methods are of order at most 2s a great deal of work must be done in comparison with other multi -derivative methods, since they attain lower order than the conjectured maximum order of an A-stable multi- derivative method, as presented in Section 2.3. Liniger and Willoughby [32] consider the two parameter class of formulas : = -y n + y n _ x + | [(l-a) y r J_ 1 + (l+a)y^] (2.3) + h 2 [(b-a) y" - (b+a) y"]. k n ~ 1 n If — $ a+b $2 and $ b - a $ — , the formula is A-stable. By picking a and b appropriately the numerical solution can be made to follow certain components of the linearized problem Enright [l6] considers formulas °=^n + ^-l + Jo "I"' f( Vi' t „-i )t 4 2> ' ,(y B'V- {Z - k) \ This provides formulas that are A-stable through order h and stiffly- stable through order 9« By requiring the equation to be autonomous of the form y_' = f(y_) ■ (2.5) Z(t ) = I<) , J (easily accomplished by adding a new variable defined by T 1 = 1, T ("t n ) = t , to the system) , the second derivative can be evaluated by rca^i = ty&ij' (2 - 6) 11 Since most methods for stiff problems require a Newton type corrector iteration to provide rapid convergence for large stepsizes, some method for finding f , the Jacob ian matrix, is usually available. However, an approximation to the Jacobian is usually sufficient. Therefore, many programs do not evaluate f exactly. Further, the procedure is time con- suming and so a particular computed value is often carried over several steps. To require exact evaluation of f at every corrector iteration is perhaps excessive. See Brown [5] for similar second derivative formulas giving higher order A-stable methods. Jeltsch [30] has developed techniques for studying mult i -derivative methods . He calls p(C) + | x u J a U) = 0, (2.7) J=-L J V - h\, the characteristic equation of a multivalue multi-derivative method (2.1), for p(^) as in 1.8, and a.U) = J 4 J) C*- 1 . (2.8) These equations are used to define the concept of e -stability. Definition 2.1 If the roots £ (h\) of the characteristic equation corresponding to the multivalue method satisfy for some finite k |hX| n max {|£.(hX)|} < k i for all n <: e as hX ■* -°°, then the method is e -stable . A sufficient but not necessary condition for a method to be stable at infinity is that the method be e -stable for e > 0. In fact, lim n -> 00 y n y n-l 12 as hX ■+ - 00 if a method is e-stable for e > 0. Jeltsch shows that e may be calculated geometrically from the Puiseux diagram corresponding to the method. Definition 2.2 The Puiseux diagram for a method (2.1) is the 2-dimensional point set P = {(i,j) : a£ ± 4 0). Puiseux diagrams for methods of Gear, Enright, and Jeltsch are shown in Figure 3. Jeltsch shows that e is the smallest slope of all straight lines having a finite slope and passing through |k, max (j) \ and any I ^4 other point of P. 2.2 Some Formulas Well-Suited to Stiff Problems A set of formulas that will solve any linear and most non-linear stiff problems should satisfy several criteria: 1) stability at infinity; 2) A-stability; 3) convergence with as high an order as possible to ensure as large a stepsize as possible and thus reduce the global error. With this in mind the following definitions and lemmas are useful. Definition 2.3 A multi-derivative multivalue formula (2.1) is implicit if at (« least one of aJ. J 4 0, j £ 1. Definition 2.k A mult i- derivative formula (2.1) is maximally implicit if 13 4 1 ? 4-1 (0) — * — * — * — *- 1 2 3 h (A) (2) 2-i * (1) i 2-i * * (o) 2-i " 1 2 (B) 2 2-i * (1) 2-i * # (0) 2-i "^ 1 2 (c) Figure 3. Puiseux diagrams of Uth order methods due to Gear (A), Enright (B) , and Jeltsch (c). «<*> * o. Lemma 2 . 1 If the stability region of a multi -derivative formula contains hX = -oo, then the formula is maximally implicit. Ik Proof Consider the characteristic equation = P(c) = pU) + Z y j ff (c) (2.9) J •*- J (s) and suppose a* = 0. Let the roots of P(c) he £ (y ) , . . . ,C ( y ) . P(c) can he expressed as pU) = S In ^ J) y j .J. (c - c.(y)). (2.10) j=0 i~l i It will be shown that requiring lim £.(y) to he bounded leads to a k— £ contradiction. Consider the terms in c; in each of the expansions (2.10) and (2.9). Setting them equal gives A a o' } y ' i -°°, obtaining (s) - n a £ - 0. (s) (s) Thus if a =0, then a £ = 0, I = 1,2 , . . . ,k, and the method cannot be an s-derivative method if hX = -°° is in its stability region. Therefore for a method's stability region to contain hX = -°°, it must be maximally implicit. Q.E.D. Since an A-s table formula must have -°° in its stability region, the requirement of maximal implicitness is a necessary condition for A-stability, Lemma 2 . 2 (s) If a method is maximally implicit and a. =0 for i f 0, then it is e-stable with e > 0. 15 Proof It is evident from the form of the Puiseux diagram, see Figure 3, (s) (si that E > if a: } and a. ' = if i / 0. Q.E.D. With this information, several restrictions are imposed on the formulas which will be considered for inclusion in a multi value method of form (2.1). First, choose a. = for i > 0, and a). S ' ± 0. Since (s ) a i- requires that the s-th derivative of y be calculated, it will probably be worthwhile to calculate the 1-st through (s-l)-th derivatives as well, so a J ^0, j=l,...,s. A further heuristic limitation involves the belief that the most recent information, i.e. closest with respect to the independent variable to t , is the most useful. Thus, if (2) (2) (2) a method is to use three values for y' ', then a , a , a^ will be non-zero . A program was written to graph the stability region boundary of a method described by ( s ,m ) where s is the highest derivative used, and m. is the number of (most recent) evaluations of y , j = 0,...,s. J This program computes the coefficients a. of the highest order method possible for a value (s,m), and then plots the locus of y for which the characteristic equation (2.7) is satisfied for any C = e , £ 6 £ 1. The following set of formulas were picked because they appear to meet the criteria of A-stability and stability at infinity, and were also of the same general form which facilitates including them in a variable order method: ■ Jo "I"' *n-l + * ^ K k " 1 ' 2 (2 ' 12) If. • - JU 0) ^ ♦ A ** #> *' K = u - 5 The formulas are of orders 1 through 11, the order being k+s-1, with e = s/k for e-stability. The stability region generated by the graphical program for several of the above formulas appears in Figure U. Since it might be possible that the boundary of the stability region crosses the imaginary axis and not be apparent on the plot, the roots |C . (y)| for jj = iy with y e{0.01n: n = 0(1)200} U {2+0. In: n = 0(l)l80} were calculated and found to be less than 1. Thus, at least a numerical indication exists that the stability region doesn't cross the imaginary axis and that the formulas (2.12) are A-stable. Genin [2k] has conjectured that the maximum- order attainable for A-stable multivalue methods satisfies = Ub^-1 s * k ( j *max V 1 3s-l s £ k. The cases (s=l, k=2) , (s=2, k=U) and (s=U, k=8) attain this maximum order. A comment on his constructive demonstration for several values of s and k of this conjecture appears in section 2.3. It should be noted that when s > 1, these orders are higher than those attainable with the Obrechkoff's methods. The coefficients a for each method are rational numbers computed i by replacing y _. by y (t _. ) in each of formulas (2.12), expanding them in a power series about t , and equating coefficients of like powers of h for as high an order as possible. This leads to a linear system of k+s+1 equations for the unknowns a. . The matrix associated with this linear system, as described in section 3.1, is a variant of the Vandermonde IT 10 in CM in i m ■ C.(0 1X10 a. 7s I Figure h. Stability regions for formulas (2.12) of order p for p equal to 3, 6, and 11. Formulas are stable on the exterior of the closed curve. 18 matrix and has a non-zero determinant, so the system has a solution. The search program described above computes the a. as double precision floating point numbers. Computation of the coefficients as rational numbers is described in section 3.1. 2. 3 Convergence and Error Estimation for Mult i -derivative Formulas Although the stability characteristics and consistency of the formulas (2.12) have been established, the standard result that stability, in the sense that p(§) satisfies the strict root condition, and consistency are together sufficient conditions to guarantee convergence has not been established for methods of the form (2.1). This section will provide such a result. Further, an asymptotic error estimate will be provided for the global error e n = y n - y(t n ) (2.H0 for small stepsize h. The theorem statement yields the proper scaling for the local truncation error d n " y „ " y n-l (t n» (2 - 15) where y n _ 1 ( t ) solves y^_ 1 (t) = f (v n _ 1 ^ t) > ^-i^n-l^ = y n-l* A note ° n the asymptotic error behavior of some formulas of Genin (197*0 used to motivate his conjecture about the maximum order of A-stable mult i derivative methods completes this chapter. All of the formulas (2.12) are implicit, thus requiring the solution of a system of non-linear equations. If this system is solved iteratively then the method of the numerical solution will be called an Iterated Corrector method. Since any iterative procedure needs some initial estimate near the fixed point of the system, the way this estimate v , » to y is obtained is the method Predictor. The predictor can be as simple as 19 ^,(0) ^n-1 or can be the result of extrapolation from the previous computed values. Its exact form has no effect on the solution as long as y • , x is in the radius of convergence of the iterative process and the corrector is iterated to convergence. If a bounded number of successive substitution corrections are used, the predictor affects the numerical solution through the truncation error. Let an iterated corrector method be applied to the single equation y' = f(y,t). Let the values used by the method be stored in a vector y^ = (y n ,hy ] J....,h y n ,. . . ,y n _ k , . . . ,h y n _ k ) where the vector elements are numbered from and T is the transpose operator. Then a successive substitution iterated corrector method with an associated predictor may be formulated as: Predictor: jr^^, = Bjr^ > (2.16) Corrector: F(^ >(m) ) = a< 0) y n>(m) + J^ \*t ( ^ (y n>(m) ) + l[a<°> y . + la'JW^ty )]. 1=1 i n-i j=l l n-i s 1 ( i-l) ^,(m+l) ^(m) -0 VJL n,(m)' j=l -J n,(m+l), f f (i) (y / n) - f (i) (y / ,0 m > o / . n VJ n,(m) VJ n,(m-l) n,(m+lj A /.. \ t -i-1 f UJ ( Y , x) - (y / rt \). -,h m = I ^ y n,(0) ; ^,(0)^+1 Here y , x is the m-th corrector iterate of y^, and y^ /x = (y^ ( m ) ) Q = T T £ y , . for £^ = (l,0,...,0), coefficients are normalized so that 20 ol = -l, and the argument t has been omitted from f(y,t) when its value is clear. If the predictor B uses extrapolation, then its first s + 1 rows contain a generalized Hermite extrapolation operator for y ~_ ,j = 0,...,s or the identity matrix otherwise, and the shift operator in the remaining rows. The i-th unit vector is I . . When iterated to — l convergence, this formulation is equivalent to any other formulation that solves the corrector exactly and will be used in proving theorem 2.1. Theorem 2.1 Given 1) a consistent p-th order mult i derivative formula (2.1) for which , r> £ (0) k-i • satisfies the strict root condition; 2) y' = f(y,t) where f(y,t) is continuous for bounded y(t), t £t£t„, and each of f^ _1 '(y,t ) , j = 1,. .. ,s, exists and is bounded by L , y n J and f (y,t ) is continuous in y, for mesh points t = t + nh, n = 1,.. . ,N; 3) a starting error e^ = y_^ - y_(t n ) bounded by i i i i n' llejl .< D»h P ; then an iterated corrector method which is iterated to convergence (by any method) at each point t or which is iterated by (2.l6) a bounded number of times (in which case, the truncation error of the predictor- 21 corrector formula must have p-th order) is convergent and there exist constants C , C , C 2 , D, h Q dependent only on the method and f(y,t) such that I le,, | I * C^'e^ + a o D h Pe C 2 b for any < h < h , b ^ Nh = t - t , provided that the solution y(t)e C _ , the space of functions with continuous (p+l)-th derivatives. Proof Express the local truncation error d as follows. Let y(t ) be the -n *- n correct value of y , and let (2.l6) be applied to v_(t ) as JL.(o)" ,B ^Vi ) ' (2 - 17) ^n = ^(M ) ' n where M is either bounded or is large enough that the corrector has converged. Then the local truncation error is Here Af / ., \ is defined for f (y / \ ). Because the method has p-th n,(m+l) •hijCm) order (if M is bounded, the predictor-corrector has p-th order), d can n -n be expressed in a Taylor series with remainder terms in h times com- binations of derivatives of f and y which can be bounded by the continuity hypothesis. Hence I Id || * Dh P+1 . 1 ' — n ' ' Let ^ = ^ -2L(t n ): 22 and ^n,(m) " ^n,(m) " y n,(m) ' From the definitions Sn,(0) = t,(0) " t,(0) = B Sa-1- = ^n " ^ U n } = ^ " ^ + ^ " L<* B > e — n ~ Sn,(M ) + 4* n By subtracting (2. IT) from (2.l6) and using the mean value theorem, defining f (y) = y and omitting the argument of 3f/3y = f for notational convenience, ^Wo =[1 + jo'o^'^.w'o (2 - l8) and h ~^Wj ■ #U) - ^U> ■ ^""'^.w'o- (2 - 19) Note that the arguments of f may not be the same in different appear- ances of f , but all we will need is that | If I i is bounded. ( 0) T Therefore, since a^ = -1 and (e . ). = %, , n x. e , ,. -n-i — (s+l)i -n,(m) e = S e + d (2.20) -n n — n-1 -n where M S = it? [I + .Lh¥ H) a.i! (2.21) n m=l j=0 y — j-0 K ( i) t — 1=0 l — (s+lji I is the identity matrix with rows through s zeroed, and SL_ is the j-th unit vector. 23 First consider the case of bounded M . Since If S = S^ + h S n n (J-D| < L , then J where S is a matrix whose elements are polynomials in h with bounded coefficients. The constant matrix s o = 8 a<°> (0) M 1 o i \ where a. is the (s+l)x(s+l) array with a. in the (0,0) position as the only nonzero element, B is the (s+l)x(k+l) (s+l) predictor operator, and represents the (s+l)x(s+l) zero matrix. Due to the form of these operators, for any M >, 1, s o = /;<°) ... ;<°> s\ (2.22) 0/ which is the companion matrix of a polynomial whose only non-zero roots are the roots of p(5). Hence, 1 is the only eigenvalue of S with norm not less than 1, and there exists a constant C such that for all j >, I S J I 6 0. 11 M v If M , the number of iterations at each t , is not bounded allowing n n iteration to convergence, it can be shown that the coefficients of h in S are still bounded for small enough h. n 2k Let oj be the total correction applied to y , _» to get y , that ^n " ^n ^n,(0) Since the iteration has converged, we must have F(y ti ) = and h J f (j-D (( } j = ( ) . . m ± B (2>23) is 'Vo J = ^ J Similarly, let oo be the change in y , x to get y. y also satisfies equation (2.23). By the mean value theorem = F(^) = F(^) ♦ Tfa - ij and -n — n h j f; j '% - ^) - (^ - ju, 3-i.... We can use equation (2.l6) to write these as = - E n - 2 n a J) h J f (j-D £ T / __ _ j i=0 j=0 l y -(s+l)i v,i 4i,(0) ^,(0) =ia ^ (2.2U) and ° = 4 (s n.(o) -i.to)*^-^' (2 - 25) j = l,...,s Note that all but the first s+1 components of to and w are zero so that — n — n l_ (w - to ) = for j > s. Premultiply the equation (2.2*0 by -l^ and the equation (2.25) by £_. 9 j = 1 9 . . . 9 s ; add, and add the equations J {. E £.£ T + K, -L .Z_ h J f^ _l) v>? 'lj A _ v.}(a) - w ) = to get j>s —j— j -0 j=0 i=l y i is+l)i -n -n 25 = (1 1 vlJWJ- 1 ^ ,,. 'i=0 j=0 -0 i y ^-(s+l)i j=i -3-3 j=i -j-o y ^,(0) ^,(0)' The matrix on the left-hand side can be inverted for all sufficiently small h because of the Lips chit z conditions on f , hence there exists an h such that for h s h the inverse of that matrix can be written as I + hQ where Q is bounded. Since the right-hand side contains polynomials in h with bounded coefficients, we can write k rp (0 ) s T A a) - w = (.Z„ l n l, ,\.a. - .£., il.il, + hQ)(y , n . - y , >) -n -n i=0 -CP{s+l)i l j=l -j-J ^ /VJL n,(0) **n,(or where Q has bounded coefficients. Hence -n ^n *-n *n *- n = J ^,(0) -*n,(0) + ^-% + 4 S rp k rp ( rv\ = [(I - .2, SL.l. + .E_ i I 1 . xl v.a u; )B + hQB]e _■ + d j=l — j— j i=0 -0— (s+l)i l -n-1 -n = S e , + d n-n-1 — n where S = S,+ hS , S is bounded for h < h„ , and S. is identical to the n n n previous S given in equation (2.22). Then for h < h , there exist C , C^, such that ||S J || * C Q , 3 5 1, (2.26) l|S n l| * C i: (2.27) These satisfy the hypotheses of lemma 2.3 below, so s n | m 1 3 R ...S || * c n e h(n - m)C C l. (2.28) n-1 n-2 m 1 ' 26 For t N = t fl Thus , % " I s - +1 A, + *" V %lhN I 1^1 I ||d.|| ♦ ||s![ +1 || I 1^1 I ( 2 .30) S C De bC 2h P * C D'e bC 2 h P ' Th e constant terms are b = Nh , H = (t -t )/h, C = C (t -t ) , C = C C . The error | |<= | | goes to zero as h -> , I |e«| I ">■ 0, and the formula is convergent. Q.E.D. Lemma 2 . 3 If S = S + hS and there exist positive constants C_, C, , and h_ n n 10 such that for h < h , 1) ||s°|| S C Q , j >. 1, 2) ||i B || t c r then iS n || = ||S n S ...S II .< C n e (n - m)hC C l. (2.31) m 1 ' ' ' n-1 n-2 m 1 ' Proof It is true for n = m. Assume it is true for m <: j < n. Then ||s n || * ||s n - m || +h. ?^||s n -J|| ||S. J I MSJ" 1 !!. (2.32) Note that hC Q C 1 s ( e hC ° Cl - 1). (2.33) Replacing norms by bounds in (2.32) and using (2.33) yields I |s n | I * C n + . ? n (e hC C l - l)e (j-m-l)hC C 1> (2>3M 1 ' m' ' J=m+1 27 Thus, S°|| < c o e (n - m)hC C l, (2.35) and (2.3l) is true by induction. Q.E.D. There are several points to note about Theorem 2.1. First, the hypotheses are more strict than those for the convergence proof of 1-derivative multistep methods. In particular, the existence of f , j = l,...,s, and the continuity in y of f (y,t ) are required to assure the existence of S . n Second, the bounds shown for the global error are not useful in a computational sense because C and thus C are dependent on knowing the bounds L on f throughout the region of integration. Also, error bounds developed by replacing summations of N terms by N times the largest term, as in (2.30), tend to be very pessimistic. An error estimate accurate for small h is more useful when trying to control the error in a computer calculation. Theorem 2.2 provides such an estimate. The following lemma will be used in proving Theorem 2.2. Lemma 2 . k If = .£ n [a< 0) y . + A lA.W f^T 15 ] i=0 L 1 J n-i j=l 1 n-i r- is a convergent method of order r 5 1, and p(c) satisfies the strict root condition, then = ija! 0) y . + ha U) f .] (2.36) i=0 1 n-i 1 n-i is a convergent method of order r 1 5 1. Proof Consider 28 y^ n » ■ jo^w + J l 1 hJ 4 J)f(J " 1, (^Vi ) 'Vi )1 (2.37) for y(t) = e . Then L h (e Xt ) - .| [(a(°» + .t^^he^-^) = e i(t - hk) [p( e hX ) + 1 (hX) 3 a (e hX )]. (2.38) Since the method is of order r L h (e Xt ) = C r+1 (hX) r+1 e X(t - hk) + 0(h r+2 ), Hence p(e hA ) + J ± (hX)J a d (e hX ) = C r+1 (hX) r+1 ♦ 0(h r+2 ). Writing hX = log (l+z) and using hX = z + 0(z ), we see that (2.39) p(l+z) + .| (log(l+z)) j a. (l+z) = C z r+1 + 0(z r+2 ) (2. hO) cJ J X J- is a necessary and sufficient condition for the method to be order r. Expand this in a power series about 1 in z; p(l) + zp'(l) + za(l) + 0(z 2 ) = C A . z r+1 + 0(z r+2 ) (2.1+1) 1 r+1 so order r^l= > p(l) = and p'(l) + a, (l) = 0. Thus, a method involving only p(C) and a -,(£) will be consistent. If p(C) satisfies the strict root condition, such a method will also be strongly stable and thus convergent . So ik < a i <0) Vi + h *i (1) f n-i> < 2 "^ defines a convergent method of order r'^l. Q.E.D. We can now state the theorem describing the behavior of the global error for a small but non-zero stepsize. The result will also yield the 29 proper scaling for the discretization error L (y(t )) to give the local truncation error (2.15), where d n=^n - y n-l (t n»' Theorem 2.2 If a strongly stable convergent pth order multivalue formula of form (2.1) is applied to y T = f(y,t), y(t ) = y Q , and f (y,t) is continuously differentiable with respect to t, and f (y,t) , j = 0,...,s-l, is continuous in y for bounded y in the range of t, and all starting values are exact, then e n = h? 6(t n ) + 0(hP+1 )' where 6(t) satisfies Q 6' = f (y,t) 6(t) + g +1 , , y (p+l) (t), y I a. i=0 1 6(t Q ) = 0. Proof Consider the discretization error which, from the assumption on f(p), ■ Vl hP+ly- The corrector iterated to convergence satisfies 0= .i[a (0) y . + 1 h j a^f ( J" l} ]. (2.UU) i=0 i J n-i j=l i n-i Subtracting (2.U3) from (2.UU) gives is 30 .l[a (0) e , + 1 h J ■o} i) (fi 4 7 l) - t^hyit .),t ,)]♦ (2.U5) i=0 i n-i j=l i n-i n-i n-i p+i ( P+ D (t j . 0(hP+2 p+1 n Assume the theorem is true for e n through e , . The mean value theorem n-1 gives f (j-D _ f (j-l) ( (t )t ) = f (j-D ((t )+e t )e (2>1|6) n-i n-i n-i y n-i n-i n-i for some 6 satisfying $ j 8 ! <: i e . . Since by theorem 2.1 iii n _ 1 i e n _. -0( h P ), (2.H6) can be written as f (J 7 D _ f (j-D ( (t _) t <)=f (j-D ((t )t )e . +0 (h 2p ). n-i n-i n-i y n-i n-i n-i (2.U7) Thus, .! n [a (0) e . + ha U) f(y(t . ) ,t ,)e . ] ♦ C +n ^ +1 y (p+1) (t ) 1=0 i n-i i y n-i n-i n-i p+1 n = .Z n .L h- j a^ ) f ( J- l) (y(t .),t .)+0(h P+2 ). (2.U8) i=0 j=2 l y n-i n-i If e = h P 6(t ) + 0(h P+1 ), n n then l[a (0) 6 . + haf l) (f (y(t . ) ,t . )6 . + -f^l- y (p+l) (t .))] i=0 i n-i i y n-i n-i n-i f (1) n-i i=0 a i = 0(h 2 ). (2.1*9) From lemma 2.k, this is a convergent method which solves 6'(t) = f (y(t),t)6(t) + g +1 (i) y (P+l) (t). (2.50) ,Z_ a. i=0 l Q.E.D. 31 Under the more realistic assumptions that the starting values are not exact, a similar result can be shown, following Henrici [26: sec 5.3]. Because of (2.50) the appropriate quantity to estimate and control is d = h? +1 V 1 y ( P +l) (t). (2.51) k (1) iio a i This will be discussed further in the next chapter in relation to an actual computer program. Genin [2k] has conjectured that the maximum order attainable by an A-stable k-step s-derivative method is f 2s + k - 1 k < s Pmax = < , n v ■ ' 2 -52» 1 3s - 1 k $ s . To motivate this conjecture, he develops canonical multi derivative formulas for the linear multistep (s = l) and linear multi derivative (k = l) formulas and demonstrates that p is attainable as in (2.52) for k-step s-derivative max r formulas that are products of these canonical formulas. The analysis depends on an algebraic manipulation of the multi derivative multistep formulas to yield a canonical polynomial H(p,q) in the variables p = (£+l)/(£-l) and q = hX corresponding to the characteristic polynomial p(£) + 1 | 1 (hA) 1 a.U) = 0. He proves that a formula (2.1) is A-stable if and only if the canonical polynomial H(p,q) = .j jZq a p k_J q 1 = (2.53) is a two variable Hurwitz polynomial in the narrow sense. He also provides an order characterization based on H(p,q). If H (p,q) is the canonical polynomial of the unique optimal highest order A-stable linear multi- 32 derivative formula, then for all k < s, the canonical polynomial \ iS (p,q) ■ H s _ (k _ l) (p,q)[H 1 (p,q)] k - 1 k < s is an A-stable k-step s-derivative formula of order 2s+k-l. If G (p,q) is an A-stable 1-derivative k-step formula of order 2, the maximum attainable order for multistep formulas, then G , (p,q) = H (p,q) is unique. For all k $ s Vs = G k-( S -l) ( 5-' l)[H l (p > 1. (hX) 1 " 1 C t t . V a (0 > "l: 1 ^ e Xt ^-J-i n i=l v+r+i m=0 m j=0 Since the difference operator for the solution e satisfies r ( e At n) = [p( e h> ) + .§ (hX) J c .(e hX ) ]e At n-k = C v+r+1 (hXr r+1 e Xt *-* + 0(h v+r+2 ), then .^[(hX)" 3 c.(e hX )]/-p(e hA ) = 1 - C v+r+1 (hX) V+r+1 (2.55) p(e ) .-#. v+r+2 , / hX\ \ +0(h /p(e )) Expanding e in a Maclaurin series and using it in p(e ) gives p(e ) = a Q (hX) E where E is an infinite sum in hX with a constant first term. Hence from (2.5*0 3U 1+ " c ^ ^n(hX) V e tn - K (hX) V At n v+r+1 n , v+l\ y n - e n - + (h ) a o z which is the expected form of the global error of a method of order r lower than what is expected from the truncation error. Thus, theorem 2.1 does not apply for the above formulas which attain p ; they appear to achieve that order only in the indicated truncation error, are not stable under the usual definition, and are very sensitive to error. In fact, the apparent error for all these formulas when integrating y' = e is 2s + k - 1 - (k-l) k < s p app ^3s - 1 - (s-1) ■ k >, s 2s , the same as for Obrechkoff's methods. Genin [2k] shows that for (k=2,s=2) (k=3,s=2), and (k=2,s=3) there are no methods with order greater than p IuclX by applying the Hurwitz criterion to H (p,q.) . However, all three of the formulas given that achieve the expected maximum order have pU) = (I - i) 2 pU) and therefore have a global error of a formula of order 2s for the equation y' = Ae . It should be noted, however, that p , as in (2.52), is attained ZlLcLa by formulas (2.12) for (k=2,s=l), (k=U,s=2), and (k=8,s=U). This suggests that the conjecture is possibly correct; only the supporting formulas in [2k] are inapplicable. 35 3. Implementation A computer program for a variable step, variable order method for the numerical solution of ordinary differential equations is based on four elements: 1) A set of formulas - to compute the next approximation, and descriptions of how such formulas are to be applied; 2) An estimator - of the local error committed at each step; 3) An acceptance rule - to determine whether to accept a value from a single step or try again with a smaller stepsize; k) A stepsize strategy ~ to determine the best next stepsize as well as how often to change it. In a variable order method, the best stepsize is picked from that used by several different orders, so the next order determination is part of the stepsize strategy. Hull [28] and Ehle [15] consider these attributes of different methods in relation to specific problems in order to prove effectiveness and compare numerical behavior, respectively. This chapter will define a method for solution of stiff ordinary differential equations by providing these four elements. To avoid notational complexity, all references will be to the equation y' = f(y,t) rather than to a system of equations . Where evaluation of a system will require greater complexity, such as for the Jacobian f , which is a, matrix for systems , such complexity will be noted. 3.1. The Iterated Corrector Method All of the formulas (2.12) are maximally implicit which makes them 36 suitable as corrector formulas in an iterated corrector method, as defined in section 2.2. By adding an extrapolating predictor formula of the same order and choosing an efficient vector representation a for the saved values of the dependent variable, an iterated corrector formulation is provided as the method calculator. A form for representing the dependent variables proposed by Nordsieck [35] provides a convenient extrapolative predictor as well as an efficient means for manipulation of the stored values. (l) s (s) The set of saved values hy ,. . . , h y , y . , i = 1, . . . ,k, used n ''n n-i to compute y in the k-step s-derivative formulas (2.12) uniquely determines a k+s-1 degree polynomial which interpolates these values. This polynomial can be equally well represented by its value y and its first k+s-1 derivatives at t . The Nordsieck method stores the elements of the Taylor n series of this polynomial . about t . If 2n = ( V y n-l' ••" y n-k+l> ^J^ > "" h \ M ^ and where q = k+s-1, T is the transpose operator, and y is the 0-th element of both a and y , then I h = Q Zn (3.1) where Q is a non -singular matrix. Further, a .. = Pa + 0(h q+1 ) (3.2) -n+1 -n where P is the Pascal triangle matrix {P..} = {(f)}. 37 Thus, (3.2) can be used as a q-th order predictor since i ! J! The efficiency of storing the dependent variables in such a form arises from two considerations. First, the predictor step can be cal- culated using only additions of elements of a in the same way as the Pascal triangle matrix can be formed. Second, when a new stepsize has been chosen, the new variable can be formed by /1 a = -n a — n for a = h/h. This essentially performs a new q-th order expansion about t which is equivalent to an approximation to a new set of data n ~ ( l) ~s (s) y , hy , ..., h y , y , ..., y . Although y(t -jh) may never II il II I1~_L Il™iv II have been approximated by some y ., it is still within 0(h ) of the correct value by the properties of the Taylor series. As will be discussed in section 3.^, this does not affect the stability of the method if stepsize and order changes are restricted by certain criteria. Define the change in y necessary at each corrector iteration by 'hi (m) n,(m) j=l w n,(m) n k (0) Z= i Z =l a i y n-i' y n,(m) = ( £n,(m)V (3 ' 3) The subscript in (y ) . denotes the i-th element of y (the initial 38 element being (y ) = y ) . The iterated corrector formula is then *n,(0) = B *n-l- «<*> f(i, ^n,( B )» " f(i> ^„, („-!)» ">° n,(m+l) \ ,.\ r^.CO)) " i * = ° • B is the matrix Q jPQ. The subscript n denotes y = v_(t ) , and (m) denotes the m-th corrector iterate. The vector e. is the i-th unit —l vector, where e^ = (1,0,...,0) T . This can be converted to Nordsieck's form as a ,_x = P a = QBQ" 1 a _ . -n,(0) -n-1 -n-1 .(J-D -n,(m+l) -11, (m) -0 v -n,(m)' j=l j n,(m+l) ,-1 G(a ) = F(Q x a , J -n,(m) -n,(m) *0 = V ij = Q( a Q j) Sq + £j)> j=l,...,s. (3.5) The derivation of the coefficients a. and I. can now be described more fully. To get a. for a formula (2.12) requires the solution of the rational system 39 -2 (-D 2 (-2) 2 (-k) 21 -k 2 (-D q (-2) (-k)« (3.6) ,(j) .(0) after scaling the a. so that a = -1. I is the s x s identity matrix. For finding the H. coefficients, note that — 1 -1 ... (-D ,-1 1 (-k+1) ... (-k+l) q 1: ... 2: ... ... is also rational. The vectors £_. , i = 1,. .. ,s, can be found by solving tf-^-c^ao + Si)- A program was written in ALTRAN , a formula manipulation language from Bell Laboratories described in [T]» to set up these rational matrices and vectors and solve for the coefficients using an algebraic linear equation solver ASOLVE which is a resident language routine. The results for (iOai 1 *, i = l,...,s, and I. appear in the comments accompanying the program DU in Appendix I for formulas (2.12) up through order 8. The rational results for larger orders could not be computed by the ALTRAN program due to overflow in the rational representation of data. Uo Another point of interest is the method for calculating E, a linear combination of previous computed values y . . If this computation were done with the exact values stored in y , then at every change in stepsize either other calculated values of y . must be available, thus restricting n-i the program to integral multiples on increasing stepsize and starting over at the lowest order on decreasing stepsize as in [l6], or else interpolating new points . The stored values of the dependent variables form the Taylor expansion of an interpolating polynomial P(t) which, for all t in t $ t £ t , satisfies = } n a? 0) P(t-ih) + L ai J> h J P (t3) (t), 1=0 l j=l p(t n -ih) = y n _ i , i = l,...,k+l, p (J) (tn _ 1) = f (J-l) (yn _ lsvl) 9 3 = 1 ,..., s . Hence the formula 1 ~ i=l a i y n-i " j=l a h y n,(0) + y n,(0) (3 " T) is exact. Formula (3.7) is used in calculating Z immediately after the predictor step in DU. Although the restriction on hX due to the absolute stability require- ment on the formula has been removed for Re(hX) < 0, the corrector may still not converge for the large values of h which the program may attempt to use. This can be overcome by considering a Newton's iteration (or quasi- Newton) on the 0-th term of the vector a . The function whose zero is -n sought is F(2L) = (-2L) + * + J-L h J a^ ) f (j - l) ((y_) , t) Ui where (y_) is the unknown y , Z and t remain constant, and the f^ -1 ^ are considered as function of y . The corrector iteration on y / \ = n n,(m) ^.(m^O then beC ° meS y n,(m+l) =y n,(m) " ^"^.(m)* (3.8) = y , x + [I - ! h J a^ } f (j-l) (y , v ,t )] _1 F(y , .). n,(m) j=l y w n,(m)' n Vi n,(m) This leads to the iterated corrector formula with extrapolative predictor a n = P a _ (3.9) -n,0 -n-1 U7 ' .(3-D ,(m+l) -n,(m) -0 -n,(m)' j=l -j n,(m+l) -n, W = [I - ,L h j a^V^Cy , v,t )] _1 . J=l y n,(m) n Because of the expense of evaluating W, which is a matrix for solving systems of equations, it is evaluated as seldom as possible, since when the resulting quasi-Newton method (3-9) converges, it converges to a solution of (3.5). In fact, W can be computed by numerical differencing with respect to y_ and still give the desired convergence. Since the proof of convergence in Theorem 2.1 assumed only that the corrector was solved, the fact that Newton's method is used does not affect the proof. 3.2 Steps ize and Order Changing Let 5 .. = C .. q+1 q+1 where C . and a\ are associated with the formula of order q in (2.12). q+1 Then the truncation error per step in the asymptotic error estimate of Theorem 2.2, k2 is used in the stepsize and order changing algorithm. An estimate E of d will be made and this will be compared to the desired error EPS set by n the user to determine whether or not the step will be accepted. This can also be used to determine what stepsize h should have been used instead q. of h to get d = EPS for the qth order method used whenever it is possible q n e or necessary to change the stepsize. In practice, h will be chosen smaller than the optimum value computed to account for the error in E . n To derive this algorithm, note that the predictor step of (3-9) does not modify the last element of the (q+l) length vector a , and thus its final corrector component r? X u.) Af^; 1 ! m=l J=l -J q. n,(m) is a backward difference h 4 U) h a (a) h q+i -S ^- — ^ +1, *o(^). (3.12, If, after the corrector converges and the error estimator satisfies E = C ._ q! || .L (A.) (h j f (j_l) - e. T P a •) || < e (3.13) n q+l " j =1 -j q -77 n j -n-1 ' ' J * for e = EPS, then the desired stepsize for a method at order q approximately satisfies h d> = e, q n and the current stepsize h nearly satisfies U3 h q+1 * = E , q n n where d> is the principal error function C y . n c q+r'n Therefore , \ = \(f] 1/q+1 - (3.1*) Thus, the asymptotic error estimate is used to decide if a value of a will be accepted, and to compute the best stepsize at order q. For — n possible order changing, the best stepsize is computed at convenient orders, in this case, q+1 and q-1, and the order corresponding to the largest stepsize is used. The computations for one order lower is E ^ ■ K »' (2 nV (3 - 15) h = h Je \ 1/4 . An estimate for h y can be obtained by differencing E to get E : « n * = $& (E Q - E^) = 8^ ^ T M ♦ 0(h^3), (3.16) Vl By dividing the calculated values h by 1.2 and h , , h _ , by 1.3, q q-1 q+1 the computed next value takes into account the inaccuracy of the approxi- mation to the asymptotic error estimate and also projects a bias toward leaving the order at q, since changing order has some overhead associated with it. Also, no stepsize increase of less than .lh is allowed because of the overhead involved. When the estimated error is too large, only the best stepsizes at orders q and q-1 are considered. uu The four elements of a variable stepsize, variable order multivalue method have now been determined. The calculator is (3.9). The estimator is (3.13). The acceptance rule (3.13) depends on the estimate of d by E . The stepsize strategy combines the stepsize and order calculations (3.lU), (3.15) 5 and (3.l6). The exact algorithm involved when changing order will be considered in section 3.H. The variable representation is the Nordsieck form a which interacts with the above error estimators to -n provide access to the higher derivatives needed, to predict a / n (, and to provide an easy means of changing a with changing stepsize. 3. 3 Program Description A program J)k has been written to incorporate formulas (2.12) and the previously presented representations and algorithms into a FORTRAN subroutine package. The organization of the program and particularly the form of the stepsize and order changing procedures are very similar to the program DIFSUB of Gear [21]. The heuristic devices employed therein are used extensively throughout D^, which should result in approximately the same proportion of program overhead in the various types of calls to DU, i.e., successful steps, steps that don't converge on the corrector iteration, and steps that converge but don't meet the error criterion. Some major features of the program are outlined below and the program is listed with comments on its use in Appendix I. The user is responsible for providing the differential equations in the form of a subroutine DIFFUN(T,H,Y,DYS,NORD) which places the first WORD elements, excluding the 0-th, of the Taylor series for Y about T into DYS. A routine is provided to calculate W by U5 numerical differencing, but the user may provide a routine to calculate W exactly if desired. In the linear case when f = A 1 for yj = Ay the exact calculation will take less execution time than the N calls to DIFFUN for an N-dimensional system needed for numerical differencing; in the non-linear case, exact f may improve the rate of convergence. Each call to DU takes one step in the independent variable T and returns. The input value of stepsize H is tried first and reduced if necessary; the output value is the best stepsize for the next step. This may be changed by the calling program but an increase in H will probably not be accepted. The program keeps internal variables from one step to the next and this prevents the use of T)k to solve more than one problem using alternating calls . 3.U Properties of Stepsize and Order Changing Methods All of the theory on stability and convergence presented in Chapter 2 and sections 3.1 and 3.2 were developed on the assumption that the order remains constant throughout the calculation and the stepsize is such that the values y . are available as calculated. It is possible that the n-i ^ interpolatory change in a when the stepsize is changed might introduce factors which will prevent the method from being stable and convergent. Wo algorithm was presented for changing the order, and factors involved in changing the order which also might disturb the stability and convergence of the total method need to be considered. This section will consider these points and develop the order changing scheme used in T)h. A concise definition of the stepsize scheme containing the stepsize and order selection algorithms will depend on several definitions, given U6 below. Definition 3.1 A step selection s cheme is a function 9 such that h - h 9(t ,h) (3.17) n n where, for some h > 0, t £ t <: t. , 1 >, 9(t,h) 5 A >' 0. Definition 3.2 A formula selection scheme is a function l(h,t) such that the formula used for the interval [t , t ■_] is F T/n \ where h is the n n+1 I(h,t ) n parameter of definition 3.1. In the following, the formulas F. are the 11 formulas (2.12) and F. is the formula of order i. Thus, l(h,t ) is an order selection i n scheme for T)h. Definition 3.3 A method is stable with respect to a step selection scheme 9_ and a formula selection scheme I_ if there exists a constant M < °°, dependent only on the differential equation, such that y - y I < M Iv - v I (3 18) |J m "^m 1 |J n °n ' ^.j.o; for all t_ $ t < t £ t„, where y. and y. are two numerical solutions n m f l i to the same problem using the same variable step, variable order method. Definition 3.^ A set of formulas {F. } is convergent with respect to 6_ and I_ if the computed solution y converges to y(t ) for any t $ t < t as h -*■ and the starting errors -> 0. UT Gear and Tu [22] and Skeel [39] investigate the effect of varying the stepsize on the stability and convergence of normal form methods, i.e. methods which store the dependent variables in Nordsieck form, and derive stepsize changing criteria which result in convergence of such methods with respect to stepsize changing schemes 9 which meet such criteria. Gear and Tu also consider variable step Adams methods that store y = [y , hy',...,hy' ] and use stepsize dependent coefficients a. , a. . These methods are compared with the normal form methods i,n i ,n and are shown to have less stringent requirements on the stepsize selection scheme. Gear and Watanabe [23] consider the effect of varying the order of normal form methods that are stable and convergent with respect to some stepsize scheme 6. Criteria for the order selection scheme I are established which result in stability and convergence of normal form methods with respect to I. The basic construct of these investigations is the stability matrix S of section 2.3 based on the matrix representation of the iterated n corrector formula (3.5) with the extrapolating predictor P. For normal form methods , this is S = t? [I - .1. U. e. T - h f (J_1) I. e n T )]P, (3.19) n m=l j=l -j -j n y -j -0 where f are the Jacobian matrices of the derivatives and P is Pascal's y triangle matrix. The stepsize changing algorithm can be represented as C = diag LI, a, . . . ,a J for a = h /h , when changing from the stepsize h , to the new stepsize n n-1 o d jr- n _2 U8 h . A number of possibilities exist for changing the order. Shampine n and Gordon (1973) have calculated the exact procedure when using the Nordsieck representation of the dependent variables to achieve the same result as when a q-th order Adams method using the representation y^ = (y^, hy n ',..., hy' n _ k ) changes from order q to either q-1 or q+1. Brown [6] develops this for Gear's backward differentiation formulas. However, these are quite complicated and it happens that the following simpler device will suffice. When the order is lowered or remains the same, the value (a ) is dropped -n q or left, respectively: when the order is raised, a new value (a ) _, is -n q+1 calculated by differencing ((a ) - (a n ) )/(q+l). Alternatively, this n q n-1 q term could be left alone and the initial estimate would be (a ) , n = 0. -n q+1 If 0. .is the matrix representation of the order changing algorithm for 1 s J going from order i to j, then = I (3.20) = i* = I - ( e e 1 ) q,q-l -q_ -q. q,q+l \ (0 n ) row } row q+1 column q y (cl) 1 /., ln-2, where A = — — ( 1 - —> — rj . n q+1 (q)' y n-l Alternatively , h9 - I^ 1 , q,q+l The method can then "be described by If we define a = S C a _ , (3.21) -n n n n -n-1 v J ' a = T a ,. -n n -n-1 S = S + h S , (3.22) n q n n u ' K ■ (I - k h ^ T)Mnp ' where S corresponds to the order q formula in (2.12), then the eigenvalues of T = S C (3.23) n q n n can be used to determine if the method will be stable and convergent since if the eigenvalues of T exceed 1 in norm then even for small h the n n computation will not converge. These matrices have been used to establish criteria for stepsize changes which allow a predictor-corrector method of form (3.9) to be convergent and stable with respect to the stepsize changing scheme 9(t , h) where n h = 6(t , h) h, n n h ^ max (h. ) , 1 S i $ N for N the number of steps taken in integrating from t fl to t f . The relevant theorems from Skeel [39] are presented in Lemma 3.1 If a variable step method is such that 50 1) The underlying formula F is strongly stable; 2) F has order a £ 1: q 3) There exist positive constants A and V such that < A $ 6 < 1 n and e, I e - e . $ V; n=l ' n n-1 ' then the method is stable and convergent with respect to 6. Gear and Watanabe sought a criterion to limit the frequency of order changes that would allow a multivalue method to be stable and convergent with respect to a formula selection scheme. The result of interest may be expressed as Lemma 3 . 2 If a variable step, variable formula method is such that 1) Each formula F is strongly stable; 2) Each formula has order q $ 1; 3) The order changing operations are exact for y' =0, and i ! 0. . II is bounded; h) The coefficients in the method M corresponding to F are q q bounded; 5) Each formula is stable and convergent with respect to its stepsize changing scheme; then, for each pair of formulas F. , F., there exists an integer p such i J -'-J that the method is stable and convergent with respect to a formula selection scheme that produces infrequent formula changes in the sense that at least p. . steps are taken without formula change after formula 51 F. has been selected before formula F. is selected, i J The values of p., of interest in DU are listed in Table 3.1. i x^ 1 2 3 h 5 6 7 8 9 10 11 1 1 2 2 1 1 3 3 1 1 2 h 1 1 2 5 1 1 U 6 1 1 2 7 1 1 3 8 1 1 9 9 1 1 k 10 1 1 3 11 Table 3.1. p. . for changing order in DU. There is now enough information about the behavior of variable step, variable order normal form methods to make a statement about their stability and behavior with respect to changing stepsize. In particular, the method M , as presented in section 3.2 and this section can be shown to have desirable stability properties. This can be seen by considering the use of IVL, on a problem y' = f(y,t), y(t ) = y , t a $ t $ t f ' ,(r) where f ' (y,t) is continuously different i able with respect to t and f (y 5 "t), j = l,...,s-l, is continuous in y for t <. t $ t and bounded 52 y. In particular, if r is the highest order formula used, s = 1 for r <: 2, s = 2 for r $ 5, s = 3 for r <■ 1 , and s = k for r $ 11. The re- straints on f(y,t) satisfy the hypotheses of theorem 2.2 so the asymptotic error estimate will be applicable. M . also requires that a minimum allowable stepsize H . be chosen to keep the error calculation from mm being affected by machine roundoff, and that a maximum useable stepsize H no greater than the interval of integration be imposed. Hull and max Enright [29] discuss the automatic choice of these values for different integrators. However, in DU the choice is left to the user. The method implemented in the program DU meets the hypotheses of Lemma 3.1 with the parameters h $ H max h = 9(t ,h)h, n n A = H min^' N*N=(t f -t )/H m . n , N I I V $ Z n ! 1 - A i n=l ' ' Further, the hypotheses of Lemma 3.2 are met if the order is never increased after less than q+1 steps, for q the current order, since p _, < q+1 for all values in Table 3.1. This is true in the program DU. q.,4+1 Hence, the method M , is stable with respect to the stepsize and order changing schemes. Since a lower bound is placed on H, the stepsize, the method cannot converge to the true solution as h, or H , is forced uLcuC to zero. However, if greater accuracy is requested on one machine inte- gration than on another, greater accuracy can be expected. Note that the method is not guaranteed to integrate to y(t ); the corrector might not 53 converge or the calculated error per step may not be small enough for some h = H . . Also, if y(t ) is computed, the global error has not been estimated, although at each step the relative error per step was near EPS, a user input. 3.5 Higher Derivative Calculation The most significant drawback to using DU is the requirement that derivatives of the first order equation solution be available through order h. A permanent disadvantage is that the calculation of each added derivative will increase the computation time as least linearly with S, probably more, for any function of interest. Further, the average user does not want to have to program these added derivatives or even calculate what they are, since this can be quite time consuming and tedious. Several means of alleviating this second problem are available. If the system of equations is autonomous and the exact Jacobian f (y ) is available through a subroutine call, then y n f (y) = f y (y n ) f(y n ) C3.2U) which can be used for up to fifth order A-stable formulas in DU. If the Jacobian is easy to calculate because the system is not highly cross coupled, or the Jacobian is well known, or if the system has a linear dependence on y as in y_' = Ay_+ f(t), (3.25) and f(t) is easily differentiable, then a similar device can be used. Enright [l6] provided (3.2U) as the justification for his use of second derivative methods that are A-stable to *rth order. Also, for linear and simply coupled autonomous systems it may be easy to calculate f = f f + f 2 , y yy y 5k which is used in the convergence matrix W. If the function f can be given as an algebraic expression of a limited number of functions, then various algebraic manipulation languages can be used to compute the necessary derivatives; f , f , ..., may also be computed. By using I/O facilities in these languages an entire FORTRAN sub- routine DIFFUN could be produced. FORMAC, a PL/I precompiler described by Fike I 17], has a function DERIV(F,T,N) that finds d F where F is an algebraic function of T and a^ any of (SIN, COS, SINH, COSH, ATAN, ATANH, LOG, LOG10, L0G2 , ERF). ALTRAN is a FORTRAN based language that handles ratios of polynomials, as well as real numbers, integers, and rational numbers. The functions DIFFN(F,T) in ALTRAN serve the same function as DERIV in FORMAC. The function TPS(F, T, N) calculates the Truncated Power Series CTaylor series) of F up to order N . A common method for solving differential equations is to calculate the Taylor series at a point (y , t ) and then find (y ,_ , t ,_) using n n n+1 n+1 specified error criteria to choose h _ = t , - t . One program that * n+1 n+1 n does this is called TAYLOR, which is written in BCPL and which accepts a description of the system involving expressions in any of the FORTRAN functions and outputs two FORTRAN subroutines . One of these output subroutines (the name is given by the user) accepts initial values and calculates the first Taylor series; the other finds a value at a given t by taking as many steps as necessary. Two additional routines, XTAY01 and XTAY02, are also produced and are called 55 by the first two. The method used was coded by Norman [36], who describes the program TAYLOR, and Barton et al [l] derived the underlying formulas. The elements of the Taylor series with h = 1 are calculated using recurrence relations based on the simplest rules of differentiation. For example, if y' = y-L y 2 > j 2 = y 2 > then M/u = Ci-2)U y (i - i y i - 2) +y < 1 - 2 > y < i - 1 >)/i, can be solved for i = 2,...,s. See Barton et al [l] for a more detailed discussion. Since the TAYLOR system was used to provide the up to Uth derivative calculator DIFFUN in the test of J)k in section 3.6, the special use of TAYLOR for that purpose is described below for a simple example. The input program for y' = y would be as in Figure 5. COMMENT Y' = Y TERMS h INITIAL SET (T, Y, TERMS) ADVANCE VAL (T, Y) EQUATIONS Y» = Y END Figure 5. Input to TAYLOR to solve y 1 = y. 56 The output would be a FORTRAN program SET (T, Y, S) which, given T and Y(T), produces the 0-th through s-th terms of the Taylor series and calls XTAY02 while doing so. SET is changed slightly to become DIFFUN (T,H,Y,DYS,S). VAL would be used to find Y(T) by the Taylor method and TAYLOR will not compile unless some minimal routine is required. TERMS k limits the terms in the series to through h and saves storage by doing so. 3.6 Numerical Tests Ehle [15] has tested a number of FORTRAN programs designed for solving stiff systems of ordinary differential equations using three values of Z on each of 2h small systems (N £ 9). Each problem is specified by p = 0^0 f max mm for y* = f(t,y), y(t ) = y Q and YMAX- is the initial value of y for a relative error calculation in max which the truncation error per step is calculated relative to max(y , y ) , ■ max n t is the final value of t , S is a set of values of t at which the global error is to be calculated, and the code specifies various parameters of the problem, such as a stiffness ratio, the linearity of the problem, and the range of a such that all eigenvalues are in S(a) as defined in Chapter 1. This last is the most important for T)h since in general none of the methods in Ehle's tests did as well for a $ 30° as for a < 30°. Since DU is never restricted by the eigenvalues as long as they have negative real parts , tests were made which show that DU did as well for 60° $ a > 30° and 90° 5 a > 60° as for 30° £ a > 0°, in the sense that no serious convergence problems or large errors occurred during the calculations. The global error 57 at points in S were monitored for those tests having closed solutions, and these were usually within a factor of 10 of EPS, the requested error per step. A subset of Ehle's problems was chosen to include as many problems with a > 30° as possible as well as some simpler benchmark tests, and these were compiled using TAYLOR. Some tests were eliminated that were inappropriate, e.g., one linear constant coefficient problem had very large values in its corresponding matrix, and since TAYLOR subroutines calculate the derivatives without scaling by H the higher derivatives produced overflow errors. The number of calls to DIFFUN by T)k and the number of matrix inversions for W were measured. W was always calculated by numerical differencing but this use of DIFFUN was not included in the recorded number of calls by DU. Results are in Table 3.2 and the test problems are listed in Appendix II. NFNS is the number of calls to DIFFUN, and NW the number of matrix inversions. 58 Dk DIFSUB TEST eps NFNS NW NFNS NW A a = 0° -3 10 -6 io_ 9 10 285 775 3186 91 207 807 122 263 380 12 . 19 17 B a = 0° 79 188 U33 12 25 Ho 12 1* 2kk 333 12 17 17 C 30° < a < 60° 86 U 210 311 76 k9 56 . 93 229 2H9 11 29 19 D a = 0° 126 3U6 120 k 22 60 202 150 35U U05 1U 30 22 E 30° < a < 60° 136 329 1U3U 23 U3 185 173 351 538 18 27 21 F 60° < a U90 1709 5631 3k 38 U8 1752 2181 3213 9 16 11 G 60° < a 375 1251 3115 23 33 U7 3577 5VT1 535fc 18 6 Ik H 60° < a 1737 5621 30 life 57 lll»7 2233 3006 8 11 11 I a = U5° 37 28 25 10 6 6 52 U8 38 9 6 5 J • a = U5° 26 26 26 12 12 12 21 21 21 2 2 2 K 60° < a 322 318 727 1+7 111 52 U26 U21 9 31 25 L 60° < a 75 75 61 25 2^ 2k 6l 71 68 20 22 2k M 60° < a Ul2 lU2 659 56 26 103 kl 108 U79 9 18 36 Table 3.2. Comparison of DU and DIFSUB on 13 tests from Ehle [15 ] 59 k. Alternatives and Conclusions The value of multi -derivative methods lies in their applicability to a wide variety of stiff problems, many of which cannot be properly solved by the less expensive to use first derivative methods such as Runge-Kutta, Adams, and backward differentiation. By using A-stable higher derivative methods, as large a stepsize as is consistent with the error criterion can be taken at any order given suitable start up and order changing procedures In this chapter, methods that require only the first derivative yet retain some of the properties of the higher derivative methods are presented. The final section considers the usefulness of the methods presented throughout. k.l Matrix-Valued Coefficient Methods When higher derivative formulas cannot be used, methods which share some of their desirable properties can be substituted; in particular, if we consider that (f y ) j y = (^"V - y (j) > (^D and in the linear constant coefficient case, y' = Ay, (i+.l) is exact, then formulas can be developed which use only y and y ' but have coefficients that include powers of an approximation Q to the Jacobian f (y ,t ). Since the approximations (U.l) are not exact for non-linear problems even if the calculated f is exact, additional non-zero coefficients must be y added to insure that the order is not lowered. Since these new non-zero coefficients may destroy the property of A-stability of a formula, new non-zero coefficients will have to be added to restore A-stability or at least stiff stability. Lambert and Sigurdsson [31] developed formulas k = 1 [a< 0) I + .!. a. (j) hJ oj] y . + (U.2) i=0 l j=l l Ti n-j * i hi\ l0)l + e i h\ i3)hi « 3]t »-y 6o By forming a Taylor expansion around y(t ), a method of a given order p can be obtained by matching terms in h Q J y J for i = 0,...,p, for each j = 0,...,s. This leads to a linear system of p+l-j equations for each value j £ , . . . ,s . If Q =X for the test equation n Y' = Ay, then the stability polynomial becomes p(£) + ^ (hX) J a (C) = (U.3) for (h.k) 0.(0=^ (a^.bp-^)^ 1 , p(5) = .So a. 5 . If the additional coefficients available as a. , b. , j > 0, l i beyond the p+l-j required to establish the order of the method are used to improve the stability of the method, then the extra matrix computations required by a method based on a formula (U.l) are justified. The stability region depends on terms (hX) J a. and J where some arbitrary £. will be called the "root of interest" and the £ , J * k £ = 2,...,k are "extraneous roots" of o . . Define ir = Z_ ( £ - ? ) for further use. If the extraneous roots of o. were made the same as the J extraneous roots of p(0 then for these roots | £ | < 1 for all values of hX in the stability polynomial (U.3). Thus, the stability region depends only on the principal root 1 of p and the roots of interest of the a.. The above can be accomplished by setting 6l i! (a i J) + ^ J_1) ) ^ k_i = o, (U.6) q = 2,...,k, for £ a root of p(5), ? 1 = 1. This adds an extra k-1 constraints to a linear system to be solved for coefficients a )b (J>, 1-0....* for each value j = l,...,s. Lambert and Sigurdsson call a formula that satisfies (k.6) stabilized . The form of formula (U.l) requires that b ' = for all i, so the i minimum value of s for a stabilized formula of a given order p can be (s) found by requiring that the k+1 coefficients a. , i = 0,...,k, have no more than k+1 constraints on them. There are p+l-s constraints to achieve (s-l) order p, and k-1 stabilizing constraints if a. ^ 0, i - 0,...,k. Thus, p + k-s<:k + l; s 5 p - 1. Since up to 2k+2 values a. , b. , can be non-zero, these values can be picked to produce formulas (U.2) that are A-stable for Q = f . In particular, if the coefficients are chosen in terms of parameters such that the order and stabilizing conditions are met for all values of the para- meters, then the parameter values for "which the formula is A-stable for Q = f can then be determined and, within those constraints, particular values may be chosen. One may choose the parameters to minimize the truncation error coefficients for accurate approximations Q to f , or to minimize the dependence of the truncation errors on Q , or to optimize the error in some other way. 62 An example of this from Lambert and Sigurdsson [31] is the third order formula = [-1 + (^a)^ - aQ^]y n + [(l+a) + ~<-15+a+3-26a(l+a) ).^ (U.7) -bQ 2 ]y . + [-(a+3) + kk-&+6&(a+Z)Q + (3a+2b)Q 2 ]y _ + n n-i z n n n— d [3 + ^|{-5-a+33-2Ua3)Q n + (2a+b)Q 2 ]y n _ 3 + [^|(23-5a-3) + (^a(l-a)+b)Q h ]f n _ 1 + [^|{2+a-3) - (ga + a(3-a+3) + 2b)Q n ]f n _ 2 + [j|<5+a+53) + (|e+a(2-3 ) +b )Qj f n _ 3 > for a = i 2 + % 3 , 3 = ? 2 C 3 , U. | < 1, i = 2,3, a * j|. The truncation error is h l4 [ 2 ^(9+a+3)y (i+) +.[^{1^+23) + ^a(lT+a-3) ]Q n y ( 3) + (3a+b)Q 2 y (2) ] + 0(h 5 ). The stability region for (U.7) for the parameter values a = 3 = 0, a = 1/6, b = -1/2 is plotted in Figure 6. The extraneous roots of the characteristic sub polynomials p(^), a ( E,) , and a {E,) for the above method are the E, , E, used in choosing a and 3. For E, - E, = 0, the characteristic polynomial becomes 5 2 [(s-D + iaU+i/2) + h 2 x 2 c] = o. (U.8) Note that if the approximation Q to f is not exact then the n y stability region R(e) for Q = (A+e), y' = Ay, is all hA which satisfy the root condition for 63 in P in rvj to CM in i m r- i TXIO 0.(0 a. 7s m i Figure 6. Stability region when Q = X for matrix valued coefficient formula (U.7). 6k p(C) + ha 1 (c) + h 2 a 2 U) = 0, (I1.9) where k a.U) ^[ap^X+e)* 3 + bp _l) ( X+e ) j_1 X ]£ k_i . This is the same as Then for it defined as for (U.5), (^-.9) becomes ir[(e-l) - h£ (j^ 4 0) £ k-i ) + hQjU-^) - he (Jo bP^" 1 )] IT 7T + hV (5-? P )] =0. (U.10) n d- In the example above, £ = -1/2, and £ = . For h and e small enough and ^ near the unit circle, the he terms are minor perturbations to all TT but the high order term of a quadratic equation in hX. The perturbed stability region R(e) is thus very close to the unperturbed region R, and converges to it as e -»■ . If a formula of form (^4.2) were used in an iterated corrector method with a quasi-Newton corrector, the corrector calculation would have the form = . h[ | r i (t (o) +b (i) v)( ^ ((m+i) .^ j(m)) , where G is the difference function (k.2). Thus, if Q^ = f afi * -I + (JU + b<°>)hQ + (a n 2) + b n l} )h 2 Q 2 . (U.12) — n n dy Since the inverse of dG is already needed, premultip lying dG by dy dy 65 before inverting yields *».(*« " *,.<») m " h(b o 0) + " b o 1> S. ) «4.(»n)-^.(-) ) ■ (U - 13) where W is the inverse of a linear combination of hQ terms up to the n * fourth power. By modifying an existing program such as DIFSUB to compute W instead of W, and premultiplying Af , > by h(b;, + hb: Q ), these n, \m) On methods can be easily implemented. Since the local truncation error is known in terms of y , Qy , 2 (2) and Q y this can be computed directly when the Nordsieck form is used to store the independent variables. Although a third order A-stable formula has been developed, the computational expense of this method relative to a third order backward differentiation formula that is stiffly stable, or a second derivative A-stable method, is large. The formula is only guaranteed to be A-stable for Q = f ; if Q is in error the stability region may be worse than for n y n J to J the third order backward differentiation method. Also, the solution of a linear system of up to the fourth power of hQ is necessary at each step, Even when an LU decomposition and backward substitution scheme is used, this is more expensive than for the mult i derivative formula which requires the solution of a linear system with only second order terms in hQ , thus saving on the computation to generate the third and fourth powers of Q . Also, since (f , , v-f , J must be multiplied by (b), + lib; Q ) at -n,(m+l) -n,(m) s la every corrector iteration, the corrector may be as expensive as forming the second derivative, and in the linear constant coefficient case it is. The error estimator also involves multiplication of derivatives by Q and 66 2 Q , making the error estimator more expensive. However, if one assumes v w = v (3> - & (2) cum (h) then the error can be estimated in terms of y only. However, a serious error may be committed when (H.lU) is not true. Thus, while matrix valued coefficient methods can give A-stable methods when the matrices are accurate approximations to the Jacobian matrix, they can be as expensive to use as corresponding mult i derivative methods and yet have inferior stability properties when the approximation to the Jacobian loses accuracy. k.2 Conclusions A set of formulas for the numerical integration of systems of ordinary differential equations has been presented all of whose elements are A-stable, stable at infinity, and convergent for order up to 11. Therefore, they are useful in solving stiff problems, especially those for which the eigenvalues of the Jacobian matrix of the system at some point (y , t ) have large imaginary parts. In return for the high order of these A-stable methods, derivatives higher than the first are required to solve first order systems of equations . It must be left to a potential user of the computer program DU to decide whether this tradeoff is feasible; however, several convenient methods for providing the higher derivatives were presented and should lighten the burden on any potential user. An alternative approach, previously analyzed by Lambert and Sigurdsson [31], was presented in the matrix valued coefficient formulas. However, it was noted that even the simplest such matrix valued coefficient formula provides a method which requires a large amount of work in the form of matrix evaluations and manipulations . Further, the methods are only 67 guarantee} t^ be A-stable for an exact Jacobian and are not necessarily stable at infinity. The program listed in Appendix I has been tested with subroutines DIFFUN based on the tests in Appendix II. The results were indicative of the power of the underlying method, but since no "fine tuning" on an extensive set of problems over a substantial period of time was undertaken, it is appropriate to list here possibilities that might improve the performance of Di+, either by decreasing overhead, computation time, or number of function evaluations. 1. Matrix inversion - DU solves linear systems of equations using the LU decomposition and back substitution subroutines DECOMP and SOLVE from Moler [33]. If the Jacobian matrix is known to be sparse, or tri- diagonal, or have some other property for which a more effective algorithm exists, then an individual user may substitute these. 2. Corrector Iteration - Hindmarsh [27] suggests that the user be allowed four choices concerning the corrector iteration. He can allow the Jacobian to be approximated by numerical differencing, or provide it exactly. In either case, a Newton-type iteration is used. He can also approximate the Jacobian with a diagonal matrix, in which case the iteration becomes a chord method. This saves time computing PW, but may increase the number of calls to DIFFUN. Or he may use straight functional iteration, so no computation time is required to find PW; however, the corrector will only converge for small stepsize h. 3. Norms - The norm in which the error is calculated and in which the corrector convergence test is applied can be modified. Currently, the error estimate E is computed with the Euclidean norm, but the max norm is n 68 used to decide if the corrector has converged. If the Euclidean norm is used instead, requiring more computation, a look ahead feature in which the corrector value at the next iteration is "predicted" by multiplying by a previous ratio of two consecutive correctors can be used. If it is "predicted" that the next correction will be much smaller than the con- vergence criterion, then that correction will not be done. Such a change to DIFSUB resulted in a 10$ reduction in calls to DIFFUN on a set of test equations . k. Heuristic tuning - numerous bookkeeping operations involve constants that were picked heuristically . An example is the values 1.2 and 1.3 used to divide next stepsizes in the order and stepsize selection algorithm. Also rules for deciding how often to update PW and how to control for repeated failures to meet the error criterion or for failure of the corrector to converge, were all picked heuristically on the basis of experience in using DIFSUB. Further fine tuning on Dk is also possible. However, even without this additional work of fine tuning the program DU can be used to solve very stiff systems of ordinary differential equations using large stepsizes due to the high order available in its A-stable formulas. 69 LIST OF REFERENCES [ l] Barton, D. , Winers, I. M. , and Zahar, R. V. M. , "Taylor Series Methods for Ordinary Differential Equations — An Evaluation," Mathematical Software , Academic Press, New York, 1971. [ 2] Bashforth, F. and Adams, J. C. , Theories of Capillary Action , Cambridge U.P., New York, 1883. [ 3] Bjurel, G. , Dahlquist, G. , Lindberg, B., Linde, S., and Oden, L. , "Survey of Stiff Ordinary Differential Equations," Dept. of Information Processing , Royal Institute of Technology, Stockholm, Report #NAT0.11, 1970. [ h] Blue, J. L., and Gummel, H. K. , "Rational Approximations to Matrix Exponential for Systems of Stiff Equations," J. Comp. Physics , 5, pp. 70-83, 1970. [ 5] Brown, R. L. , "Second Derivative Methods for the Solution of Stiff Ordinary Differential Equations," Dept. of Computer Science , University of Illinois at Urb ana-Champaign, File No. 877, 1973. [ 6] Brown, R. L. , "Recursive Calculation of Corrector Coefficients," ACM SIGNUM Newsletter , 8, No. k, 1973. [ 7l Brown, W. S., ALTRAN User's Manual , Bell Laboratories, Murray Hill, New Jersey, 1973. [ 8] Cavendish, J. C, Culham, W. E. , and Varga, R. S. , "A Comparison of Crank-Nicolson and Chebyshev Rational Methods for Numerically Solving Linear Parabolic Equations," J. Comp. Physics , 10 , pp. 35^368, 1972. [ 9] Coppel, W. A., Stability and Asymptotic Behavior of Differential Equations , D. C. Heath, Boston, Massachusetts, 196 5. [10] Curtiss, C. F. , and Hi rsch f elder, J. 0., "Integration of Stiff Equations," Proc. Nat. Acad. Science , U.S. , 38 , pp. 235-2^3, 1952. [ll] Dahlquist, G. , "Numerical Integration of Ordinary Differential Equations," Math. Scandinavica , h_, pp. 33-50, 1956. [12] Dahlquist, G. , "Stability and Error Bounds in the Numerical Integration of Ordinary Differential Equations," Trans. Roy. Inst . Tech . , Stockholm, No. 130, 1959- [13] Dahlquist, G. , "A Special Stability Problem for Linear Multistep Methods," BIT, 3, pp. 27-^3, 1963- TO [Ik] Ehle, B., "High Order A-Stable Methods for the Numerical Solution of Systems of Differential Equations," BIT , 8, pp. 276-278, 1968. [15] Ehle, B., "A Comparison of Numerical Methods for Solving Certain Stiff Ordinary Differential Equations," Dept of Mathematics , University of Victoria, Victoria, Report #70, 1972. [16] Enright, W. , "Studies in the Numerical Solution of Stiff Ordinary Differential Equations," Dept. of Computer Science , University of Toronto, Toronto, Report #H6, 1972. [17] Fike, C. T. , PL/I for Scientific Programmers , Prentice-Hall, Englewood Cliffs, New Jersey, 1970. [18] Gear, C. W., "The Automatic Integration of Stiff Ordinary Differential Equations," Information Processing , 68 , ed. A. J. H. Morrel, North Holland Publishing Company, Amsterdam, pp. 187-193, 1969. [19] Gear, C. W., Numerical Initial Value Problems in Ordinary Differential Equations , Prentice-Hall, Englewood Cliffs, New Jersey, 1971. [20] Gear, C. W. , "The Automatic Integration of Ordinary Differential Equations," Comm. ACM lU , pp. 176-179, 1971. [21] Gear, C. W. , "DIFSUB for Solution to Ordinary Differential Equations," Comm. ACM , lh_, pp. 185-190, 1971. [22] Gear, C. W., and Tu, K. W. , "The Effect of Variable Mesh Size on the Stability of Multistep Methods", submitted to SLAM J. Num. An. , 1973. [23] Gear, C. W. , andWatanabe, D. S . , "Stability and Convergence of Variable Order Multistep Methods," submitted to SIAM J. Num. An. , 1973. [2U] Genin, Y., "An Algebraic Approach to A-Stable Linear Multistep- Multi-derivative Integration Formulas," MBLE Research Laboratory Report R2kk, 197U. [25] Graves-Morris, P.R. ed, Pade Approximants and Their Applications , Academic Press, London, 1973. [26] Henrici, P., Discrete Variable Methods for Ordinary Differential Equations , John Wiley and Sons, New York, 1962. [27] Hindmarsh, A. C, "GEAR: Ordinary Differential Equation System Solver," Report UCID- 30001, Lawrence Livermore Laboratory, University of California, 1972. 71 [28] Hull, T. E. , "The Numerical Integration of Ordinary Differential Equations," Information Processing , 68 , ed. A. J. H. Morrell, North Holland Publishing Company, Amsterdam, pp. UO- 53, 19 69. [29] Hull, T. E., and Enright , W. H., "A Structure for Programs That Solve Ordinary Differential Equations," Dept. of Computer Science , University of Toronto, Report #GG\ 19 TU. [30] Jeltsch, R., "Multistep Methods Using Higher Derivatives and e -St ability," submitted to SIAM J. Num. An. [31] Lambert, J. D. , and Sigurdsson, S. T., "Multistep Methods with Variable Matrix Coefficients," SIAM J. Num. An. , 9, pp. 715- 733, 1972. [32] Liniger, W. , and Willoughby, R., "Efficient Integration Method for Stiff Systems of Ordinary Differential Equations," SIGNUM, 7_, No. 1, pp. U7-66, 1970. [33] Moler, C. B. , "Matrix Computations with FORTRAN and Paging," Comm. ACM , 15, pp. 268-271, 1972. [3h] Moulton, F. R. , New Methods in Exterior Ballistics , University of Chicago, Chicago, Illinois, 192 6. [35] Nordsieck. A. , "On the Numerical Integration of Ordinary Differential Equations," Math. Comp . 16 , pp. 22-1+9, 1962. [36] Norman, A. C, "TAYLOR User's Manual," University of Cambridge, Computer Laboratory, Cambridge, 1973. [37] Ralston, A., A First Course in Numerical Analysis , McGraw-Hill, New York, 1965. [38] Shampine, L. F. , and Gordon, M. K. , "Local Error and Variable Order, Variable Step Adams Codes," 1973. [39] Skeel, R. D. , "Convergence of Multi value Methods for Solving Ordinary Differential Equations," Dept. of Computing Science, University of Alberta, Report TR73-16, 1973. [UO] Varga, R. S., "Some Results in Approximation Theory with Applications to Numerical Analysis," Numerical Solution of Partial Differential Equations: II , ed. B. E. Hubbard, Academic Press, New York, pp. 623-6U9, 1971. [ Ul] Widlund, 0., "A Note on Unconditionally Stable Linear Multistep Methods," BIT, 1, pp. 65-70, 19 67. 72 APPENDIX I Dk FOR THE INTEGRATION OF STIFF SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS 73 SU«cCUTINfc D4(N,T, Y, SAVE , H ,HMIN ,HMA X, DYS, OY, D4 OOl ♦ PPStYMAX.EBSV, SIGMA, PW , IP, K.F LAG, JST ART ,NQM*X ) D4 002 IMPLICIT PEAL*8(A-H,Q-n D4 003 q************ ************************************** ****m**************** 04 004 C* * D4 005 C* THIS SUBROUTINE INTEGRATES A SET OF N ORDINARY DIFFERENTIAL FIRST * 04 006 C* ORDFR EQUATIONS OVER ONE STEP OF LENGTH H AT EACH CALL. H CAN BE * 04 007 C* SPECIFIED BY THE USER FOR EACH STEP, BUT IT MAY BE INCREASED OP * D4 008 C* DECREASED BY D4 WITHIN THE RANGE HMIN TO HMAX IN ORDER TO * D4 009 C* ACHIEVE AS LARGE A STEP AS POSSIBLE WHILE NOT COMMITTING A SINGLE * D4 010 C* STEP ERRCP WHICH IS LARGER THAN EPS IN THE L-2 NORM, WHERE EACH * D4 Oil C* COMPONENT TF THE FPROR IS DIVIDED BY THE COMPONENTS OF YMAX. * 04 012 C* * D4 013 C* THE PPOGPAM REQUIRES FOUR SUBROUTINES NAMED * D4 014 C* DIFFUVJ(T,H,Y,DYS,NJRDI * D4 015 C* DErOMP(N,N,PW, IP) * D4 016 C* SOLVE(N,N,PW,DY, IP) * D4 017 C* MATRIX(PW,T,N, Y, EPS,H,B,DYStNORD) * D4 018 C* THE FIRST, OIFFUN, EVALUATES THE DERIVATIVES OF THE DEPENOENT * D4 019 C* VARIABLES STORED IN Y VALENCE ( EC(i),E), ( EC( 2 ) , EUP ) , ( EC ( 31 , EDWN ) , ( EC ( 41 , ENQl ) , ( 5I.ENQ2), (EC(6),ENQ3) , ( EC ( 7 ) ,BND) START. GT.OJGO TO IQO ********** ****** ************* ********************************* FIRST CALL, THE ORDER IS SET TO 1 AND THE INITIAL * IVES ARE CALCULATED. * ************************************************************** MEANINGS.. THE FIRST STEP VALUE OF JSTAPT CAN INITIALIZE D4 062 D4 063 04 064 D4 065 DA 066 D4 067 04 068 04 069 DA 070 D4 071 04 072 04 on D4 074 DA 075 04 076 04 077 D4 078 D4 07 9 D4 080 D4 081 D4 082 DA 083 04 084 D4 085 D4 086 D4 087 04 088 04 089 D4 090 04 091 D4 092 D4 093 04 094 D4 095 D4 096 04 097 04 098 D4 099 04 100 04 101 D4 102 04 103 D4 104 D4 105 04 106 D4 107 04 108 04 109 04 110 04 ill D4 112 04 113 04 114 04 115 04 116 04 117 D4 118 D4 119 04 120 04 121 04 122 75 10=2 04 123 CALL DIFFUN(T,H,Y,DYS,N0R0) 04 12* IWEVAL*1 04 125 DC 10 1*1, N 04 126 10 Y(2 f II«DVS( Iti) 04 127 HNEW*H 04 128 £***** ******** ** **** ** ******** ***************** ************************* 04 12 9 C* BEGIN BY SAVING INFORMATION FOR POSSIBLE RESTARTS AND CHANGING * D4 130 C* H BY THE FACTOR H/HNEW IF THE CALLER HAS CHANGED H. ALL VARIABLES * D4 131 C* DEPENDENT ON H MUST ALSO BE CHANGfcD. * D4 132 C* E IS A COMPARISON FOR ERRORS OF THE CURRENT ORDER NO. EUP IS * 04 133 C* TO TEST FOP INCREASING THE ORDER, EDWN FOR DECREASING THE ORDER. * 04 134 C* HKEW IS THE STEP SIZE THAT WAS USED ON THE LAST CALL. * D4 135 Q ******************* ******************************* ********************* Q4 i3 6 100 IF(H.NE.HNEW»f ALL SC AL E< Y , Y , H/HNEW, K , N 1 04 137 KFLAG*1 04 138 DO 110 1*1, N D4 139 DO 110 J«l,K 04 140 110 SAVE( J,I )«Y( J,I) 04 141 HOLD=H D4 142 TOLD=T 04 143 IF( JSTAP.T.GT.OJGO TO 200 D4 144 C ******************************************************* **************** 04 145 C* CALL COEF. THIS WILL ... * 04 146 C* SET THE COEFFICIENTS THAT DETERMINE THE ORDER AND THE METHOD * D4 147 C* TYPE. CHECK FOR EXCESSIVE ORDER. THE LAST TWO STATEMENTS OF * D4 148 C* THIS SECTION SET IWEVAL .GT.O IF PW IS TO BE RE-EVALUATED * D4 149 C***** ******** ********** ************************************************ 04 150 150 CALL COEF( A,B,NORD,NQ,K,EC,IWEVAL»EPS) D4 151 C***** ************************************************* ***************** 04 152 C* THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY * D4 153 C* MULTIPLYING THE SAVED INFORMATION BY THE PASCAL TRIANGLE * D4 154 C* MATRIX. * D4 155 C** ************************************************************* ******** 04 156 200 T*T+H 04 157 DO 210 J«2,K 04 158 DO 210 Jl«J,K D4 159 J2-K-J1+J-1 04 160 DO 210 I-l,N D4 161 210 Y( J2,I»*Y(J2,I l*Y( J2*i,I ) 04 162 C* PERFORM UP TO 3 CORRECTOR ITERATIONS. CALCULATE SIGMAUJt THE * D4 163 C* LINEAR COMBINATION OF PREVIOUS VALUES VIII AT T-J*H, J*l,2,..., ON * 04 164 C* THE FIRST PASS. IF THE CORRECTOR DOES NOT CONVERGE, GO TO 300, * 04 165 C* REDUCE H, AND SET IWEVAL-1 TO CAUSE PW TO BE REEVALUATED. IF * 04 166 C* SUCCESSFUL, GO TO 350. • D4 167 DO 300 L«l,3 04 168 230 NT«N D4 169 IF(L.GT.UGO TO 250 04 170 DO 240 I»1,N D4 171 SIGMA( I )»Y<1, I) D4 172 DO 240 J«1,N0RD 04 173 240 SIGMA(I)*SIGMA(II-B( J»*Y( J + l , I J D4 174 250 CALL DIFFUN(T,H,Y,OYS,NORDI D4 175 !F( IWEVAL. LE.OIGO TO 255 D4 176 CALL MATRIX(PW,T,N,Y,EPS,H,B,DYS,NORD) D4 177 CALL DECOMP(N,N,PW,IP» 04 178 IF( IP(NI.EO.O)GO TO 640 D4 179 IWEVAL— 1 D4 180 255 CONTINUF 04 181 DO 260 I«i,N 04 182 DY(II«SIGMA(II-Y(1,I» 04 183 76 DO 260 J-l.NOPD 260 DY< It«OYU >*DYS(I, JI»BC J) CALL SOLVE(N,N,PW,DY,IPI on 270 I«l,N Y(l,I)-Y(l,II*DY*XK-Y< J*l» I I 280 Y< J*1,I )«XK GO TO 350 300 CONTINUE H»OMAX1(HMIN,.25DO*H) IF|H.LT.1.1D0*HMINIG0 TO 630 301 T-TOLD ID-K IWEVAL-l P-H/HOLD CALL SCALF(Y,SAV6,R,K,NJ GO TO 200 350 D»0. IF(NQ.EQ.1)G0 TO 361 355 00 360 I»i,N DD«0.00 DO 358 J«1,N0RD 358 DD-OD*DYS( I,J)*A| K,J J 360 D«n»(DD/VNAXm>**2 GO TO 365 361 DO 362 I»1,N 362 0-0*((Y(2,Il-SAVE«2,I»*H/HOLDl/VMAXCni**2 365 IWEVAL'O * IF *00 KFLAG*1 ID=ID-l IF( ID.NE.11G0 TO 390 OD 380 1*1, N 380 ERSV( I)-DYS( I,NORD> 390 IFUD.GE.DGO TO 575 395 PR2-(D/E)**EN02*1.2 PR3«1.E*20 IF(KFLAG.LE.-1.0R.K.GT.NQMAXIG0 TO 420 0*0. AJ2,l)*l.D0 DO 410 1*1, N DD*A(K,NORD»*(UYS( I ,NORD»-ERSV( 1 1 ) IFCNORD.EQ.DGC TO 410 00 405 l*2,N0RD 405 DD*DD*A< K , L- 1 ) *DYS ( I ,L ) *0FL0AT ( I ) 410 D*0*(DD/YMAX< I l>*»2 PR3»«D/EUP)**ENQ3*1.3 420 PR1*1.E*20 0*0. IFU.E0.2l GO TO 460 00 430 1*1, N 430 D»DMY(K,I>/YMAXKFLAG-2 IF(H.LT.1.1D0*HMIN)G0 TO 620 T»TOLD IF(KFLA6.LE.-51G0 TO 600 GO TO 395 460 IF(PR2.LE.PR3IGO TO 510 If (PR3.LT. PRIIGO TO 520 470 R»l./AMAXHPR1 ,1.E-41 NEWQ>NQ-1 480 IF(KFLAG.EQ.l.ANO.R.LT.l.l)GO TO 575 IFINEWO.LE.NQIGO TO 500 IF1N0.EQ.11G0 TO 495 DO 490 1*1, N Y(NEW0>1,I)«0.D0 DO 485 L-lfNORO 485 Y(NEWQU,I»*Y(NEWQ*1 ,II*0YSU,L1*A(K,L) 490 Y = SAVE(1,I J Y<2,I)=SAVE<2,I)*R 610 CONTINUE KFLAG«1 GO TO 150 620 KFLAG=-1 GO TO 650 630 KFLAG=-2 GO TO 650 640 KFLAG=~3 650 CALL SCALE'Y,SAVE,1.D0,K,N> h=HOLD T=TOLD GO TO 585 END D4 306 04 307 04 308 04 309 04 310 04 311 D4 312 04 313 D4 314 D4 315 04 316 D4 317 D4 318 04 319 D4 320 04 321 D4 322 04 323 04 324 D4 325 04 326 D4 32 7 04 328 D4 329 D4 330 D4 331 04 332 D4 333 04 334 04 335 04 336 D4 337 04 338 04 339 04 340 D4 341 04 342 04 343 04 344 04 345 D4 346 04 347 SU8R0U c*********** C* THIS SUBR C* FORM OF T C* MI, J) B(J) NORO NQ K IWEVAL EPS EC C* C* C* C* C* C* C* C* C* C* C* C* THE RATIO C* SET OF AS C* GIVEN IN TINE COE ******** OUTINE P HE MULT I THE COE THE COE THE CUR THE CUR SET TO SET TO THE ERR PARAMET AT OP IN 04 N£L VALU SIGNMENT THE SAME F< A, B, NORO, NQ,K, EC, IWEVAL, EPS) **************************************************** LACES THE CORRECTOR COEFFICIENTS FOR THE NOROSIECK * DERIVATIVE METHOD OF D4 INTO A AND B. * * FFICIENT OF DVS(L,J) TO UPDATE YII.Li. * FFICIENT OF 0YS5 c**** C* 10 r* c* c* 156 157 C* 52 C* C* c* c* c* c**** S( 1 ) = NOPD = 8(2»» P ( 3 ) ■ A< 5,1 A( 6,1 M 7,1 4(5,2 M6,2 A( 7,2 A( 5,3 A(6,3 4(7, 3 r,c TO » *. *** 3020/1 22854 -23 I ****** <3( 1 ) = NOP.D- B( 21* P( 3» = A(5,l A(6, 1 A(7,l A( 8,1 A( 5,2 A(6,? A( 7,2 A(8,2 A(5,3 A(6,3 A( 7,3 A(8,3 GO TO NQRD* **< *** 491180 135270 ■♦925 28 fco201l)2651 B3 92U0 61591 1 956372 .25523 .122/5 .01550 - .6212 -.2694 -.0322 .99820 .32934 .03550 90 03507 96834 95209 44910 4 7048 5 74 85 61077 92557 35928 131 73 04277 2711 7D0 901600 580800 1796400 75962D0 02 99 00 8443D0 74166D0 1437D0 652700 159900 2 49 7 7,-493200/ 8 7485 3, 2 16000 Z87 4853, 356895/999832 t 3/9 998 32 ,512 2 5/9998 32,390 5/999832,-40 6285/499916, 42 85/499916,-3483 3 7/3499412,-36 75/499916,29422 5/2499 58, 280 4 5/2*9958, 172065/1749706,1717/24995 8 82429848214500 5o37518 4689862 .356954 .228581 .051233 .003905 -.81270 -.46864 .09954 -.00735 1.17709 .51 5466 .098339 .006869 90 531684 182 560 96 84 34 401675 60 7246 o515J2 653469 873298 158012 123500 775242 598388 378158 154017 800 857D0 687D0 48200 01700 332400 782700 714D0 8318100 74812D0 24D0 529D0 3878D0 8 7 5 00 /58067611 ,43268400/58067611,29592000/58067611, 12960000/5806/6 9 75/9290 817 76,3178 5985/92 90817 76,428418347/116135222 000, 48490/92 90817 76,32 38 7 44 74/92 90817 76,1038 652097/11613522200, 4 3 995/46 4 540 888,162 79665/2 32270444,12 7918665/2322 70444, 760 0419/2 322 7 0444,20595 75/2322 70444,6 4038605/116135222, 12401655/116135222. ****** B( l) = 9(2) = B( 3J = B(4)« A(6, 1 A(7,l A(8,l A(9,l A(6,2 A(7,2 A(8,2 A(9,2 A(6,3 A(7,3 .903966584745495D0 -.7451382837155099DO 09612837352647800 2231881039500 -.21010653318422800 .14559641410941400 .034212257544057100 - .0036889613643656500 .53014546482721100 .3485963C483163400 .08 9434 71 92 10 3401 DO .0061221629214266900 -.93115448214325700 -.5 50 72 8981256026 00 .5 04 428 04 429 04 430 04 431 04 432 04 433 04 434 04 435 04 436 04 437 04 438 04 439 04 440 04 441 04 442 04 443 04 444 04 445 04 446 04 447 04 448 04 449 04 450 04 451 04 452 D4 453 04 454 04 455 04 456 04 457 D4 458 D4 459 04 460 04 461 04 462 04 46 3 04 464 D4 46 5 04 466 04 467 04 468 04 46 9 04 470 04 471 04 472 04 473 04 474 D4 475 04 476 04 477 04 478 04 479 04 480 04 481 04 482 D4 483 04 484 04 485 04 486 D4 487 04 488 81 A( 8, 3 1 — . 11 882 8803 71986400 A (9, 3)*-. 00886 7 14 195 974 150 J A (6, 4 1*1.2 1 8972 76 5902 1600 A( 7, 4)*. 55141415237490900 A (8, 4) =.106 78633 73 9556 300 A( 9 ,4) =.00 75 3 3 0548 72 8797 600 GO TO 190 158 NCRD=4 *( l) = . 88879254103287900 R< 2) =-.71 2 394495 589035900 B< 3)=. 46 58 106 19845461200 B(4)»- . 190126783610400 A (6, I » =-.2 8 17900 730 1499400 A( 7,1) =-.2 3558395893 83800 A (8,1 I — .0 751100484661 83 9D0 A (9,1) =-.0106 75267662 145900 A( 10, I )=-. 0005665640 8668 386200 A( 6,2)=.684830C488991200 A (7, 2 1 =.542 778 762 52 71 5D0 A( 8,2)=. 16768 7305413 30 800 A( 9, 2 I = .023355665492022200 M 10,2)=.00122257815818899D0 A (6, 3) =-1.133 8080075796600 A (7, 3) — .81049190 840915900 A (8, 3 I =-.2 368865849852 5900 A( 9,3)=- .0319208468362 96400 A( 10,3) — .00 16354746 1881 029D0 A( 6,4) = 1.3851 5 740055 29600 A( 7, 4)«. 74747969050837600 A (8, 4) = .195894 75 124 1 68500 A (9, 4) -.024933678980902600 At 10, 4)>. 001234434084782200 GO TO 190 159 N0RD=4 B(l)».8754775810l935600 8(21 — .684913223 756951900 B< 31-. 43131650690545400 «M 41 — .166347964095 7473100 A (6, I) =-.35636313612 07600 A( 7,11— .34119753667435700 A (8,1) — .132 882 15203663600 A (9,1) =-.02 5950780 I 8 8853 200 A( 10, II— .002522714361042200 A( 11,1)— .0000972666 74731368700 A (6,2) -.83 87443 3902 66600 A (7, 2 1 -.760758796752 07200 A (8, 2) -.28692543 162002400 A(9,2 )-. 054883399085987800 A( 10,21-. 0052599539655092700 A( 11, 2)». 00020075253142186100 A (6, 3 1 — 1.331271167503300 A ( 7,3) — 1.08409743 74 3743 700 A (8, 31 — .38655263455005800 A( 9, 3 I — .071494022705914600 A( 10, 3)— .0067031329 83 7448400 A( 11,31— .00025198180543352200 A (6,4) =1.51 83 35 3 752 845600 A( 7,4)-.9360920634U8)8D0 A (8, 4) -.2990683 5614792 9D0 A (9, 4) -.052213 7951 05650500 A( 10,41-. 004727868858163400 04 489 04 490 04 491 04 492 04 493 D4 494 04 495 D4 496 04 497 D4 498 04 499 D4 500 04 501 04 502 04 503 D4 504 04 505 04 506 04 507 04 508 04 509 04 510 04 511 04 512 D4 513 04 514 D4 515 04 516 04 517 04 518 04 519 04 52 D4 521 04 522 04 523 04 524 04 525 04 526 04 527 D4 528 04 52 9 04 530 04 531 04 532 04 533 D4 534 04 535 04 536 04 537 D4 538 D4 539 04 540 04 541 04 542 04 543 04 544 04 545 D4 546 04 547 04 546 04 549 82 A( 11, 4)=. 0001 73705869253642DO 04 550 r,C TO 190 D4 551 loO N0P0=4 04 552 R| 1 )=.8636388197J0276D0 04 553 IM 2)=-.6fcl38b334459858D0 04 554 »M ">)= .<*u33025338O76807O0 0* 555 B<4)=-. 1483898941735602900 D4 556 A(6, 1) =-.43283540780964800 04 557 A( 7,1 )=-.460011012203615D0 04 558 A(8,l )=-. 207181319492200 04 559 MS, 1 )=-. 049911494634830600 04 560 A< 10,1 )=-. 006733124734215300 04 561 A< 11, 1)=-. 0004802856918797900 04 562 A( 12, 1)=-. 0000141207653807700 04 563 A(6,2)=.99071586935629D0 D4 564 M 7,2 )=. 99687399643370800 04 565 A(8, 2)=. 4345784034885800 D4 566 A(9,2I=. 1024999581561700 04 567 A( 10, 2)=. 013627202558920700 04 568 A( 11,2 )=. 00096191713259789600 04 565 A( 12, 2) = . 0000280618618618618600 04 57C A(6,3)=-l. 5122269372386700 D4 571 A(7,3)=-l. 36524488400800 04 572 A(8, 3)=-. 56236620384223300 04 573 A(9, 3)=-. 128192083 71779600 04 574 M 10, 3)=-. 016666195844873900 04 575 A( 11, 3)=-. 0011583168670151300 04 576 A( 12,3)=-.0000334138624662642DO 04 577 A(6,4)=l. 6343352141625800 04 578 A(7,4)=l. 1163187483711800 04 579 A(8,4)*.411718402099D0 D4 580 A(9,4)=.0885595051740368D0 04 581 A( 10, 4>=. 01111458804538900 D4 582 A( 11, 4)=. 00075470265794556500 04 583 A( 12,4)=.000021419613577347D0 04 584 190 K=NQ+1 04 585 196 EC(5)=.5/DFLCAT=»)**2 04 593 EC<7)=(EPS*EC<6)) 04 594 IWEVAL=1 04 595 200 PETUPN D4 596 END 04 597 SUBROUTINE SCALEiY, YS,P,K,N) 04 598 IMPLICIT BEAL*8 (A-H,Q-Z) D4 599 DIMENSICN Y<12,1) ,YSU2,l) > D4 600 P1=1.D0 04 601 DO 10 J=1,K 04 602 DO 9 1=1 ,N 04 603 9 Y( J, I )=YS( J, I )*Rl 04 604 10 Rl=Rl*R 04 605 RETURN 04 606 FND 04 607 83 SUBROUTINE MATPIX(PW,T,N,Y,EPS,H,B,DYS,N0RD1 ^m* ************************************************************* ******** C* THIS SUBFCUTINE 7*KES N CALLS TO DIFFUN TC FIND THE PAFTIAL C* OEPIvrTIVE WITH RESPECT TO Yd), 1 = 1, N, OF THE JTH DERIVATIVE IN T C* OF EACH CF THE VARIABLES , FOR J*l,NOPD. THE CHANGE IN Yd) FOR C* FCPWAPO DIFFERENCING IS R = MAX< EPS*Y < I ) , EPS*EPS 1 . r* THESE A°F USED Tf CALCULATE THE MATRIX PW. C* C* NCFD C* PW I - SUM H**J * 8E T A(J) * DF|(J-1))/DY C* . J=l C* WHERE ((J)) IS THE JTH DERIVATIVE IN T. C* T THE INDEPENDENT VARIABLE. C* N THE NUMBER IF EQUATIONS; THE CPDER OF PW. C* Y DEPENDENT VARIABLE ARRAY AS IN D4. C* EPS ERRCR CPITERI IN C* H STEPSIiE C* eETA(I) COEFFICIENT USED IN D4 = R ( I ) / FACTOR I AL ( I ) . C* DYS(I,J) THE JTH DERIVATIVE ELEMENT OF THE TAYLOR SERIES FOP Y(I) C* N1RD THE CURRENT HIGHEST DERIVATIVE BEING EVALUATED. Q*m ************************************************ ********************* IMPLICIT PEAL*8 (A-H.Q-Z) COMMON NFNS.NS DIMENSION PW(N,N),Y( 12, 1 J , B< 4) ,DYS( N, 1 ) NFNS=NFNS-N DO 5 1=1, N 0^ 5 J=1,N PW( I, J)=0. IF(I.EQ.J)PW( I ,J)=l. 5 CONTINUE DO 20 1=1, N R = -EPS*DMAX1(EPS,DABS(Y( 1,1) ) ) F=Y(1,I) Y(l,I)=Y(l,I)~P CALL OIFFUN(T,H,Y,OYS< 1,5) ,NORD) DO 10 K=1,N DO 10 J=1,NCPD 10 PW(K,I )=PW(K,I ) ♦ Bl J)*(DYS(K, J+4)-DYS(K, J) )/R 20 Y(1,I)=F RETURN END 04 606 * 04 009 * 04 610 * 04 611 * 04 612 * 04 613 * 04 614 * 04 615 * D4 616 * 04 617 * 04 616 * 04 619 * 04 620 * D4 621 * 04 622 * D4 62 3 * 04 624 * 04 625 * 04 626 * D4 627 * D4 628 D4 629 04 630 04 631 D4 632 04 633 D4 634 04 635 D4 636 D4 637 04 638 D4 639 04 640 04 641 04 642 04 643 D4 644 D4 645 04 646 04 64 7 04 648 8U SUBROUTINE DEC^Pdj, NUIM.A, IP) D4 CtC TO 5 04 671 KP1 = K+l 04 672 M = K 04 673 DC 1 I = KPl,N D4 67<» IF(ABS(A(1 ,K)) .GT.ABS(A(M,K) )) M= I D4o75 1 CONTINUE D4 o76 IP(K) = M 04 677 IF(M.NE.K) IP(N) = -IP(N) 04 678 T = A(M,K) D4 679 A(M,K) = t (K,K) 04 680 A(K,K) = T 04 681 IF(T.EO.O) GO TO 5 04 682 DO 2 I=KP1,N 04 683 2 A(I ,KI = -A( I,K) / T 04 684 DC 4 J=KP1,N 04 685 T = A( M, J) 04 686 A(M,J)=MK,J) 04 68 7 A(K,JJ = T D4 688 IF(T.EO.O.) GO TO 4 04 689 DO 3 I=KP1,N 04 690 3 A(I,J) = A(I,J) *■ A(I,K)*T D4 691 4 CONTINUE D4 692 5 IF(A(K,K».EQ-0. J IP(N) =0 04 693 6 CONTINUE 04 694 RETURN 04 695 END D4 696 85 SUBROUTINE SOLVE ( N, NOI M, A, B, I P ) 04 697 IMPLICIT REAL*8 (B-H,Q-il D4 698 C D4 699 C SOLUTION OF LINEAR SYSTEM, A*X = B. 04 700 C 04 701 C INPUT... 04 702 C N = ORDER OF MATRIX.. D4 703 C NDIM = DECLARED DIMENSION OF ARRAY A. D4 704 C A = TRIANGULAKIZED MATRIX OBTAINED FROM 'DECOMP'. 0* 705 C B = RIGHT HAND SIDE VECTOR. 04 706 C IP ■ PIVOT VECTOR OBTAINED FROM «DECOMP«. 04 707 C OUTPUT... 04 708 C B = SOLUTION VECTOR, X. D4 709 DIMENSION A(NDIM,NOIM) , B(NDIM), IP(NDIM) 04 710 IF(N.EQ.l) GO TO 9 04 711 NM1 * N-l D4 712 DC 7 K=1,NM1 D4 713 KP1 = K ♦ 1 D4 714 M « IP(K) 04 715 T=B(M) D4 716 D4 717 04 718 04 719 04 720 D4 721 04 722 K = KMl +1 04 723 B(K) = B(K)/A(K,K) D4 724 T = -B(K) D4 725 DO 8 1=1, KMl 04 726 8 B(I) * B(I) ♦ A(I,K)*T D4 727 9 B(l) « B