033390-2-T THE UNIVERSITY OF MICHIGAN COLLEGE OF ENGINEERING Departments of Aerospace Engineering Meteorology and Oceanography High Altitude Engineering Laboratory ehnical Report.': AMPLCT A NUMERICAL INTEGRATION ROUTINE FOR SYSTEMS OF STIFF DIFFERENTIAL EQUATIONS By. L Stuart W. Bowen ORA Project 033390 under contract with: NATIONAL AERONAUTICS AND SPACE ADMINISTRATION Grant NGR 23 -005 -360 WASHINGTON, D. C.

......

CONTENTS Page Introduction. Derivation of Equations Accuracy and Stability 6, Calling Requirements of AMPLCT 13 FOR AN IV Listing of AMPLCT and MINV 17 vi Comments for Subroutine AMPLCT 27 VI Modifications to AMPLCT for Eigenvalue Display 29 VIII. Example Solution 30 IX Comparison to Other Methods 35 X ~ References

I. INTRODUCTION The numerical integration of the ordinary (nonlinear) differential equations describing the steady flow of a reacting' gas by the usual Runge-Kutta, Adams-Moulton or Hamming's methods typically are found to require a vanishingly small stepsize because of an extreme stability property of the differential equations known as "stiffnesst. Physically this situation arises from the very large spread in the time scales associated with the various hydrodynamic and chemical rates in the flow. The step size required by these usual explicit methods is always determined by the fastest rates even though the overall solution may not be determined by these fastest solution components. One would like to use a step size more in keeping with the changes in the overall solution. Rather than discussing this phenomena further here, we refer the interested reader to Ref. 1-7 and the references contained therein. Of particular interest are Chapters 3 and 6 in Ref. 7. In the follow ing sections the equations required to integrate large sets of stiff differential equations which may have physically meaningful parasitic eigenvalues are developed following the procedure given in Ref. 2, 4-6. Since an algorithm which implements this very powerful

requirements and helpful explanatory comments concerning the computer program. Some remarks are made regarding the program modifications needed if the matrix eigenvalues are to be computed and displayed, and an example solution is provided. Finally, comparisons between this method and other methods are given for several systems of equations. II. DERIVATION OF EQUATIONS A set of coupled differential equations can be written in vector notation as ds,= d= F) (1) The equation is assumed to have been written in the automorphic form so that the vector F does not explicitly contain the independent variable s. We then locally linearize the right hand side by performing an expansion of the ith component of the vector F (denoted by F.) in a multidimensional Taylor series about the point sn = nh. This gives, after substitution in (1) for the ith component of 'WI -(denoted by wi) w iI Fin + (w1 - win) (aF /awp1+.. + (Win - wmn) (a Fi/awm)n (2) +O0(!W-Wn12) i I to m Thenumer f ompnens f ~is akn t em ec hr r

Note that when W -=W(s + nh) W w -w12 1W ~ n n+n = h2lWn' ( Define the Jacobian matrix [A] = (aij) = (aFi/awj) (3) In many problems, explicit expressions for the matrix elements, ij cannot be conveniently given. In such cases the elements can be umerically computed by the approximate formula aij = aF i/awj [Fi(1.005 wj) - Fi(. 995 01 w (3a) which is accurate to third order. In terms of [A] evaluated at s n and denoted by [An], Eq. (1) can be written W'= Fn+ [An]( w) + 0(h ) =[A] W+ FK-An] W+O0(h2) or W 4A]W+f +0O(h2) (4) where~~~~~~n

Note although Eq. (4) represents the local linearization of Eq. (1) (referenced to sn = nh), it is the derivative W' that is represented as linear in the step-size h. We use the modified Euler implicit difference scheme which is consistent with this error and is unconditionally stable to advance the solution w n+lW + (h/2)(w w) (6) The local linearization has converted the set of differential equations into a set of m algebraic equations that must be solved by the methods Of linear algebra. Inserting Eq. (4) into (6) gives n+1 wn+ (h/2)([An] Wn+ + fn+ [An] W + f ) which can be rearranged to form ([I] -(h/2)[An](i-)=h[An Wn + hf =h([ An] Wn + FnK-[An] Wn) =h Fn where [I] is the unit matrix. The increment in the dependent argument vector is given by &W= W1 w =([I] - (h/2)[A]1)1 hFn (7)

where ([I] - (h/2)[An] is the inverse of the matrix I] (h/2)[A] A measure of the "stiffness" of Eq. (1) is given by the ratio of the greatest to least eigenvalues of the matrix [A]. Standard explicit methods such as the fourth order Runge Kutta typically require a step size less than 2. 78/ max Ifor stability, where X is the largest eigenmax max' value of [A]. This small step size is governed not by the behavior of the solution as a whole but by the rapidly varying transients, i.e. largest eigenvalue, which may have ceased to be numerically important to the overall solution at large times. The method by which the step size is controlled is quite important to the success and usefulness of a method. This will be illustrated in Section IX in which the same formula as Eq. (7) is used to advance the solution under two other step size control techniques. The technique used here is described as follows:0 1.Compute the increment in the independent argument iAw.1 wi.n)Win for each i Ito m. 2. If any 1w. IK<E, ignore. it in the following tests. The value of E used is 1 x 10 3 3. Find the largest I& ~W/w I = I &w/w J: max ~ i ax th4acuae incre ments> Aw, adivd rheaculaete sthep stzep. aldicr

5. If IA&w/W max < 6/b, double the current step size, discard the calculated increments Awi and recalculate the step. 6. If 6/b < I Aw/w ma <, the calculations are considered satismax factory, and the calculated step is accepted. Values for b and 6 that have proved satisfactory are b=5 6 =.05 This step size control works well for implicit techniques. It decreases the step size in regions where the largest eigenvalues are important in numerically determining the solution. Since implicit techniques are stable for all negative eigenvalues, the value of the maximum eigenvalue itself should not be used to control step size. All implicit techniques require the inversion of a matrix, usually at each step, to advance the solution. This matrix inversion typically requires of order N3/3 operations for a matrix of order N and so becomes the operation requiring the most time for very large N, (say of order 50-100). For small N the evaluation of the derivative vector F is usually the most time consuming operation. IIIL ACCURACY AND STABILITY.

properties of the trapezoidal rule used in Section II. For a fuller discussion Ref. (5,7,8 and 10) have proved helpful. The accuracy of the numerical solution to an ordinary differential equation depends on two things; the truncation error at each step due to the finite step size and the retention of only a finite number of terms in the equivalent power series solutionaand the stability of the finite difference method used. Stability of the finite difference method refers to the manner in which unwanted solutions that are inevitably introduced, are propagated from one step to another. We shall not be concerned Fo8 here with what calls inherent instability, wherein a small amount of a strongly increasing complementary solution to the differential equations, whose coefficients should be zero for the given initial conditions numerically, is later introduced and swamps the wanted solution, or strong instability in which spurious finite difference complementary solutions are introduced through the use of higher order difference equati~ons than the differential equations to which they are applied. Fox refers to the stability of interest to us as partial instability since in principle it can always be overcome by selecting a sufficiently small step size. To test the stability of a numerical method we consider a linear set

'- [A] W(x) (1) where [A] is a constant matrix which may be identified with the local Jacobian matrix of Section II for a nonlinear set of equations. The initial conditions are presumed to be specified by W (O). 7 It can be shown that if the matrix [A] is diagonalizable, that is, there exists a matrix [S] such that [S] A] [ S] = [I] X, then we can write yIs] W and transform Eq. (1) to the form '?: [I]3A y. (2) If [ S] [ A] [ S] 1 = [I] X, then multiplying from the left on [S] - 1 and transforming the right hand side yields [A] = [I] I so that the elements of the vector X are the eigenvalues of the matrix [A]. Expanding Eq. (2) we see that it is therefore sufficient to discuss a set of linear differential equations of the form y = Xi yi (3) for each i. Note it is only the eigenvalues Xi of the matrix [A] that are important in determining the behavior of the solution. The analytic solution to Eq. (3) is y = e y(O) (4) 8

for each component. We apply the modified Euler finite difference method h I Yn+1 =n + 2 (Yn+1 + Yn) to Eq. (3) and find hX hX Yn+ 1 Yn 2 Yn+ 1 2 - n or hX n+2\ Yn+ i 1 hX n (6) The quantity in the brackets is known as the characteristic root and for the trapezoidal rule hX 2 The characteristic root when the simplest first order Euler method ~n+1 =+ hyl is used can easily be shown to be Yn =( + hX) yn 8 and

To understand the importance of the p let Y = the true value of the omputed solution y which also satisfies Eq. (3), i. e. Y' = XY If e = Y - y is the error, then the error at the nth step is e= Y - Yn Now e' Y'- yt _y (Y- Y) or e' = Xe for small values of e. Thus the error e satisfies the same differential equation as does y. Applying the modified Euler method to the error equation yields hXk 2 e =+ hX en and clearly if the absolute error is not to grow with n we must have 1X I1112 Inheentinstbilty i avide byrqiigal< 0 fahge 2~~~~~~~

characteristic root p. If some of these p are greater than unity then the corresponding difference solution will grow faster than the true solution. This is the case of strong instability. For the simple Euler method we have Ip < 1 only for -2 < hXi < 0 while the modified Euler or trapezoidal rule allows for all i -oo< hXi < 0 in order that 1,<1 Thus the trapezoidal rule is at least stable for any negative eigenvalue no matter how large. The characteristic roots of Eq. (8) and (6) are special cases of what are known as rational Pade approximations to the true solution ehX giver by Eq. (4). In general any implicit method generates a characteristic root pi(x) (where x = hX) that can be identif ied with a particular Pade approximation to exp(x) which is written in the form n m PnmZ ax/Z bxk i=O k= Equation (8) is Pt owhile Eq. (6) is P1 Even though the trapezoidal rule is stable for an arbitrarily large

an accurate solution for some X's whose absolute value is closer to zero and which may be more important in numerically determining the local solution. Note as hX - -co, P — 1 for the implicit trapezoidal rule while exp(hX) - O. Even though the "stiff" components are not then accurately calculated, their numerical contribution to the overall solution is small. The virtue of the implicit method is that these stiff components no longer important in determining the solution, do not influence the stability and accuracy of the overall solution, and the step size can be adjusted to the remaining components to give accurate results. In other words the step size h can be freely adjusted to give accurate resultswith a practical value tailored to fit just those components most important in determining the local solution. This cannot be done using P, O i. e. the simple Euler method. It 'is clear that once any component of a solution has gone unstable, subsequent values of the whole solution are in error. This may happen even if the truncation error per step may be unimportant. If X. > 0, Eq. (7) indicates IMj.I> I and hence the absolute error per step will increase. In this case however the true solution is also increasi ng in value and it is the relative error that is important in determining the step size. The solution will be valid as long as no spurious component grows faster than the true solution.

For positive eigenvalues, Xi > 0 it is clear that P 1 is a good approximation to exp(hX) only for values of hX considerably less than 1. For hX -. 4, the relative error of P 1 is 5. 48 x 10. The allowable step size is therefore limited by the largest positive eigenvalue for the case of the parasitic eigenvalues. Local eigenvalues of real physical problems 9 can be locally parasitic9, i. e., have opposite signs. This is the strong 8 instability case of Fox8. In summary we see that for the modified Euler method h must be chosen to have I hX < 1 for the most important eigenvalues that determine the local solution but that the presence of very large negative X. will not 1 affect the stability of the solution. IV. CALLING REQUIREMENTS OF AMPLCT Name: AMPLCT Purpose: To obtain a solution of a system of first order ordinary differential equations dY./dx= F (X YEQilIto NEQ which show extreme numerical "stiffness", and have parasitic eigenvalues. Calling Sequence, CALL AMPLCT (VI NEQ, DERW,) QMAT, LI M) FORTRAN:

Parameters: V is a real vector of dimension > 8*NEQ + 3 containing the following information: V(1) is an accuracy parameter controlling step size changes. V(1) is the maximum relative change in any of the independent variables. A value of V(1) -. 05 has proved satisfactory in past use. V(2) = x, initial and current independent variable. V(3) = h, the initial and current step size. V(4), V(5),... V(4 + NEQ - 1) contain the NEQ values of the Yvector, both the initial and current values. V(4 + NEQ),... V(4 + 2NEQ - 1) contain the NEQ values of the dY(I)/dx vector which are supplied by the derivative subroutine i.e., DERWI). NEQ is the number of differential equations, written in the form dY(I)/dx = F(I). These must be written in the automorphic form ifx appears explicitly in the dYV/dxr vecto nr. This. is- donea by havring,

the form dx/dx = 1 and using the dependent variable x in evaluating dY/dx. As written NEQ < 50. DERIV is the name of a subroutine supplied by the user that evaluates the derivatives dY/dx and stores them in V(4 + NEQ) through V(4 + 2*NEQ - 1). The name of this subroutine must be declared EXTERNAL. This subroutine has no arguments. The variable V must be declared in a block COMMON between the main program and DER IV. In the dY/dx evaluation in DER IV, the dependent value of x (if it appears) must be used, not the independent value of x. QMAT is a working matrix, which must be dimensioned exactly NEQ by NEQ in the main program. L, M are integer work vectors which must be dimensioned exactly NEQ in the main program. QNA MA T.L N% area used by the matrix inupverion routine.

Usage: Prior to calling, the user must set V() 05 recommended) V(2) = initial value of x, V(3) initial value of step size and V(4) through V(4 + NEQ - 1) to the NEQ values of the initialY vector. Each call to AMPLT advances the solution one step. Upon return the new value of the independent variable x is stored in V(2) and the corresponding values of the solution vector Y are stored in V through V(3 + NEQ). Only 8*NEQ + 3 elements of V are needed. The auxiliary vectors Q and C used internally in AMPLCT need dimension> NEQ**2. As written, the vectors Q and C limit NEQ < 50. Successful solutions have been obtained to systems of equations having the ratio. of the absolute value of the greatest to least eigenvalues on the order 8 of 10.Situations in which physically meaningful positive and negative eigenvalues occur (parasitic case) pose no special problem although the step size must be reduced inversely with the magnitude of the

The subroutine can be converted to double precision by removing the C from the IMPLICIT REAL*8 statement in AMPLCT, line 2, and in MINV, line 43, changing the function references ABS, AMAXI in lines 50 and 51 of AMPLCT and line 69 of MINV to DABS and DAMAX1 respectively. In addition in AMPLCT, lines 12, 18, 25 and 50 the single precision power E should be changed to the double precision D power. The subroutine MINV called from AMPLCT is the standard IBM routine, supplied from SSP. V. FOR hAN WV LISTING OF AMPLCT AND MINY The leading carat > indicates the beginning of a line, while the leading star * in MINV signifies the continuation of a wrapped around line. The leading number at the left is the line number and is not part of the FORTRAN statement.

I ~~SUBROUTINE AMPLCTCVPNEQ*DERIV*QMAT*L*t4).2 C, IMPLICIT REAL*8CA-H*0-Z) 3 DIMENSION Q(2500)PCC2SOO)*V(1)D DIP 4 1 QMAT(NEQpNEQ)pL(1)*M(1) lip 5 CALL DERIV 6 DO I JEQlI*NEQ > ~7 V(2*NEQ+3+JEGQ=V 3+JEQ) 8 1 V(S*NEQ+3+JEQ)=V(NEQ+3+JEQ) 9 DO 2 JTERM1=*NEQ 10 DO 9 JT=!*NEQ 11 9 V(3+JT)-V(2*NEQ+3+JT) 12 V(3+JTERM)=1.005*V(2*NEQ+3+JTERM)+1.E-30 13 CALL DERIV 14 DO 3 JEQ=l~NEQ 15 3 V(6*NEQ+3+JEQ)=V(NEQ+3+JEQ) 16 DO 109 JT=1*NEQ 17 109 V(3+JT)- V(2*NEQ+3+JT) 18 V(3+JTERM)-.995*V(2*NEGQ+3+JTERM)-1.E-30 19 CALL DERIV 20 DO 103 JEQ=lpNEQ 21 103 V(7*NEQ+3+JEQ)= V(NEQ+3+JEQ) 22 15 DO 5 JEQI=,NEQ 23 JEQTER=C(JTERM- 1 )*NEQ+JEQ 24 5 C(JEQTER)=(V(6*NEQ+3+JEQ)-V(7*NEQ+3+JEQ))/ 25 1 (0o. 010*V(2*NEQ+3+JTERM)+2. E-30) 26 2 CONTINUE 27 JMAT= 28 20 DO 17 ICOL=1,NEQ 29 DO 17 IROW=INEQ 30 JEQTER= (ICOL- )*NEQ+I RO W 31 QMAT ( I ROW I COL)=-. 5*V ( 3)*C(JEQTER) 32 IF(ICOL.EQ.IROW) QMAT(IROW, ICOL)=QMAT(IROWIC0L)+ 33 17 CONTINUE 34 CALL MINV ( QMAT, NEQ, DETERM, L, M) 35 DO 16 ICOL= 1,NEQ 36 DO 16 IROW=INEQ 37 JEQTER= (ICOL- 1)*NEQ+IROW 38 16 Q (JEQTER) = QMAT ( I ROW, ICOL) 39 D0 7 JEQ=Il3NEQ 40 V ( 3*NEQ+3+JEQ) =o 0 AMPLCT LISTING 18

41 DO 8 JTERM1lsNEQ42 JEQTER=(JTERM-1 )*NEQ+JEQ 43 8 V( 3*NEQ+3+JEQ)=V( 3*NEQ+3+JEQ)+V( 3)*QCJEQTER)* 44 1 V(5*NEQ+3+JTERM.) 45 7 CONTINUE ~46. JMAT=JMAT+ 1 ~47 ~ IF(JMAT-10)21, 19 19 ~48 21 AM=O. 0 497 DO 10 JEQ=1,NEQ ~50 IF(ABS(V(2*NEQ+3+JEQ)).LE.I.E-30) GO TO 10.~51 ~ AMA=ABS ( V ( 3*NEQ+3+JEQ)/V ( 2*NEQ+3+JEQ ) ) ~52 ~ AM=AMAX I (AM., AMA) ~53 10 CONTINUE IF(AM-V(I)) 11,11X12 55 11I IF(AM-0.-2*V (1)) 18.18,19 ~56 12 V(3)=0* 5*V(3) ~57 GO TO 20 5~8 18 V(3)=2.0*V(3) ~59 GO T0 20 ~60 19 DO 14 JEQ=I1NEQ ~611 4 V(3+JEQ)=V(2*NEQ+3+JEQ)+V(3*NEQ+3+JEQ) ~62 V(2)=V(2)+V(3) ~63 RETURN ~64 END AMPLCT LISTING, CONT. 19

I C DE SURCKUTN MINVNV C C PU RPO SE C INVERT A MATRIX C C USAGE C CALL MINV(AN, DSLUM) 6 C *~~~~ C DESCRIPTION OF PARAMETERS C A - INPUT MATRIX, DESTROYED IN COMP BY C RESULTANT INVERSE C N ORDER OF MATRIX A C D - RESULTANT DETERMINANT C L - WORK VECTOR OF LENGTH N C M - WORK VECTOR OF LENGTH N C C REMARKS C MATRIX A MUST BE A GENERAL MATRIX C 0 C SUBROUTINES AND FUNCTION SUBPROGRAMS. R MINV LIADTING 21

25 C NONE 26 C 27 C METHOD 28 C THE STANDARD GAUSS-JORDAN METHOD IS SED THE D *ETERM I NANT 29 C IS ALSO CALCULATED A DETERMINANT OF ZERO INDIC *ATES THAT 30 C THE MATRIX IS SINGULAR 31 C 32 C.................... > 33 C 34 SUBROUTINE MINVCA,NDLM) 35 DIMENSION A(1),L(1),M(M) 36 C 37 C - -, --- —*.* *,,**,,0o.O*@O e** 38 C 39 C IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS D *ESIRED, THE 40 C C IN COLUMN I SHOULD BE REMOVED FROM THE DOUBLE. PR *ECISION > 41 C STATEMENT WHICH FOLLOWS. 42 C 43 C DOUBLE PRECISION A*D*BIGAPHOLD > 44 C > 45 C THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION S * TATEMENTS~: 46 C APPEARING IN OTHER ROUTINES. USED IN CONJUNCTION WI *TH THIS > 47 C ROUTINE.P > 48 C 49 C ~~~THE DOUBLE1 *PRECISION VERSION~k MOF THI URUT INEQ0llTkC MUi

> 51 C 10 MUST BE CHANGED TO DABS. * > 52 C * > 53 C.. ~00O000 0 O0 00 > 54 C * > 55 C SEARCH FOR LARGEST ELEMENT * > 56 C * > 57 D=1.O * > 58 NK —N * > 59 DO 80 K=IN * > 60 NK=NK+N * > 61 L(K)=K * 62 MC(K)'=-K * 63 KK=NK+K * > 64 BIGA=ACKK) > 65 DO 20 J=KN > 66 IZ=N*(J-1) * 67 DO 20 I=KN * > 68 IJ=IZ+I > 69 10 IFC ABS(BIGA)- ABSCA(IJ))) I5S20.,20 > 70 15 BIGA=ACIJ) > 71 LCK)=I 72 M(K)=J IC 73 20 CONTINUE I:.. 74 C k MINV LISTING, CONT. 22

> 75 C INTERCHANGE R0WS * > 76 C > 77 J=LCK) 78 IF(J-K) 35, 35, 25 > 79 25 KI=K-N * > 80 D0 30 I=1,N * > 81 KI=KI+N > 82 HOLD=-A(KI) * > 83 JI=KI-K+J * > 84 A(KI)=A(JI).* 85 30 ACJI) =H0LD * > 86 C * 87 C INTERCHANGE COLUMNS * > 88 C > 89 35 I=MCK) > 90 IF(I-K) 45.45,38 * > 91 38 JP=N*(I-1) * 92 DO 40 J=I.N * > 93 JK=NK+J * > 94 JI=JP+J * > 95 HOLD=-A(JK) > 96 A(JK)=ACJI) * > 97 40 A(JI) =H0LD > 98 C > 99 C DIVIDE COLUMN BY MINUS PIVOT (VALUE 0F PIV0T ELEME *NT IS > IlO C0 CONTAINED IN BIGA) * MINV LISTING, CONT. 23

> 101 C 102 45 IF(BIGA) 48.46,48 103 46 D-O 0 104 RETURN 105 48 D0 55 I=I*N 106 IF(I-K) 50,55.50 107 50 IK=NK+I 108 A(IK)=A(IK)/(-BIGA) 109 55 CONTINUE 110 C 1> 11 C REDUCE MATRIX 112 C 113 DO 65 I=1.N 14 IK=NK+I 1 15 HOLD=A(IK) 116 IJ=I-N > 117 DO 65 J=IN > 118 IJ=IJ+N > 119 IF(I-K) 60,65,60:> 120 60 IF(J-K) 62,65,62 > 121 62 KJ=IJ-I+K > 122 A(IJ) =HOLD*A(KJ) +AC IJ) ~> 123 65 CONTINUE >:124 C MINV LISTING, CONT. 24

125 C DIVIDE ROW BY PIVOT * 1> 26 C * > 127 KJ=K-N 128 DO 75 J=I,N > 129 KJ=KJ+N * > 130 IF(J-K) 70,75,70 > 131 70 ACKJ)=ACKJ)/BIGA * > t132 75 CONTINUE 133 C > 134 C PRODUCT OF PIVOTS > 135 C > 136 D=D*BIGA > 137 C 138 C REPLACE PIVOT BY RECIPROCAL 139 C ]140 A(KK)=1 ~ O/BIGA > 141 80 CONTINUE > 142 C > 1143 C FINAL ROW AND COLUMN INTERCHANGE > 144 C * > 145 K=N o> 146 100 K=(K-1) * > 147 IF(K) 150150,105 148 105 I=LCK) 149 IF(I-K) 120,120,108 150 108 JON*C(K-1) MINV LISTING, CONT. 25

151 JR=N*(I-1) 152 DO 110 J=t,N 153 JK=JQ+J 154 HOLD=A(JK) 155 JI=JR+J 156 A(JK)=-A(J.I) 157 110 A(JI) -HOLD 158 120 J=M(K) 159 IF(J-K) 100O 100.125 160 125 KI=K-N 161 DO 130 I=I.N 162 KI=KI+N 163 HOLD=A(KI) 164 JI=KI-K+J 165 A(KI)=-A(JI) 166 130 A(JI) =HOLD 167 GO TO 100 168 150 RETURN 169 END MINV LISTING, CONT.. r.26

VI. COMMENTS FOR SUBROUTINE AMPLCT LISTING (F, w refer to Eq. (1-7)) Line 5 DERIV is called to evaluate F's from current w's. Line 7 The independent values, i. e., the w's JEQ = 1 to NEQ read into V(3 + JEQ) are stored in V(2*NEQ + 3 + JEQ). Line 8 The values of the derivative functions F's, JEQ = 1 to NEQ, evaluated in DERIV and passed through in V(NEQ + 3 + JEQ), are stored in V(5*NEQ + 3 + JEQ). Line 11 The w's are reset in V(3 + JT), JT = 1 to NEQ. Line 12 The current value w. only is multiplied by 1. 005 in preparation to call DERIV. Line 15 The new F's for w. = 1. 005 w. are stored in V(6*NEQ + 3 + JEQ) JEQ = 1 to NEQ Line 17 The w's are reset Line 18 The current value of wi only is multiplied by. 995 in preparation for call to DERIV. Line 21 The new F's for w. =. 995 w. are stored in V(7*NEQ + 1 1 3 + JEQ) JEQ = 1 to NEQ. Line 24 The value of the matrix [A], whose elements are.. = aFi/aww is calculated and stored column by column in the linear array C. 27

Line 27 JMAT is a counter incremented by one each time a trial solution cycle has been completed. No more than 10 successive doublings or halvings of the step size are permitted. Line 31, 32 The matrix [[I] - (h/2)[A]] is made up and stored in the two dimensional matrix QMAT. Line 34 The subroutine MINV inverts QMAT and returns the inverted elements in QMAT itself. Line 38 The two dimensional array QMAT is loaded column by column into the linear array Q. Line 43 The increment for each independent argument Awl = (n+l) w (n) = hFi E (-ij (h/2)aij) is calculated and stored in V(3*NEQ + 3 + JEQ), JEQ = 1 to NEQ. Line 47 No more than 10 step size changes and subsequent recalculations of the whole step are allowed. If JMAT > 10, the test for the maximum relative change is deleted and the results of the step are accepted. Line 48-59 The largest value of IAw&/wil is found and stored in AM; to avoid division by zero, this test is deleted for any Iwi < 1 x 10-30. The step size h is changed as follows: If IAw/wi max > V(1), divide current step size h (stored in V(3)) by 2. 0, discard the results of 28

the calculated increments Aw. and recalculate. If lwi/w.ilma < (1/5) V(1), double the current step size, discard the calculated results and recalculate. If (1/5) V(1) < A/wilw max < V(1) the calculations are considered satisfactory. Line 61 The dependent arguments, w's are advanced. Line 62 The independent argument is advanced. VII. MODIFICATIONS TO AMPLCT FOR EIGENVALUE DISPLAY A principle advantage of the scheme given here is that the eigenvalues of the matrix [A] do not need to be explicitly determined. There may be situations however wherein the actual values of the eignevalues are wanted. To do this create a storage matrix A dimensioned NEQ X NEQ. This can be easily done by including A as an additional argument in the argument list for AMPLCT and dummy dimensioning A(NEQ, NEQ) in AMPLCT. The actual dimensioning of A is then done in the main program where the value of NEQ is known. Between lines 30 and 31 in AMPLCT, load A from the array C as follows: A(IROW, ICOL) = C(JEQTER) Then at preselected intervals set by a counter (which may be added as an additional argument in AMPLCT) after a step has been successful, the eigenvalues of A can be obtained for display by performing a decomposition 29

of the matrix A into a lower triangular matrix having l's on principle 10 diagonal and an upper triangular matrix. The values on the principle diagonal of the upper triangular matrix are the eigenvalues of A. Such decompositions, done by Gaussian elimination and partial pivoting, -are usually available as library subroutines. In Ref. (11), a FORTRAN IV program that does this decomposition by the Rutishauser method is given. VIII. EXAMPLE SOLUTION To illustrate the use of AMPLCT, a complete FORTRAN IV program and the resulting computer output are given for the following equations dy~ d Y; -2000 Yl + 1000 Y2 + 1 dy2 dx 2 = Y1 -Y2 with Y1(0) - 0, Y2(0) = 0, from x - 0. to 4. The eigenvalues are Xl = -2000. 5 and XA2 = - 5 and are constant because the equations are linear. These equations are already in automorphic form since the independent variable x does not appear on right sides so that it really was unnecessary to include the third equation, dx- Y3 1., 3() = ) but doing so effectively includes the independent argument in the step size test. 30

Due to the step size test on I y/y i, the solution was not started at x = 0, but at x = 1 x 10 to avoid an unnecessarily small initial step size, since all initial conditions were zero. The 360/70 CPU time required for this example was. 991 sec for 65 steps using a rather coarse value of the accuracy parameter G(1) =. 50. The CPU times and number of steps for G(1) =. 10 and. 05 were 4. 072 sec and 8. 050 sec for 625 steps respectively. In each case about 30% of the steps were required to reach x =. 01.

> EXTERNAL DERA 3> COMMON/D/G(50) >C INPUT NEEDED: G(I)=ACCURACY, G(2)=XSTART, G(5 O) XMAX >C AND G(3)=INITIAL STEP SIZE >C WITH G(4),G(5),G(6) AS INITIAL CONDITIONS > NAMEL I ST/INPUT/G > DIMENSION QMAT(3,3),L(3),M(3) > I READ(5, INPUT) > WRITE(6,110) G(1),G(2),G(3).,G(4),G(5),G(6) > 110 FORMAT "1 ACCURACY PARAM... G( 1 )=', PE15 6/ > t ' INDEPENo VARIABLE. G(2)=',E15.6/ > 2' STEP SIZE, G(3)=',E15.6/ > 3' DEPEN. VARIABLE. G(4)= 'E15.]6/ ~ 4' DEPENo VARIABLE, G(5)= ',EIS.6/ > 5 ' DEPENDENT VARIABLE G(6) ' E15.6) > WRITE(6. 101 ) > N=37 > WRITE(6 100) G(2). G(4)G(5),G(6) > 100 FORMAT(' '.1P4E16.6) > 101 FORMAT( ' I 'XX'. 1 5X, 'Y( ) ',12X, 'Y(2) ', 1 2X, 'Y(3) ',/) 6 CALL AMPLCT(GNDERA, QMATL,,M) WRITE(6, 1 00)G(2), G(4) GC5), G(6) ~ rIFCG(2).LE.G(50)) GO TO 6 END SUBROUT INE DERA COMMON/D/G(50) > G(7)=-2000o*G(4)+1000.*GCS)+1 > G(8)=G(4)-G(5) > G(9)=I. ~> RETURN > END EXAMPLE PROGRAM LISTING > &INPUT G=-50 1 E-3,.01,5.0012E-4,2. 4993E.-7,1E-3,G(50)=4. ~ &END INPUT DATA FOR EXAMPLE PROGRAM 32

> ACCURACY PARAM., G()= 5.O000000OOOE-01 * INDEPEN* VARIABLE, G(2)= 9,.999999E-04 ' STEP SIZE G(3)= 9. 9999998E-03, DEPEN. VARIABLE, G(4)- 5.001200E-04 * DEPEN. VARIABLE, G(5)= 2. 499300E-07 DEPENDENT VARIABLE G(6) 9.99999E-04 > X Y(I) Y(2) Y(3). 9.999999E-04 5.001200E-04 2.499300E-07 9.999999E-04 l 156250E-03 5.00126SE-04 3.280294E-07 1 156250E-03 I. 312500E-03 5.001419E-04 4*061180E-07 I. 312500E-03 1. 468749E-03 5.001635E-04 4.841975E-07 1. 468749E-03.624999E-03 5,001901E-04 5.622686E-07 1.624999E-03 > 781249E-03 5.002199E-04 6.403317E-07 1 781249E-03 1. 937499E-03 5. 002522E-04 7. 183875E-07 1. 937499E-03 2.093749E-03 5.002865E-04 7.964362E-07 2.093749E-03, 2. 406249E-03 5.003582E-04 9.525139E-07 2.406249E-03 2. 718749E-03 5. 004329E-04 1. 108565E-06 2. 718749E-03 3.031249E-03 5.005093E-04 1 264591E-06 3.031249E-03 3. 343749E-03 5*005863E-04 1.420593E-06 3. 343749E-03 3.656249E-03 5.006639E-04 1. 576569E-06 3.656249E-03 4.~281245E-03 5.008194E-04 1.888448E-06 4.281245E-03 4. 906245E-03 5.009754E-04 2.200231E-06 4.906245E-03 5.531244E-03 5.011314E-04 2.511916E-06 5*531244E-03 6. 156243E-03 50.12872E-04 2. 823504E-06 6. 156243E-03 6781243E-03 5.014429E-04 3.134995E-06 6*781243E-03 8031242E-03 5.017542E-04 3.757686E-06 8.031242E-03 9.281240E-03 5.020658E-04 4.379988E-06 9.281240E-03 1.053124E-02 5.023766E-04. 5S001902E-06 1.053124E-02 1. 78124E-02 5.026872E-04 5.623426E-06 1. 178124E-02 1. 303124E-02 5.029982E-04 6.244561E-06 1.I 303124E-02 1.553123E-02 5.036192E-04 7.485672E-06 1.553123E-.02 1.803123E-02 5 042385E-04 8.725228E-06 1.803123E-02 COMPUTER OUTPUT 33

2.053123E-02 5.'048583E-04 9.96 32411E-06 2.053123E-02 2.o 303123E-02 5.054765E-04 1. 119971 E-05 2.303123E-02 2.553122E-02 5,060937E-.04 1. 243.463E-05 2.553122E-02 3.053122E-02 5.073270E.-04 1.489984E-05 3.053122E-02 3.553122E-02 5.085573E-04 1.735889E-05 3 553122E-02 4.053122E-02 5.097836E-O4 1 981 182E-05 4. 053122E.-02 4.553122E-02 5. 110074E-04 2.225861E-05 4. 5531 22E- 02 5.053122E-02 5 122274E-04 2.469930E-05 5.053122E-02 6.053122E-02 5. 146612E-04. 2.956245E-05 6.053122E-.02 7.053119E-02 5. 1 70791E-04 3. 440131E-05 7.053119E-02 8. 0531 1 8E-02 5. 94889E-04 3. 921604E-05 8. 053118E-02 9.053117E-02 5. 218831E-04 4. 400678E-05 9*053117E-02 1. 00531E 2-01 5.242690E-04 4. 877364E-05 1. 00531 2E-0 1 1 2053 1 E-01 5. 289924E-04. 823614E-05 1 2053 1 1E-01 1.40531 1 E-O 5336907E-04 6. 760456E-05 1.40531 IE-01 1.605311E-01 5.383205E-04 7.687985E-05 1.60531 E-01 1.805311E-O1 5.429175E-04 8.606281 E-05 1.805311E-01 2.005311E-01 5.474603E-04 9.515442E-05 2. 005311E-01 2. 405310E-01 5564245E-04 1 1 30676E-04 2.405310E-01 2. 805310 E-01 5.651971E-0O 1 306262E-04 2. 805310E-01 3.205310E-01 5.738169E-04 1.478372E-04 3.205310E-01 3.605309E-01 5.822424E-04.647075E-04 3.605309E-01 4.005309E-01 5 905265E-04 1 ~ 812439E-04 4. 005309E-0 1 4. 805309E-0 1 6.065653E-04 2, 133440E-04 4, 805309E-01 5.605308E-01 6.220061E-04 2.441856E-04 5.6053.08E-01 6. 405308E-01 6.368090E-04 2.738182E-04 6*405308E-0O 7 7 205308E-01 6. 510646E-04 3. 022887E-04 7. 205308E-01 8. 005308.E-01 6 647294E-04 3. 296430E-04 8.005308E-01 9.605308E-01 6.905228E-04 3.811962E-04 9.605308E-01 io 1 20530E 00 7. 1 43167E-04 A.4287849E-04 1 120530E 00.280530E 00 7.362934E-04 4.727137E-04 1.280530E. 00 I. 440530E 00 7.565685E-04 5. 1 32644E-04 1.440530E 00 1.600530E 00 7.752941E-04 5.506966E-04 1 600530E 00 1*920529E 00 8.085738E-04 6.172450E-04 1.920529E 00 2. 240529E 00 8. 369293E-04 6.739367E-04 2.240529E 00 2 560529E 00 8.610785E-04 7*222313E-04 2*560529E 00 2.880528E 00 8.816605E-04 7.633730E-04 2880528E. 00 3o200528E 00 8.991791E-04 7.984205E-04 3.200528E 00 3.840528E 00 9.269945E-04 8.540163E-04 3.840528E 00 4. 480527E 00 9.471201E-04 8.942788E-04 4.480527E 00 COMPUTER OUTPUT, CONT. 34

IX. COMPAR ISON TO OTHER METHODS It is of interest to compare AMPLCT with other methods devised to integrate stiff equations. Lapidus and Seinfeld7, p. 286, have tabulated the running time and accuracy of numerical experiments on several examples of stiff differential equations. Their results were obtained with an IBM 7094 while the results reported here for AMPLCT were obtained using an IBM 360/67. IBM 7094 and 360/67 machine times are not directly comparable, so to make the run time comparisons meaningful the results for each example on each machine have been normalized by their respective CPU times. required for the classical fixed step 4t order Runge Kutta technique using the same step size over the same range of independent variable. If only a portion of the Runge Kutta solution was computed due to excessive run time, the full run time was estimated on the basis that the time per step was approximately constant. Two of Lapidus and Seinfeld's four examples were chosen, the example numbers being theirs. Example 1 y' = -200(y - F(x)) + F'(x) where F(x) = 10 - (10 + x) e-x with y(0) = 10. The solution is desired between x = 0 and x = 15. The exact solution is - 200x y = F(x)+ 1Oe 35

Example 1 has two components, the eigenvalues being -200 and -1 and is only moderately stiff and not parasitic. The solution component e decays very rapidly compared to e x Example 4 y -.041 + 1 X 10Y2Y3 y7 2 y =04Y1 - 1X10 Y2Y3 - 3x 10 Y y= 3 x 107 2 with y1(0) = 1, Y2(0) = 0, Y3(0) = 0. The solution is wanted between x = 0 and x = 40. The eigenvalues are X1 = 0 and, 2 4 7 7 11 2 X + (. 04+ lxl104y + 6x107Y2) A+ (. 24x10 Y2+ 6x0 y) 0 At x = 0 the eigenvalues are (0, 0, and -. 04) while for large x they become (0, 0, and -1 x 104 ) with the absolute value of the largest eigenvalue changing from. 04 to 2405 in the range x = 0 to. 02. The system is fairly stiff. The various implicit methods use the following basic equation to advance the step Yn+l = n + Wlkl + W2k2 36

where kl = h[f(yn) + al A(Yn) kl] k2 h[ f(Yn + bl kl) + a2 A(Yn + Clkl) h2] for the equation y' = f(y) with A(y) = af/ay being either a scalar or the Jacobian matrix. The values of al, a2, bl, cl, wI and w2 depend on the particular method. The trapezoidal rule uses a1= 1/2, w1=1, w2 =0 Some of the better methods cited by Lapidus and Seinfeld were: RK4 Classical, fixed step fourth order Runge Kutta, used for the time normalization. TR Trapezoidal procedure in which from yn at xn, the trapezoidal rule is used twice to give two successive values y n+l' Y2 and then the first new value Yn+1 is replaced by (Yn + 2Yn1 + Yn2)/4 TR-EX Trapezoidal rule using the input filtering described above and global extrapolation. In global extrapolation after integration over the whole interval, the step size is decreased by a factor of 2 and the whole problem is redone. The values at a particular point were then obtained by Romberg extrapolation to the limit after three such cycles of halving. 37

CAL Calahan's method using al = a2 =. 788675134, b1 = -1. 15470054, cl = 0, w1 =3/4, w2 = 1/4. The characteristic root.t is given by 1 -. 578 hX -. 45622 (hX)2 1 - 1. 578 hX +. 622 (hA) LW1 Liniger and Willoughby's method using exponential fitting. The step is advanced using Yn+1 Y + h[(l - A)Y + BY] which has a characteristic root 1+ P3hX 1 - (1 - ( ) hX The constant 1 is determined by requiring,p to approximate the exponential p(hX) = ehX which yields 1. 1 hAX -hXe -1 The largest negative eigenvalue is chosen to make the exponential fit. The backwards Euler method uses j3 = 0 which corresponds to a fit at hA - -0. ARK A classical fixed step fixed step fourth order Runge-Kutta method, set up so as to use the same calling conventions as AMPLCT and used for the time normalization. 38

The numerical results using these methods are summarized in Tables 1-4. The relative accuracy Rn is given by IYn - Y(xn) I/Y(xn) where yn is the computed value and Y is the true value. The errors are shown at two values of x. The CPU times, normalized by the classical fourth order Runge-Kutta method using the step size noted, are given in the last column. For example 4, the solution obtained with RK4 or ARK using a step size h =. 001 was assumed exact. The parameter V(1) noted with AMPLCT is the relative accuracy parameter which controls the step size. In addition to the previous methods, a recent stiff method developed 7y 13G 14 by Gear' 1' was run with the 360/67 using example 4. Gear's method is a predictor and iterated corrector method in which both the order and step size are automatically selected to give the largest step and hence reduce the overall calculation time. A FORTRAN IV listing of the alogrithm is given in Ref. (14). The result of Gear's method is given in Table 5. Table 2 shows the almost linear increase of computing time with ARK that occurs with decreasing step size. The largest relative error in example 1 occurs near x -=. 016 for all methods. Even though a step size of h =.01 was small enough so that a stable solution was possible using ARK, the accuracy at x =. 016 is unacceptable for this step size. For example 1, AMPLCT appears to offer no marked advantage over the stiff methods cited by Lapidus and Seinfeld. 39

Step Size Relative Accuracy(Rn) Relative Time Method h x =.4 x = 10. t/K4 RK4.01 1. 0-5 2. 09 1. 0 TR.2 1.85-2 4. 3-5.182 TR-EX.2 1.44 1. 08 3. 26 CAL.01/. 2** 1. 7-2 4. 0-8.091 LW1.2 1. 1 5.08.273 TABLE 1. Example 1 Accuracy and Relative Run Times from (7). tRK4 = 11 sec to x = 15 using h =.01. * 1.0o-5.1.0x 10-5 ** h changed from. 01 to. 2 at x =.1 40

Step Size Relative Accuracy (Rn) Relative Time Mlethod h x=.016.4 10 t/t.001 1.91 <1-7 <1-9 8. 85.002 7. 04 <4.65 ARK 2-7 -K.005 4.662 6. 55 < I 1. 95.010 2. 57 2. 35-5 1. 03-9 1. 0 AMPLCT V(1) =. 05 1. 13- 1. 14- 2. 98-.713 AMPLCT 4 5 7 v(1) =. 10 2.61 4.58 5.72.332 TABLE 2. Example 1 Accuracy and Relative Run Times Using ARK and AMPLCT, tARK = 12. 2 sec to x = l5 using h =.01 41

Step Size RIn R 2n R3n Relative Time Method o h | x=.4 10.4 10.4 10 t/tRK4 RK4.001 0 0 0 0 0 0 1.0 TR.2 1.35-3 1.05 2.12-1 2. 4-1 9.0 -2 1.5-2.0674 TR-EX.2 1. 72-5 3. 6-4 3. 5 4. 3 6. 8-4 1 2.246 CAL.005/. 02* 2.4-3 1. 01-1 2. 5+0 6. 01 1. 62-1 5.4-1.0725 LW1.02 1. 6 4. 9- 2. 4 1. 3 3. 2-3 4. 4.145 TABLE 3. Example 4 Accuracy and Relative Run Times from (7). tRK4 = 138 sec to x = 40 using h =.001 (estimated). * h changed from.006 to.02 at x = 1. 42

R Step Size elative Accuracy Rel. Time Step Size Method R n R2n R3n t/tARK h x=.40 7. 40.40 7. 40.40 7. 40 ARK.001 0 0 0 0 0 1.0 x =.40 7. 40.40 7. 40.40 7. 40 AMIPLCT 4. 53. 3- -3............. 00933 V(1)=.05 6.68x104 6.85x10'5 3. 50x10 1. 12x10-3 4.58x10-2 2.15x10 x=.645 7. 34.645 7. 34.645 7. 34 AMPLCT 00516 V(1) =.10 2. 59-2 6.42x10-4 1.84-2 2.68x10 3 1. 111 6. 68 x 103 TABLE 4. Example 4 Accuracy and Relative Run Times using ARK and AMPLCT. tARK = 382 sec to x = 40 using h =.001 (estimated). 43

Relative Accuracy Rel. Time Method R t/t In ~~~2n R3AR x=.40 7. 4.40 7. 4.40 7. 4 GE ARP- 4 4 4 2 4 EPS = 1.9.1 4.5 4. 95 1.4 2. 94 1.69 h I x 10 TABLE 5. Example 4 Accuracy and Relative Run Time of GEAR's Method. ARK = 382 sec to x = 40 using h =.001 (estimated) Error control, Euclidian norm of estimated relative errors for all components per step. initial step size, automatic step size contnil used. 44

In example 4 AMPLCT shows up quite well compared to the methods discussed by Lapidus and Seinfeld even considering that ARK took about 3 times longer than RK4 to integrate this example. It was found that. the run time of AMPLCT depends almost entirely on the relative accuracy parameter V(1) and is independent of the initial step size. In example 4 the errors at 7. 40 and 10 were similar; the point of error evaluation at x = 7. 40 being a matter of convenience. The results of RK4 or ARK using h =. 001 were taken as exact for example 4. ARK and RK4 both become unstable using h =.001 for x > 16 where maIXm 2780. The generally higher accuracy of TR-EX has been obtained at the expense of an increased run time. Gear's method is faster than AMPLCT by about a factor of 5 for example 4 using EPS =. 0001. The order varied between 1st and 2nd up to x =.70 at which point it selected 3rd order for the remainder of the problem. Using EPS =. 001, Gear's method was unstable. 45

REFERENCES 1. M. P. Sherman, '"Radiation Coupled Chemical Nonequilibrium Normal Shock Waves, " J. Quant. Spect. Rad. Transf. 8, 1968, p. 569. 2. H. Lomax, H. E. Bailey, "A Critical Analysis of Various Numerical Integration Methods for Computing the Flow of a Gas in Chemical Nonequilibrium, " NASA TN D-4109, 1967. 3. H. E. Bailey, "Numerical Integration of the Equations Governing the One Dimensional Flow of a Chemically Reactive Gas, " Phys. Fluids 12, 1969, pp. 2292-2300. 4. H. Lomax, "On the Construction of Highly Stable Explicit, Numerical Methods for Integrating Coupled Ordinary Differential Equations with Parasitic Eigenvalues, " NASA TN D 4547, April 1968. 5. H. Lomax, "Stable Implicit and Explicit Numerical Methods for Integrating Quasi-Linear Differential Equations with Parasitic-Stiff and Parasitic-Saddle Eigenvalues, " NASA TN D 4703, July 1968. 6. H. Lomax, H. E. Bailey, F. B. Fuller, "On Some Numerical Difficulties in Integrating the Equations for One Dimensional Nonequilibrium Nozzle Flow, " NASA TN D 5176, April 1969. 7. L. Lapidus, J. Seinfeld, Numerical Solution of Ordinary Differential Equations, Academic Press, 1971. 8. L. Fox, Numerical Solution of Ordinary and Partial Differential Equations, Pergammon Press, 1962, Ch. 4. 9. S. W. Bowen, C. Park, "Computer Study of Nonequilibrium Excitation in Recombining Nitrogen Plasma Nozzle Flows, " AIAA Piper 70-44, New York, 1970. Also, AIAA J., 9, No. 3, March 1971, p. 493. 10. A. Ralston, A First Course in Numerical Analysis, McGraw Hill, 1965, p. 510. 11. B. Carnahan, H. A. Luther, J. O. Wilkes, Applied Numerical Methods, Wiley, 1969, pp. 236,249. 46

12. M. E. Fowler, R. M. Warten, "A Numerical Integration Technique for Ordinary Differential Equations with Widely Separated Eigenvalues," IBM J. of Res. and Dev., 11, 1967, p. 537. 13. C. W. Gear, "The Automatic Integration of Ordinary Differential Equations," Comm. ACM, 14, 1971, pp. 176-179. 14. C.W. Gear, "Algorithm 407, DIFSUB for Solution of Ordinary Differential Equations," Comm. ACM, 14, 1971, pp. 185-190.

UNIVERSITY OF MICHIGAN 3 9015 02519 623211111 I!9d1L 02519 6232