LIBRARY OF THE 
 
 UNIVERSITY OF ILLINOIS 
 
 AT URBANA-CHAMPAIGN 
 
 510.64 
 
The person charging this material is re- 
 sponsible for its return to the library from 
 which it was withdrawn on or before the 
 Latest Date stamped below. 
 
 Theft, mutilation, and underlining of books 
 are reasons for disciplinary action and may 
 result in dismissal from the University. 
 
 UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN 
 
 DEC 1 4 WT6 
 
 |EC 1 2wEC0 
 
 IP i o T .|82 
 
 AUG 2 9 1997 
 
 L161 — O-1096 
 
A) 
 
 ■2- 
 
 UIUCDCS-R-73-575 
 
 /yuJC-i 
 
 coo-11469-0225 
 
 f 
 
 DOCUMENTATION FOR DFASUB— A Program for the Solution 
 of Simultaneous Implicit Differential and Nonlinear Equations 
 
 ty 
 
 R . L . Brown 
 C. ¥. Gear 
 
 July 1973 
 
 DEPARTMENT OF COMPUTER SCIENCE 
 UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 
 
 URBANA, ILLINOIS 
 
 THE LIBRARY OF THE 
 
 AUG 6 1973 
 
 UNIVERSITY OF ILLINOIS 
 AT URRArw r»"*MOAiGN 
 
UIUCDCS-R-T3-5T5 
 
 DOCUMENTATION FOR DFASUB--A Program for the Solution 
 of Simultaneous Implicit Differential and Nonlinear Equations* 
 
 by 
 
 R. L. Brown 
 
 C. W. Gear 
 
 July 1973 
 
 DEPARTMENT OF COMPUTER SCIENCE 
 
 UNIVERSITY OF ILLINOIS AT URBAN A- CHAMPAIGN 
 
 URBANA, ILLINOIS 6l801 
 
 Supported in part by the Atomic Energy Commission under contract 
 US AEC AT(ll-l)lU69. 
 
Digitized by the Internet Archive 
 in 2013 
 
 http://archive.org/details/documentationfor575brow 
 
TABLE OF CONTENTS 
 
 Page 
 
 1. INTRODUCTION 1 
 
 2. THEORETICAL BACKGROUND 2 
 
 3. REPRESENTATION OF VARIABLES; SUBROUTINES 10 
 
 k. USER PROGRAM FOR DFASUB 15 
 
 5. PROGRAM LOGIC 19 
 
 LIST OF REFERENCES 25 
 
 APPENDIX I 
 
 SAMPLE PROGRAM 26 
 
 APPENDIX II 
 
 DFASUB LISTING 30 
 
1. INTRODUCTION 
 
 The purpose of the FORTRAN program DFASUB is to integrate a 
 system of equations — described "by a mixture of ordinary differential 
 equations, nonlinear equations, and linear equations — using a discrete 
 variable method. The equations are written in the form 
 
 f(y_» y_' , t) + Pv = (1.1) 
 
 where the vectors y_ and y_' = dy_/dt are of length p, v is of length 
 q-p , P is a q x (q-p) matrix, t is the independent variable, and f is 
 a vector function of length q. 
 
 DFASUB uses the variable values to integrate the system using 
 a multistep method designed for use with stiff ordinary differential 
 equations which also works well with non-stiff equations and which will 
 handle nonlinear equations. This paper deals with both the theory behind 
 the integration method and the program itself. Section 2 presents 
 the theory. 
 
 At present, the routine calls a number of other subroutines which 
 are compiled by a general -purpose system [l]. The function of these will 
 be described so they can be replaced by user-supplied FORTRAN subroutines. 
 The form of the internal variables of the program is described in Section 3 
 to allow the user to write subroutines equivalent to those above as described 
 in that section. Section h details calling parameters to allow a user to 
 write a complete integration package built around DFASUB. Section 5 describes 
 the detailed program logic of DFASUB. 
 
 This introduction, together with Sections 3 and k, serve as a 
 user's guide to DFASUB. Reading Section 2 is sufficient to understand the 
 theory behind the routine. Someone interested in the programming aspects 
 of the routine may read Sections 1, 3-5 to understand how DFASUB works. 
 
2. THEORETICAL BACKGROUND 
 
 Suppose that the values of y_ are known (approximately) at a 
 
 number of equally spaced points t . , i = l,2,...k, where t. jn > t.. 
 
 n-i l+l 1 
 
 The backward differentiation formula gives the relation 
 
 h K = ~ 3^ ( Vn + a i^n-l + - ■ ' + Vn-k } 
 
 where h = t , - t. > 0. (The coefficients a. and B„ can be found in 
 i+l l l 
 
 Gear [2], p. 217.) If this is substituted in (l.l) for t = t we get 
 
 (2.1) 
 
 F ( 
 
 ) = t(y 9 - 
 
 n ^-n' -i. 
 
 
 + E,t)+Pv =0 
 n n n— n 
 
 (2.2) 
 
 where 
 
 E = - -=— (a y +, 
 
 n h3 Q Pti-1 
 
 + \ ^n-k } 
 
 is known. Hence, (2.2) is a system of q equations in the q unknowns 
 y and v . In general, they are nonlinear. If these are solved by i 
 Newton-like iteration, we get 
 
 where 
 
 ^n,(m+l) 
 
 v 
 
 -n,(mKL) 
 
 ^,(m) 
 
 -n,(m) 
 
 r-1 
 
 J " F n ( ^n,(m)' ^n,(m) 
 
 3F 
 
 J = 
 
 3(i,v) 
 
 3f 
 
 a ° £ i p 
 
 3y_ h6 n 3v_' , 
 
 (2.3) 
 
 (2.U) 
 
 and y , » and v , v are the q iterates for the solution y and v of (2.2) 
 ti,[i) -n,(m) "hi -n 
 
 If accurate initial guesses y ,- N and v ,< s are used, very few iterations 
 
 -n,(0) -n,(0) 
 
 of (2.3) are necessary. In fact, accuracy is not needed in v 
 
 as we 
 
 -n,(0) 
 
 show below that the iterates y ,_, N and v ,_, s are independent of v /_.v 
 
 *4i,(l) -n,(l) -n,(0) 
 
 if round-off errors are ignored. To see this, compare the first iterate 
 
y , N and v ,_x with the first iterates 
 •41,(0) -n,(0) 
 
 f / v and v / x starting with 
 
 y /_, x and v ,_ * starting with y , _ N and v , n * . From (2.3) 
 «4i,(l) -n,(l) s ^.(O) -n,(0) 
 
 z = 
 
 4,(1) "Zn.CL) 
 
 ^,(1) " ^n,(i: 
 
 o 
 
 •4i 9 (0) 
 
 _: 4i,(0) 
 
 - J" 1 P (v 
 
 - V 
 
 n v n,(0) '" -4i,(0) 
 
 ) 
 
 Hence , 
 
 Jz = J 
 
 ^1,(0) 
 
 ^1,(0! 
 
 - P n ( V(0) V(0) } 
 
 3f 
 
 a 31 
 
 3y_ hS Q 3y_' 
 
 P - P 
 n n 
 
 o 
 
 j4i,(0) " ^n,(0)_ 
 
 = 
 
 Since we assume J is non -singular (or the Newton iteration will not work) , 
 
 we see that z_ = 0; hence, f , v = y /_ \ and v , <. = v , s are independent 
 
 of v /«\. Note that the only condition on J for this to be true is that 
 -n,(0) 
 
 the right-hand p columns of J be exactly P . The left-hand q-p columns can 
 contain anything, so we need only approximate J in those positions. 
 
 Initial Guess for y , >. 
 
 Since past values are known, a good initial guess can be found 
 using polynomial extrapolation. We assume that we know y . , i <_ i <_ k. 
 d y" . (The latter was found at the last step or could be evaluated 
 
 an 
 
 from (2 . $ with n-1 replacing n.) With those values we can use a Hermite 
 interpolation formula in the form 
 
Since equal intervals are used, the a. and 3 are independent of n and h. 
 
 The Basic Algorithm 
 
 The initial guess is calculated from (2.5) — it is called the 
 
 predicted value, v ,„ N is set to v _, . Equation (2.3) is iterated 
 * -n,(0) -n-1 
 
 until v / \ and v / \ appear to have converged. (if they do not, the 
 ■Hi,(m) -n,(m) 
 
 step size h is reduced as discussed later and the process is repeated.) 
 
 The error introduced by one step of (2.1) is known to be of 
 
 the form 
 
 ,k+l ( k+1 ) 
 
 k+1 K ^ J 
 
 If the error in the numerical solution is a sufficiently smooth function 
 
 of h and t, h y (^-i) + '-'(h can ^ e computed using a difference 
 
 formula where £ is not necessarily the same point as £ , but both £ and 
 
 £ are in the interval (t , , t ). This enables the error to be estimated 
 
 n-k n 
 
 under the assumption that y changes slowly over the interval. 
 
 If the error estimate is too large, the result is rejected, and 
 a smaller step is tried. Otherwise, the step is accepted and a suitable 
 step size and order (k) is chosen for the next step. 
 
 Internal Representation of the Algorithm 
 
 We have used a Nordsieck [3] vector to represent past values of 
 the data because of the ease of order changing and error estimation. We 
 
 will describe this vector below for a single variable y and its associated 
 
 B — n 
 
 past values y ,...y , . This allows us to use the vector notation 
 n-1 n-k 
 
^ = [y n> hy n' K-l'-" y n-k + l ] 
 
 for the set of saved information concerning y at t . (T is the transpose 
 
 operator.) When we step from t to t , we use (2.5) to predict y , _. N 
 
 n-1 n n,(0j 
 
 from y At the same time, we save the values of y , . . . .y , , from 
 
 J 4i-1 n-1 n-k+1 . 
 
 y . We can represent this process by 
 
 n,(0) 
 
 n-1 
 
 n-2 
 
 n-k+1 
 
 6 1 a 2 
 
 
 
 
 
 
 
 
 Vi 
 
 \ 
 
 
 
 
 
 
 
 
 
 
 
 
 
 i 
 
 
 
 
 .. 
 
 *n-l 
 
 (2.6) 
 
 Next we iterate (2.3). To do this we must substitute y , N and v / >, 
 
 " , n,(in) -n,(m) 
 
 into F (y , v ). Thus, from (2.2) we must evaluate 
 n n n 
 
 a 
 
 F (y , s , v , v ) = f(y ( v, - — — y , x + E , t ) + P v . . (2.7) 
 n n,(m) n,(m) J n,(m) h3 n,(m) n n n n,(m) 
 
 We will normalize (2.l) so that a = -1. Let us write 
 
 
 
 J n,(m) 8 Q n,(m) n 
 
 =hy n,(m-l) + 3^ (y n,(m) " ^.(m-lp if m > ° 
 
 = + 37 [a i y n-l + Vn-2 + "- + Vn-k + hB l hy n-l ] 
 
 (2.8) 
 
 ~-t [a i y n-l + --' + Vn-k ] if m= ° 
 
Hence, 
 
 hy n,(0) = " 
 
 a k" a k P l 
 
 n-1 3 Q y n-k B Q y n-l 
 
 = yv + . . . + y, y , + 6 ^ hy ' n 
 'rn-1 ■ Vn-k 1 °n-l 
 
 Thus, (2.7) can be written as 
 
 F n (jr n,(m)> v n,( m )> = f(y n,(m)' h K,U)' V + Vn,(») 
 
 If ve vrite Aj(m) = [y n>(m) . ^ j(m ). Vl y n-k + l ]T ' ve see that 
 
 we can combine (2.9) with (2.6) to get 
 
 (2.9) 
 
 (2.10) 
 
 ^,(0) 
 
 TT 3 1 a 2 
 
 Y l 6 1 Y 2 
 
 10 
 
 
 
 
 a k-l a k 
 
 Y k-1 Y k 
 
 
 
 
 
 Zn-1 = %-l 
 
 (2.11) 
 
 This is called the predicted value of y. Now we note from (2.8) that 
 
 2n',(mM) = *n,(m) + ^ (y n,(irH-l) " y n,(m) } (2 ' 12) 
 
 1 T 
 
 where c = [l, — , 0,. . . , 0] and y / . -, \ - y / \ is given by a component 
 — ^n n,(m+lj 'n,(mj 
 
 of (2.3). 
 
 The final step is to perform the transformation to a Nordsieck 
 
 vector. The approximations y , y' , y _,,... y , ., determine a unique k-th 
 ^ J n' J n' °n-l' J n-k+l 
 
 k (k) 
 
 degree polynomial. Let its scaled derivatives at t he y , hy' ,... h y ' /k! 
 
 n n n n 
 
 There is a unique non-singular k+1 by k+1 matrix Q such that if 
 
 ^[ V ^>, y M 2 >/ 2 *<*w 
 
then 
 
 z = 
 — n 
 
 %L- 
 
 With this we restate (2.1l) as 
 
 an 
 
 d (2.12) as 
 
 where 
 
 and 
 
 Vo 
 
 ^1,(0) 
 
 z 
 — n 
 
 = QBQ z , = Az _ 
 -n-1 -n-1 
 
 ,(m+l) " QZ n,(m+l) 
 
 = ^n,(m) + I ( ^n,(m+l) 
 
 I = Qc 
 
 ^n.dn) 
 
 A = QBQ 
 
 -1 
 
 1 
 
 1 
 
 1 
 
 1 
 
 1 1 
 
 
 
 1 
 
 2 
 
 3 
 
 U k 
 
 
 
 
 1 
 
 3 
 
 6 
 
 
 
 
 
 1 
 
 
 1 
 
 o \ k 
 \ 
 \ 
 
 \ 
 
 1 
 
 (2.13) 
 
 (2.1U) 
 
 The fact that A is a Pascal triangle matrix follows from the fact that the 
 predictor performs a k-th order exact approximation z , v to z . The 
 
 advantage of this form is that the step size can "be changed easily by 
 
 2 k 
 
 multiplying the vector z by C(a) = diag[l, a, a ,..., a ] where 
 
 a = h /h , .. The order can be reduced by dropping a column from A, a 
 new old J ** ° 
 
row from A and z , and changing l_. The order can be increased by adding 
 
 — n 
 
 k+l (k+l) 
 a row and column to A, changing &_, and estimating h y /(k+l)!. This 
 
 is estimated by the last element of (z - z _.)/(k+l). (This is also 
 J -n — n-1 
 
 used in computing the error estimate.) The values of l_ can be found in 
 Gear [2], p. 217- 
 
 The choice of step size and order depends on several things. 
 First, if the error estimate after the corrector converges is too large, 
 the step size is reduced to h/k and the step is taken over; if this fails 
 to work after three such attempts, the order is reduced to 1 and a final 
 attempt is made before giving up. 
 
 If the convergent corrector value is within error bounds , then 
 every k steps at order k the step size at the original order, order k-1, 
 and k+l is calculated with a bias toward the lowest order (with the least 
 work required) to get the largest step size. This is only done every k 
 steps due to stability considerations [h]. 
 
 The error estimator 
 
 h "4i M _ n ,i ,, n-.i ,i n 
 
 j_0 J n_ J J n ~J 
 
 is the local truncation error for the method in use if a = -1. If 
 y(t) has a continuous (k+l)-th derivative, then 
 
 L h (y(t)) = C k+1 y (k+1) U) h k+1 
 
 if the method is of order k and t-h < 5 < t. We can define 
 
v* k+1) 
 
 c. 
 
 k+1 h k+1 (k + l)! 
 
 since for a k-th order method applied to the (k+l)-th order polynomial 
 t , the error term L (t ) is known and 
 
 (t w ) (ktl) - (k*l)l 
 
10 
 
 3. REPRESENTATION OF VARIABLES; SUBROUTINES 
 
 DFASUB requires several variable arrays for integrating a 
 system; these arrays are described below. The remaining arguments in 
 the calling sequence are briefly described following those. For 
 notational convenience N is the number of equations q, NY is the number 
 of nonlinear equations p, and NL = N-NY. The notation A(b)C means 
 "from A to C in steps of B" . 
 
 Y(J,I) for I = l(l)NY and J = l(l)T contains the current value 
 of the differential and nonlinear variables in its first row. Each 
 subsequent row J contains the (J-l)-th derivative of the variable y. 
 multiplied by H**( J-l)/ ( J-l) ! . Thus, Y contains all the elements 
 necessary for a Taylor series expansion without requiring the derivatives 
 to be multiplied by H /k! before use. Thus, the predictor step consists 
 of multiplying Y(l,J) by the Pascal triangle matrix as described in 
 Section 2. 
 
 YL(l) for I = l(l)NL holds the value of the linear variables, 
 corresponding to v in Section 2. 
 
 DY(l) holds the result of calling the evaluation subroutine 
 DIFFUN(T,G,DY,Y,YL,HINV) ; for I = l(l)NY it contains the difference 
 y! - f(y.) for y. = Y(l,l), y! = Y(2,l)*HINV where HINV is the inverse 
 of the step size H. For I = NY+l(l)N the error in the linear equations 
 stored in YL(l) is in DY(l). Thus, DY(l) is essentially the correction 
 to the values in Y(l,l) and YL(j) before the Newton iteration. 
 
 SAVE(J,I) for J = 1(1)7, I = l(l)NY holds the initial values 
 of the nonlinear variables at the beginning of an integration step. If 
 the step fails, the independent variable T(l) is restored to its original 
 values and Y(j,l) takes on the values in SAVE(J,l). YLSV(l) for 
 I = l(l)NL serves the same purpose for the linear variables in YL(l). 
 
11 
 
 ERROR(l) for I = l(l )NY stores the sum of the corrections 
 computed during each integration step. It is initialized to zero after 
 the predictor evaluation and is incremented by Fl(l) — the improved 
 corrector — after each corrector step. If the integration step is 
 successful, all elements of Y(J,l) are updated by Y(j,l) = Y(J,l) + 
 A(J)*ERR0R(I). The difference for ERROR(l) between two integration steps 
 is used to calculate the next higher derivative element when the best 
 step size for order k+1 is being determined. When it is known that on 
 the next step a new value of H will be calculated, ERROR(l) is saved in 
 ERSV(I) for I = 1(1)NY. 
 
 YMAX(I) for I = l(l)NY stores the maximum of 1.0 or the 
 largest value of |Y(l,l)| computed. YMAX(l) is used in the error 
 evaluation for an absolute error test when Y(l,l) has never exceeded 
 i.O and a relative error test when it has. 
 
 A(j) , J = l(l)NQ+l, stores the updating factors for the higher 
 order derivative terms in Y. It corresponds to the vector 2_ in Section 2, 
 
 G(J) is storage for global variables, parameters that may be 
 easily changed from one integration to the next. They allow a user to 
 experiment with different parameters in attempts to optimize given 
 systems under simulation. 
 
 T(J) contains variables that are dependent only on G(*) and 
 the independent variable T(l). Thus, they need to be evaluated only 
 for each value of the independent variable. 
 
 PW(I,J) is the matrix J. 
 
 The remaining calling arguments to DFASUB follow: 
 
12 
 
 EPS is the error per step criterion that nust be met by the 
 L2 norm of the error. 
 
 EQN and VAR are two arrays that are required by routine 
 MATSET (see following) when MATSET is compiled by the general-purpose 
 system [l]. When this is not the case, they can be dummy arrays. 
 
 HMAX is the largest step size that DFASUB will take. It 
 should be chosen less than the period of any known periodic solution. 
 Otherwise, no special care need be taken in choosing it. 
 
 HMEN is the minimum step size taken. H is the initial step 
 size used unless it causes convergence problems , in which case a new 
 step size is chosen. 
 
 JSTART is a start indicator. When it is zero on a call to 
 DFASUB, the routine initializes itself. When it is 1, it assumes that 
 initialization has taken place and continues integration. On exit, 
 it contains the current order. 
 
 KFLAG is the output completion code. Its meanings are 
 
 +1 successful integration 
 
 -1 error test failure for H > HMEN 
 
 -2 too many floating point errors 
 
 -3 the corrector failed to converge for H > HMIN 
 
 -k the corrector failed for even the first order method 
 
 It is possible to write into the main program a routine to 
 try possible remedies when DFASUB quits with an unsuccessful completion 
 code making use of KFLAG and the computed GO TO in FORTRAN. After these 
 remedial measures are taken, JSTART should be set to 1 and DFASUB called again. 
 
13 
 
 M is the dimension of DY , normally the same as N. 
 
 MAXDER is the maximum order method used, never more than 6. 
 
 ML is the number of variables in Y(l,J), J = l(l)ML, which are 
 included in the error test. If the error in any group of variables is 
 not critical, then they should be placed at the end of the Y array and 
 ML set to exclude them in the error test. N is the total number of 
 equations. NL is the number of linear equations. 
 
 TEND is the final value of T(l). Integration stops when this 
 value is reached. 
 
 Subroutines 
 
 The dependent variables in DFASUB fall into three classes : 
 those dependent only on global terms or parameters in the G(*) array, 
 variables described by linear equations which have Jacobian elements 
 dependent at most on the arrays G(*) and T(*) so these are the elements 
 of the matrix P of Section 2, and variables whose entries to the 
 Jacobian matrix change whenever some other variable changes its value. 
 To handle these differences efficiently, five subroutines called by 
 DFASUB are concerned with computing the Jacobian matrix and performing 
 the corrector iteration (2.3). 
 
 Since the elements in the Jacobian have different dependencies, 
 they can be recomputed at different times by 
 MATSET (A(2) ,DY ,EPS ,EQN ,G,HINV,0,M,MF,N ,NY ,IT,PW,F1 ,T ,VAR,Y,YL) . 
 
 When IT = 1, only those values in the Jacobian that are fixed 
 for one set of parameters G(*) are calculated. This need be done only 
 once per integration. When IT = 2 , only those values that are dependent 
 on T(*) and G(*) are computed. This is done whenever a new step is 
 
Ik 
 
 started, causing a new value T(l) = T(l) + H to be computed, and 
 
 3 new O-Lu. 
 
 whenever the step size is reduced and the integration step restarted. 
 When IT = 3 , all remaining elements are computed. This must be done 
 at each corrector iteration. 
 
 The corrector iteration (2.3) does not actually multiply by J 
 but rather performs a symbolic inversion of the sparse matrix PW. Again, 
 the symbolic inversion of PW is broken up according to the dependence of 
 each element on T , G, and Y. MATINl(PW) inverts just those elements 
 that are dependent only on G. MATIN2(PW) inverts elements dependent on 
 T and G. MATIN3 inverts all remaining elements. 
 
 MATMUL(PW,DY,Fl) performs the equivalent of solving 
 (PW)(Fl) = DY for Fl; DY holds the values of F(y, v) = F(y , y' , t) + P(t)v; 
 and Fl is a vector used in the correction step. 
 
 DIFFUN(T,G,DY,Y,YL,HINV) is the differentiating subroutine which 
 places y! - f . (y_) = HINV*Y(2,l) - f(y(l,l)) into DY(l) for I = l(l)NY, 
 and the correction to the linear equations into DY(l), I = NY+l(l)N. 
 
 The remaining subroutines are not connected with the 
 differential algebraic equation system, but are used for convenience. 
 
 KNTSPl(l) is an assembly language function that keeps track of 
 floating point overflow and underflow, and returns a value of three or 
 more if enough such errors are encountered to invalidate the results of 
 th e co mput at i on . 
 
 COPYZ(S,X,L) copies the L values of single precision array X 
 into single precision array S. The actual variables used are double 
 precision arrays so L is twice the size of the actual arrays. If 
 X = Y, S = SAVE, we have L = lU*NY. For X = YL , S = YLSV, we have L = 2*NY. 
 
 S2(T,G) evaluates those variables stored in T(*) that are 
 dependent on T(l), the independent variable; and G(*), the global parameters 
 
15 
 
 h. USER PROGRAM FOR DFASUB 
 
 To write a program to use DFASUB all the subroutines called 
 by it must be present in some form. The subroutine DIFFUN , which is 
 supplied by a system compiler in the general -purpose system, must be 
 provided to solve the system 
 
 f(y s y 1 , t) + P(t)y_ = 0. 
 The calling arguments are the T(*) and G(*) arrays, the Y array containing 
 y. in Y(1,I) and H*y! in Y(2,l), HINV = l/H, and v. in the array YL. The 
 system can be put into the proper form by expressing the equations as 
 
 = y l -(f 1 ( Y » x L,T,G) + P^) 
 
 ■ y NY "( VY,YL,T,G) t ? m Z) 
 
 ° 7 W (Y ' YL ' T ' G) + P NY+ 1^ 
 = f N (Y,YL,T,G) + P N v 
 
 for P the I-th row of P(t). 
 
 The zeros are then replaced by 
 
 DY(I) , I = 1(1)N. 
 An example system is to be found in [5]. It is 
 
 = y! - s + (r-y. ) + E b. . y., i = l(l)U 
 
 i=l 
 
 ij J 
 
 for 
 
 r = (y 1 + y 2 + y 3 + y u )/2 
 
 s = E (r-y )"/2 
 i=l 
 
 [b ]-B 
 
 i+UT.5+e -1+52. 5+e -^7.5+e -52.5-e 
 
 ■i+52.5+e UUT.5+G 52.5+e U7.5-e 
 
 -HT.5+e 52.5+e UU7.5+e ^52. 5+e 
 
 -52.5-e UT.5-e ^52 . 5-e UUT-5+e 
 
 for e = .00025. 
 
16 
 
 = *5 + y l y 6 + y l y 6 
 
 = 2 H = y 6 " y i + v i " x " e_t 
 
 = v l + F 6 
 = v 1 - v 2 + y x y 6 
 
 = v x - v 2 + F T 
 
 = v + v + 5v v 
 1 2 ^1 y 2 
 
 = v 1 + v 2 + F 8 
 
 for y <°> = y<°> - 1, v[°) = -2, v^ = -3, y[ 0) = -1, i = l(l)U. 
 
 Appendix I contains the FORTRAN routine DIFFUW(T,G,DY,Y,YL,HINV) 
 that is used to evaluate this system as veil as a sample routine S2(T,G) 
 to solve T(2) = EXP(-T(l)). 
 
 The computation of the Jacobian can be split up into the three 
 subsets, each of which is evaluated by MATSET depending on the parameter 
 IT. Three different inversions are also required. If the user would 
 rather not write these four routines, a simpler routine MATSET can be 
 written so that nothing is done when IT is equal to 1 or 2 , and the entire 
 Jacobian matrix is generated when MATSET (...,3,...) is called. The 
 
 evaluation of the Jacobian can be done by a special routine that computes 
 
 9f. 9f. 
 
 - — and - — ; — according to an analytic formula, or else numerical 
 y j y J 
 
 differencing as in [6] requiring N calls to DIFFUN may be used. MATIN 3 
 can be replaced by a call to a standard Gaussian elimination subroutine, 
 while MATIN 1 and MATIN2 are made dummy subroutines. MATMUL can be replaced 
 by a back substitution subroutine written to work with the routine in MATIN3. 
 
IT 
 
 Two implementations of this procedure are possible. The 
 easiest one is to write a subroutine which has a place in its calling 
 sequence for each of the parameters in the calling sequence in DFASUB , 
 and then manipulate these parameters, possibly preparing them for a call 
 to a second subroutine. This has been done in the remainder of 
 Appendix I. For example, MATSET has dummy variables for EQN and VAR , 
 two arrays MATSET as compiled by the general-purpose system use but are 
 not needed in this version, which performs numerical differencing to 
 compute PW when IT is 3 and does nothing when IT is not 3. MATIN 3 (PW) 
 simply calls Moler's routine [7] DECOMP. Since DECOMP requires the 
 number of variables N, this has been provided in the unnamed COMMON block 
 in MATMUL and DFASUB. MATMUL calls Moler's routine SOLVE which also 
 needs N. Both of these routines communicate through an array IP, which 
 is provided in a named COMMON/IPP/ and dimensioned at least as large 
 as the system. 
 
 The second implementation requires that the user remove all 
 
 superfluous subroutine calls and change the calling sequence of the 
 routines he does use, thus saving the overhead of dummy subroutine calls 
 and extraneous preparation of data. This has the disadvantage that if 
 DFASUB is updated, all this work must be repeated on the new version; 
 whereas if the subroutine calls are left alone, those that work with 
 one version of DFASUB will probably work with any new version. 
 
 Usable versions of COPYZ and KNTSPI are also included in 
 Appendix I. If the computer in use has some way to count suppressed 
 underflow and overflow under program control, KNTSPI can be programmed 
 to return a value greater than 3 when too many errors occur. 
 
18 
 
 Finally, DFASUB must provide the initial values of Y(l,l), 
 Y(2,l), and YL(J). The general -purpose system calls a program DIFMF3 [8] 
 which, if given either Y(l,l) or Y(2,l) for each I and YL(j) , will find 
 all necessary initial values. If all initial values Y(l,l) are known , 
 then by setting Y(2,l) = for I = 1,NY, calling DIFFIM once and setting 
 
 Y(2,I) = -DY(I), 
 good approximations can be given to all necessary initial values. 
 
19 
 5. PROGRAM LOGIC 
 
 An anthropomorphic view of the actions of DIFSUB while solving 
 a system of equations follows. A main program specifically written for 
 the PDP -8/ GRAPHICS system [l] provides the calling sequence that enters the 
 subroutine DIFSUB and the subroutines which DIFSUB calls are compiled by 
 the system loaded with DIFSUB. A detailed description of what each of 
 these do is to he found in Section 3. 
 
 The program is entered with JSTART = 0. The bookkeeping variables 
 that record the number of steps, NS , and the number of matrix evaluations, NW, 
 are initialized to 0. Both y and y' are expected to be in array Y(l,l) and 
 Y(2,l). This can be accomplished by a call to the steady state problem 
 package DIFMF3 [8]. Y(2,l) is scaled by H, YMAX(l) is set to 1.0, the order 
 NQ is set to 1, and the variables l_. are placed in A(j) , J = l(l)NQ+l. 
 
 J 
 
 Various error test criteria are computed in E, EUP , and EDWN ; ENQ1 , ENQ2 , 
 and ENQ3 held l/(2*NQ), l/(2*(NQ+l), l/(2*(NQ+2) the power to which the 
 ratio of computed error to desired error must be taken to find the ratio 
 of new to old step size. 
 
 MATSET is called with IT = 1 to evaluate any constant elements of 
 PW; MATIN1 inverts these elements. 
 
 The value of Y(J,l), J = l(l)NQ+l, I = l(l)NY are saved in 
 SAVE(J,l) in case they are needed to restart the program; YL(l), I = l(l)NL 
 is saved in YLSV(l). H, T, and NQ are also saved. DIFSUB then calls S2 
 which evaluates any variables there are functions only of T and G. MATSET 
 is called with an argument of 2 to evaluate and update any part of the 
 Jacobian matrix PW which would change as a result of the call to S2. DIFSUB 
 then calls MATIN2 which does an inversion of parts of PW that change as a 
 result of calling S2 and MATSET ( . . . ,2 , . . . ) . 
 
20 
 
 DIFSUB multiplies Y(J,l) by the Pascal triangle matrix as 
 described in [6]. This replaces y. with y. , s, the predicted value 
 of the nonlinear variables. 
 
 The corrector step is a long loop that starts by initializing 
 ERROR to 0. DIFSUB calls DIFFUN (T, G,DY ,Y ,YL,HINV) which places in DY 
 the updated y! , v value. MATSET (...,3,...) is now called to update 
 the Jacobian as a result of having the predicted value of y. and the 
 corrected value of y! available. MATIN3(PW) is called to make any 
 changes in the LU decomposition as a result of the call to MATSET 
 (...,3,...). IWEVAL is set to -1 so that all of this re-evaluation of 
 PW need not be done again unless it seems necessary. MATMUL(PW,DY ,Fl) 
 is called which does the back substitution that computes the amount by 
 which Y(j,l) and YL(l) are to be changed to hold the corrected value. 
 YL(I),Y(1,I) and Y(2-,I) are corrected by A(J)*Fl(l ) ; ||Ay, J| is 
 then computed. If | |Ay| | is less than BND, the computed error bound 
 for the order, then the program continues. If not, MATMUL is again 
 called, YL(l), Y(l,l), Y(2,l) are corrected, Ay , \ is added to Ay/-, \ 
 and placed in ERR0R(*). 
 
 Actually, not only the present corrector term but also an 
 estimate of the next corrector term is compared to BND. Since the L2 
 norm of the i-th corrector term | | ERROR | | is approximately 
 R*| | ERROR I | for R constant throughout each small region of 
 integration, the predicted next value of the corrector can also be 
 compared to BND and the method is considered to have converged if 
 
 MIN(| | ERROR 1 | | , R*| |ERR0R 1+1 | | *2) <BND. 
 R is initially set to | | ERROR | \/\ | ERROR | | and updated at each step. 
 
21 
 
 Assume that the corrector term did converge in two corrections. 
 Then the estimated error 
 
 | |error(i)/ymax(i)| I 
 
 is computed; and if this is less than E, the maximum error hound allowed 
 to have the real error per step less than EPS, the program continues. 
 Otherwise, a better value of H for this or one lower order is computed 
 as described below, KFLAG is reduced by 2, the original values for this 
 step are returned, and the step is repeated for the new step size. Should 
 KFLAG reach -5, indicating three steps taken with three different step 
 sizes yet none gave an acceptably small error, the order is reduced to 1 
 and another attempt is made. If this is unsuccessful, KFLAG is set to 
 -k and DIFSUB makes an error return to its calling program. If at some time 
 H becomes less than HMIN, KFLAG is set to -1 and DIFSUB makes an error return 
 
 Otherwise, the rest of the Y(J,l) array, if NQ > 1, is updated 
 using the contents of ERROR(l) by Y(J,l) = Y(J,l) + ERROR( I ) *A( J) , then 
 DIFSUB returns to the beginning of the program. JSTART has been set to 
 1 so no initialization occurs. The integration step procedes as in the 
 first step, except that after NQ successful steps ERROR(l) is stored in 
 ERSV(l) and the next step (if successful) uses this and other information 
 
 to check for a better step size and order. 
 
 (k) 
 Since we also have an estimate for y and y /k! stored in the 
 
 n n 
 
 array Y(7,l), and since ERROR(l) stores the current step's accumulated 
 
 corrections for Y(l,l), the backward difference of the last component of 
 
 Y(K,N) is ERROR(l)*A(K) which is an estimate for 
 
 ,k+l (k+l) 
 h y 
 _n 
 
 k! 
 
22 
 
 Also, since ERROR(l) corresponding to the last y are saved in ERSV(l) , 
 
 we know that 
 
 ERSV(I) - ERROR(I) = h (k+2) y (k+2 } (t Q ) /k! . 
 These values can he used to compute the "best step size and order in 
 the following ways . 
 
 After the predictor-corrector method converges, we can check 
 to see if the truncation error was within hounds by seeing if 
 
 2 Ml C k+1 b ^ y <^> g 
 EPS 1 £ ' YMAX(I) > 
 
 Ml C. *ERROR(l)*K!*A(K) 
 
 = y f k+1 \^ 
 
 * v YMAX(I) y 
 
 by comparing 
 
 to 
 
 E = (EPS/C. ._ A(K)*K!) 2 
 
 Ml 
 _ v ^ , ERR0R(lK 2 
 
 D " z z 4max(i) j ' 
 
 When we want to see about changing order and step size (after 
 K+N*9, N = 0,1,..., successful steps at order K and step size H) , we can 
 compute the best step size to use at the current order, order 1 higher 
 and order 1 lower. Since 
 
 D ^ f H NEW s2 (p+1) 
 E ' [ H ' 
 
 by computing PR2 = l^D/E) 1 ' 2 ^"*" 1 ' , 
 
 MEW = H/PR2 
 
 will give the best step size for order K (1.2 is a heuristic constant 
 
 k+2 
 to compensate for ignoring the 0(h ) terms). 
 
23 
 
 letting 
 
 and 
 
 If one order lower method is used, then 
 
 C h k v' k ^ C *v(k l)*r>' 
 
 EPS ^ Z( YMAX(l) } =E( YMAX(l) } 
 
 EDWN = (EPS/C *k! ) 2 
 
 then 
 
 D " M YMAX(l) J 
 
 PR1 = 1.3(D/EDWN) 1/2p 
 
 gives the factor for the "best step size for order k-1 methods 
 Finally, if order one higher is used, we need 
 
 n >k (k) 
 2 C k+2 h y I s2 
 
 ^ ^ mx(i) } 
 
 C. Q (ERSV(I)-ERR0R(I))A(K)*K!) 
 v YMAX(I) 
 
 Letting EUP = (EPS/C *A(K)*p!) 2 and 
 
 k+2 
 
 then 
 
 and 
 
 _ . w D/ERR0R(lh 2 
 D " M YMAX(I) j 
 
 D _ ,HNEW,2 (p+2) 
 E " ^ H ' 
 
 pes = i.m^) 1/2 (p+2) 
 
 gives the best step size. The 1.3 and l.k terms bias the method toward 
 not changing the order or picking one order lower, respectively, since 
 these would minimize overhead calculations. 
 
24 
 
 Once a new step size has been computed, the array Y(J,l) must 
 be changed to reflect the new step size since Y(J,l) = Y h /(J-l)!. 
 
 This is accomplished "by multiplying the elements of SAVE ( J, I ) (if the 
 single step error were too large and the step is being repeated) or of 
 Y(J,l) (if a new step size is being prepared for the next step) by 
 (H/HOLD) (J_l) . 
 
 This is done in a DO loop after statement 750 or 800 depending 
 on whether SAVE or Y is supplying the values to be reached. If no 
 significant improvement in step size can be made, IDOUB is set to 9 (or 
 other large integer) and the step size is re-evaluated again if the 
 steps continue to be successful. 
 
 If, contrary to our assumption, the corrector steps do not 
 converge after two tries, then H is reduced to E/k unless IWEVAL = 
 indicating that PW should be re-evaluated first. Y(J,l) is scaled to 
 reflect this new step size and a new integration step is attempted. If 
 H cannot be reduced, KFLAG is set to -3 and DIFSUB makes an error return. 
 
 Each time a new step is entered, DIFSUB checks for T > TEND 
 and KFLAG < 0. If either of these is true, the program exits to the 
 calling routine with JSTART set to the present value of the order NQ , 
 HNEW set to the best H for the next step, H set to the present H, and 
 Y(j,l) set to its last successful value. 
 
 The entry point REDSUB is used to re-enter the program after 
 problems with the matrix processing routines. 
 
 This section, while describing the general program logic, leaves 
 out many details of programming that would be unwieldly to describe here. 
 The theory and some ideas about the programming techniques are in the other 
 sections, and after reading these and using the program copy in Appendix II, 
 this section should be sufficient to understand the program logic. 
 
25 
 
 LIST OF REFERENCES 
 
 [l] Gear, C. W. , "An Interactive Graphic Modeling System," Department 
 of Computer Science Report No. 318 , University of Illinois at 
 Urb ana-Champaign, 19&9- 
 
 [2] Gear, C. W. , NUMERICAL INITIAL VALUE PROBLEMS IN ORDINARY DIFFERENTIAL 
 EQUATIONS, Prentice -Hall: New Jersey, 1971. 
 
 [3] Nordsieck, A., "On the Numerical Integration of Ordinary Differential 
 Equations," Math. Comp. , l6, 1962, pp. 22-^9. 
 
 [h] Gear, C. W. and Watanabe, D. S., "Stability and Convergence of 
 
 Variable Order Multistep Methods," submitted to SIAM Journal 
 on Numerical Analysis. 
 
 [5] Gear, C. W. , "Simultaneous Numerical Solution of Differential- 
 Algebraic Equations," IEEE Trans. CT , 18, #1, 1971, pp. 89-95. 
 
 [6] Gear, C. W. , "DIFSUB for the Solution of Ordinary Differential 
 Equations," CACM , Ik, #3, 1972, pp. 185-190. 
 
 [7] Moler, C. B. , "Matrix Computations with FORTRAN and Paging," 
 CACM , 15, #U, 1972, pp. 268-270. 
 
 [8] van Melle, B. , "An Operating Manual for DIFMF3," Department of 
 Computer Science File No. 870 , University of Illinois at 
 Urban a- Champaign , 1972. 
 
26 
 
 APPENDIX I 
 SAMPLE PROGRAM 
 
27 
 
 c**** 
 
 c* 
 
 c**** 
 
 c**** 
 
 c* 
 
 c**** 
 
 10 
 
 20 
 
 c**** 
 
 c* 
 
 c**** 
 
 30 
 35 
 
 999 
 50 
 
 IMPL 
 DIME 
 
 1 G(l 
 DATA 
 
 «• -45 
 
 ♦ 52. 
 
 ♦ 452 
 ***** 
 
 SUPP 
 
 ***** 
 
 DATA 
 
 ♦ 5.D 
 CALL 
 CALL 
 
 ***** 
 SET 
 
 ***** 
 nr i 
 
 Yd, 
 Yd 
 
 Yd, 
 
 DO 2 
 Y(2, 
 
 YL( 1 
 YLI2 
 
 M = N 
 ***** 
 FIND 
 ***** 
 CALL 
 DO 3 
 Y(2, 
 JSTA 
 CALL 
 
 1 JST 
 
 2 YL 
 WRIT 
 FORM 
 STOP 
 FND 
 
 AL*8( A-H.Q-Z) 
 
 (7,6),YL(2),SAVE(7,8),YLSV(2),PW(R4), 
 
 , 0Y(8) , ERSV(6) .ERROR (6) ,F1( 8) ,EQN(1 t ,VAP (1) ,YMAX( R) 
 
 5002 5,-452.49975,-47.499 75,-52.5002 5, 
 
 ,44 7.50025,52.50025,4 7.49975,-47.49975, 
 
 47.5002 5,4 52.499 75,-52.5002 5,47.49975, 
 
 447.50025/ 
 
 ********************************* 
 
 ERFLOW AND UNDERFLOW 
 ********************************* 
 
 ,EPS,HMAX,HMlN,H,T,TENP/8,2,8,l.D-4, 
 0, l.D-4,0. rO, L. DO, 1.03/ 
 (208,256,-1,1 ) 
 (207,256,-1,1 ) 
 
 ********************************* 
 
 VALUES AND ZERO FIRST DERIVATIVES. 
 ********************************* 
 
 ICIT RE 
 NSIPN Y 
 6),T(2) 
 G/447. 
 2.49975 
 50025,4 
 .49975, 
 ******* 
 
 KESS OV 
 ******* 
 
 N,NL,M 
 2,1.D-1 
 
 ERPSET 
 
 EPRSET 
 
 *** **** 
 
 INITIAL 
 ******* 
 
 1=1,4 
 I)=-l. 
 ,5)=1. 
 6) = 1. 
 1=1,6 
 I »=0.D0 
 )=-2. 
 )=-3. 
 
 **************************************** 
 
 FIRST DERIVATIVES OF N1NLINEAP VARIABLES 
 **************************************** 
 
 DIFFUN(T,G,DY,Y, YL , L.I 
 1=1,6 
 I )=-PY( I ) 
 RT = 
 
 DFASUB(DY,FPS,EQN,ERPOR,FRSV,Fl ,G, H, HMAX, HMIN , 
 ART,KFLAG,M,6,6,N,NL,PW,SAVE,T,TEND,VAR,Y,YL, 
 SV, YMAX) 
 E(6,999)KFLAG 
 ATI ' KFLAG =• ,16) 
 
 SUBROUTINE S2IT.G) 
 REAL*8 T«2) 
 T(2)«DEXP(-T(1I I 
 RETURN 
 ENO 
 
 SUBROUTINE DIFFUNJ T, G ,DY, Y, YL , HI NV ) 
 
 IMPLICIT REAL*8 (A-H.Q-Z) 
 
 DIMENSION G(4,4),DYI8),Y(7,6),YL(2> ,T(2> 
 
 R=»(Y(l.l)*Y(l,2)*Y(l,2>*Y(l,4)l/2. 
 
 S = 0. 
 
 DO 20 !«lt* 
 20 S = S*(R-Y(l, I ) » **2/2- 
 
 DO 3 1=1,4 
 
 DY( I l=HlNV*Yt2,II-S*(R-Y( 1, I) )*»2 
 
 00 2 5 J =1,4 
 25 DY< I )=DY(I )*G( I , J > *Y ( 1, J) 
 30 CONTINUE 
 
 0Y«5)=HINV*(Y( 2,5)*Y( 1, ll*Y( 2,6)*Y(2,l)*Yl 1,61 I 
 
 DY(6)*2.*Y( l,6)*Yd,6)»*3-Y( I, ll*Ylt 1I-1.-TC2I 
 
 0Y<7)=YL( 1)-YL(2)*Y( l,l)*Y( 1,6) 
 
 0Y(8) = YLd )*YL (2)»5.*Y(1,1)*Y( 1,2) 
 
 RETURN 
 
 END 
 
28 
 
 10 
 
 SUBROUTINE CDPYZ(S,Y,L) 
 DIMENSION Y( 1) ,S ( 1 ) 
 DO 10 J = 1,L 
 S(J)=Y(J) 
 
 RETURN 
 END 
 
 FUNCTION KNTSPim 
 
 KNTSPI=I 
 
 RETURN 
 
 END 
 
 SUBROUTINE 
 
 RETURN 
 
 END 
 
 MATINKP) 
 
 SUBROUTINE MATIN2IPI 
 RETURN ' IK| 
 
 END 
 
 C* 
 
 C* 
 
 c* 
 
 c* 
 c* 
 c* 
 
 c* 
 
 c + 
 
 SUBROUTINE M AT SET ( A2 , CY , EPS , D4, G,HI NV, D5 ,06 , D7, Nt NY, I CHK , 
 1 PW,F1,T,D9,Y,YL) 
 
 IMPLICIT REAL*8( A-H.Q-Z) 
 
 DIMENSION OYU), G(l), PW 11), TU). Y(7,1),YL(1),F1(1) 
 
 SEE IF THIS IS A CALL TO BE FULLY HANDLED. 
 
 IFI ICHK.NE.3) GO TO 100 
 
 NL=N-NY 
 
 NN=N*N 
 
 DO 20 1=1, NN 
 20 PW( I ) = 0. 
 
 DO 40 J=1,NY 
 
 SAVE COLUMN ELEMENTS 
 
 F=Y( 1, Jl 
 
 E=YI2,J) 
 
 R=FPS*DMAX1 (EPS, DABS ( F) ,DABSIG) ) 
 
 FIND A DIFFERENCE ELEMENT = MAX( EPS**2 , EPS*Y( 1 , J ) , EPS*Y ( 2 , J ) ) 
 
 Yd, J) = Y(1, J|«-R 
 
 Y(2, J)=Y<2, J)-A2*R 
 
 CALL DIFFUN(T,G,F1,Y,YL,HINV) 
 DF A(2) DF 
 
 ASSIGN VALUES TO -- 
 
 DY H DY» 
 
 DO 30 I=1,N 
 30 PW< I*(J-l)*N)= IFK I)-DY( I»)/R 
 
 RESTORE COLUMN ELEMENTS 
 
 Y(2,J)=E 
 40 Y< 1, J |=F 
 
 IF ANY LINEAR ELEMENTS, FIND DF/DYL. 
 
 IF(NL.E0.0)GO TO 100 
 
 DO 80 J=1,NL 
 
 F=YL( J) 
 
 R*EPS*0MAX1IEPS,DABSI F) » 
 
 YL( J)=YL(J)*R 
 
 CALL DIFFUN(T,G,F1 ,Y,YL,HINV) 
 
 DO 70 1 = 1, N 
 70 PW( I*JJ*NY-1)*N)»IFH II-DYI I) l/R 
 80 YL<J)*F 
 100 RETURN 
 
 END 
 
 SUBROUTINE MATIN31PW) 
 
 IMPLICIT REAL*8 (A-H.Q-Z » 
 
 DIMENSION PW(1 ) 
 
 COMMON N 
 
 COMMON /IPP/ IP(40» 
 
 CALL DECOMP(N,N,PW,IP) 
 
 RETURN 
 
 END 
 
 in 
 
 SUBROUTINE M AT MUL( PW , DY , Fll 
 
 IMPLICIT REAL*8 (A-H.O-Z) 
 
 COMMON N 
 
 COMMON /IPP/ IPUO) 
 
 DIMENSION PWI1),DY(1) ,FIC 1) 
 
 DO 10 I-1,N 
 
 Fl( I )*DY(I) 
 
 CALL S0LVE(N,N,PW,F1, IP) 
 
 RETURN 
 
 END 
 
29 
 
 SUBROUTINE CECOMPI N , NCI M, A , IP ) 
 
 IMPLICIT REAL*8<R-H,Q-Z) 
 C 
 
 C MATRIX TR IANGULARIZATICN BY GAUSSIAN ELIMINATION. 
 C 
 
 C INPUT... 
 
 C N = ORDFP OF MATRIX 
 
 C NDIM = DECLARER DIMENSION OF ARRAY A. 
 
 C A = MATRIX TO RF TR I ANGULAR I ZED. (FOR STIFF METHODS, A IS SINGLE 
 
 C PRECISION; ALL OTHER VARIABLE AP F OOUBLE PRECISION.) 
 
 C OUTPUT... 
 
 C A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U. 
 
 C All, J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I-L. 
 
 C IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW. 
 
 C IPIN) = «-l)**INUMBEP OF INTERCHANGES) OR 0. 
 
 C USF 'SOLVE' TO OBTAIN SCLUTION OF LINEAR SYSTEM. 
 C DPTERM(A) = IP(N)*A(1,1)*A(2,2)*...*AIN,N). 
 C IF IPIN) = 0, A IS SINGULAR, SOLVE WILL DIVIDE BY ZERO. 
 
 C 
 
 DIMENSION AINDIM.NDIM), IP(NDIM) 
 
 DIMENSION BI4,*) 
 
 COMMON NFNS,NW,B 
 
 NW = NW + 1 
 
 IPIN) = 1 
 
 DO 6 K=1,N 
 
 IFIK.FQ.N) GO TO 5 
 
 KP1 * K+l 
 
 M = K 
 
 DO 1 I=KP1,N 
 
 IF(ABS(A(I,K)).GT.ABS(A(M,K)I) MM 
 
 1 CONTINUE 
 IP(K) = M 
 
 IFIM.NE.K) IPIN) = -IPIN) 
 
 T = A(M.K) 
 
 A(M,K) = AIK.KI 
 
 AIK.K) * T 
 
 IFIT.EQ.O) GO TO 5 
 
 DO 2 I*KP1,N 
 
 2 A( I,K) = -All ,K) / T 
 
 ?° \,i ! n 1,N SUBROUTINE SOL VE IN.ND IM , A, B, IP ) 
 
 AIM,;! ^AIK.J, c ,MPUCIT REAL * 8 'B-H.O-II 
 
 Jf^EoIoI) GO TO 4 I S ° LUTI0N ° F L,N6AB SYSTEM ' A * X ' B « 
 
 DO 3 I-KPUN c , NPUT ... 
 
 I rn^JM,^' " AII ' JI * ACI ' K, * T C N - ORDER OF MATRIX. 
 
 * CONTINUE c N0 , M , DECLARED DIMENSION OF ARRAY A. 
 
 5 IF(A(K,K).EQ.O.) IPIN) -0 C A « TR I ANGUL ARI ZED MATRIX OBTAINED FROM 'DECOMP' 
 
 6 CONTINUE C B ■ RIGHT HAND SIDE VECTOR. 
 
 RETURN C IP - PIVOT VECTOR OBTAINED FROM 'DECOMP'. 
 
 END C OUTPUT... 
 
 C B ■ SOLUTION VECTOR, X. 
 
 DIMENSION AINDIM.NDIM), BINDIMI, IP(NDIM) 
 IFIN.EO.l) GO TO 9 
 NM1 - N-l 
 DO 7 K"1,NM1 
 KP1 ■ K ♦ I 
 M ■ IPIK) 
 T-BIM) 
 BIM) - P(K) 
 BIK) - T 
 DO 7 I-KP1.N 
 
 7 BID - B(I) ♦ AII,K)*T 
 DO 8 KB-l.NMI 
 
 KM1 - N - KB 
 
 K ■ KMl ♦ I 
 
 BIK) ■ BIKJ/AIK.K) 
 
 T * -BIK) 
 
 DO 8 I-1.KM1 
 
 8 Bl I ) « BID ♦ A(I,K)*T 
 
 9 BID - BU)/AI1, 1) 
 RETURN 
 
 END 
 
30 
 
 APPENDIX II 
 DFASUB LISTING 
 
31 
 
 SUBROUTINE OFASUB ( DY , EPS , EON, ERP OR , EP S V, 
 + Fl ,G,H,HMAX,HMIN, 
 
 ♦ JSTAPT,KFLAG,M,MAXDER,M1 ,N, 
 «■ NL, PW, SAVE, T, TEND, VAR,Y, 
 «■ YL.YLSV, YMAX) 
 
 IMPLICIT PEAL** (A-H.O-Z) 
 
 C* 
 
 (-* THE PARAMETERS to THE SUBROUTINE PIFSUB HAVE 
 
 C* THE FOLLOWING MEANINGS: 
 
 C* 
 
 C* N THE NUMBER TF VARIABLFS. 
 
 C* NL THE NUMBER OF LINEAR VARIABLES 
 
 C* (NY = N-NL IS THE NUMBER OF VARIABLES WITH DERIVATIVES) 
 
 C* Ml THE NUMBER OF EQUATIONS TO TAKE PART IN THE ERROR TEST, 
 
 C* TFND END CRITERION 
 
 C* T THE INDEPENDENT VARIABLF. 
 
 C* G AN ARRAY OF GLOBAL VARIARLES 
 
 C* Y A 7 BY NY ARRAY CONTAINING THE DEPENDENT VARIABLFS 
 
 C* AND THEIR SCALED DERIVATIVES. Y(J*l,I) CONTAINS 
 
 C* THE J-TH DERIVATIVE OF Yd) SCALFD BY 
 
 C* H**J/FACTPRI AL( J) , H THE CURRFNT STEP SIZE. 
 
 C* ON THE FIRST ENTRY, THE CALLER SUPPLIES 
 
 C* Ylltl) ANT Y(2,I), UNSCALEO. (IF THE CALL TO 
 
 C* DIFSUB WAS PRECEDED BY A CALL TO PIFMF3, THE 
 
 C* USER NEED NOT TOUf"H Y AT ALL). THE PROGRAM 
 
 C* WILL SCALE Y(2,I) BY H. ON ANY SUBSEQUENT 
 
 C* ENTRY, THE PROGRAM ASSUMES THAT THE Y VALUES 
 
 C* HAVE NOT BFEN CHANGED SINCE THE LAST EXIT 
 
 C* FROM DIFSUB, AND WILL SCALE THESE VALUES IF 
 
 C* THE CALLER HAS CHANGED H. 
 
 C* IF IT IS DESIRED TO INTERPOLATE TO NON-MESH 
 
 C* POINTS THESE VALUES CAN BE USED. IF THE CURRENT 
 
 C* STEP SIZE IS H AND THE VALUE AT T*E IS NEEDED, 
 
 C* FORM S = E/H AND THE COMPUTE 
 
 C* NO 
 
 C* YCIMT+E) = SUM Y(J*1,I )*S**J 
 
 C* J*0 
 
 C* YL AN ARRAY OF NL VARIABLES THAT APPEAR ONLY LINEARLY. 
 
 C* CALLER MUST SUPPLY VALUES FOR THESE VARIABLES. 
 
 C* SAVE AN ARRAY OF LENGTH AT LEAST 7*N. 
 
 C* H THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP. 
 
 C* H MAY BE AOJUSTED UP OR DOWN BY THE PROGRAM 
 
 C* IN ORDER TO ACHEIVE AN ECONOMICAL INTEGRATION. 
 
 C* HOWEVER, IF THE H PROVIDED BY THE USER DOES 
 
 C* NOT CAUSE A LARGER ERROR THAN REQUESTED, IT 
 
 C* WILL BE USED. TO SAVE COMPUTER TIME, THE USER IS 
 
 C* ADVISED TO USE A FAIRLY SMALL STEP FOR THE FIRST 
 
 C* CALL. IT WILL BE AUTOMATICALLY INCREASED LATFR. 
 
 C* HMIN THE MINIMUM STEP SIZE THAT WILL BE USED FOR THE 
 
 C* INTEGRATION. NOTE THAT ON STARTING THIS MUST BE 
 
 C* MUCH SMALLER THAN THE AVERAGE H EXPECTED SINCE 
 
 C* A FIRST ORDER METHOD IS USED INITIALLY. 
 
 C* HMAX THE MAXIMUM ALL3WABLE STEP SIZE 
 
 C* EPS THE ERROR TEST CONSTANT. SINGLE STEP ERROR ESTIMATES 
 
 C* DIVIDED BY YMAX(I) MUST BE LESS THAN THIS 
 
 C* IN THE EUCLIDEAN NORM. THE STEP AND/OR ORDER IS 
 
 C* ADJUSTED TO ACHEIVE THIS. 
 
 C* YMAX A VECTOR OF LENGTH NY WHICH CONTAINS THE MAXIMUM 
 
 C* OF EACH Y SEEN SO FAR. ON THE FIRST CALL, THESE 
 
 C* WILL BE INITIALIZED AS YMAX(I) - MAX( 1 , | Y ( 1 , I ) | ) 
 
 C* ERROR A VECTOR OF LENGTH NY. 
 
 C* KFLAG A COMPLETION CODE WITH THE FOLLOWING MEANINGS: 
 
 C* *l THE INTEGRATION WAS SUCCESSFUL. 
 
 C* -1 THE ERROR TEST FAILED FOR H > HMIN. 
 
 C* -3 THE CORRECTOR FAILED TO CONVERGE FOR 
 
 C* H > HMIN. 
 
 C* -2 TOO MANY FLOATING-POINT EXCEPTIONS 
 
 C* OCCURRED DURING LAST STEP. 
 
 C* -4 THE CORRECTOR FAILED THREE TIMES WITH 
 
 C* EVEN THE FIRST-ORDER METHOD. 
 
 OF AS 
 
 001 
 
 DFAS 
 
 00? 
 
 DFAS 
 
 003 
 
 DFAS 
 
 004 
 
 OFAS 
 
 005 
 
 DFAS 
 
 OOh 
 
 DFAS 
 
 007 
 
 DFAS 
 
 008 
 
 DFAS 
 
 009 
 
 DFAS 
 
 010 
 
 DFAS 
 
 311 
 
 DFAS 
 
 012 
 
 DFAS 
 
 013 
 
 DFAS 
 
 014 
 
 DFAS 
 
 015 
 
 P C AS 
 
 016 
 
 DFAS 
 
 017 
 
 DFAS 
 
 010 
 
 DFAS 
 
 019 
 
 DFAS 
 
 020 
 
 DFAS 
 
 021 
 
 DFAS 
 
 022 
 
 DFAS 
 
 023 
 
 DFAS 
 
 024 
 
 DFAS 
 
 025 
 
 DFAS 
 
 026 
 
 DFAS 
 
 027 
 
 DFAS 
 
 02 8 
 
 OFAS 
 
 029 
 
 DFAS 
 
 030 
 
 DFAS 
 
 031 
 
 DFAS 
 
 032 
 
 DFAS 
 
 033 
 
 DFAS 
 
 034 
 
 DFAS 
 
 035 
 
 DFAS 
 
 036 
 
 DFAS 
 
 037 
 
 DFAS 
 
 038 
 
 DFAS 
 
 039 
 
 DFAS 
 
 040 
 
 DFAS 
 
 041 
 
 DFAS 
 
 042 
 
 DFAS 
 
 043 
 
 DFAS 
 
 044 
 
 DFAS 
 
 045 
 
 DFAS 
 
 046 
 
 DFAS 
 
 04 7 
 
 DFAS 
 
 048 
 
 DFAS 
 
 049 
 
 DFAS 
 
 050 
 
 DFAS 
 
 051 
 
 DFAS 
 
 052 
 
 DFAS 
 
 053 
 
 DFAS 
 
 054 
 
 DFAS 
 
 055 
 
 DFAS 
 
 056 
 
 DFAS 
 
 057 
 
 DFAS 
 
 05 8 
 
 DFAS 
 
 059 
 
 DFAS 
 
 060 
 
 DFAS 
 
 061 
 
 DFAS 
 
 062 
 
 DFAS 
 
 063 
 
 DFAS 
 
 064 
 
 DFAS 
 
 065 
 
 DFAS 
 
 066 
 
 DFAS 
 
 067 
 
 DFAS 
 
 068 
 
 DFAS 
 
 069 
 
 DFAS 
 
 070 
 
 DFAS 
 
 071 
 
32 
 
 c* 
 
 r * 
 C* 
 
 r* 
 C* 
 r * 
 
 C* 
 C* 
 C* 
 C* 
 C* 
 
 c* 
 
 r* 
 C* 
 C* 
 c* 
 r* 
 c* 
 c* 
 c* 
 c* 
 
 JSTAPT AN' INPU T INCICATOR WITH T 
 
 PERFORM THF FIRS 
 
 INITIALIZES IT 
 AT IVES IN Y<2, 
 THF INTEGRATIO 
 ANY SUBSEQUENT 
 WITH JSTART = 
 
 1 CONTINUE FROM TH 
 
 UNTIL T > TFN'D 
 JSTART IS SET TO NO, TH 
 THE METHOD, AT EXIT. 
 
 MAXPFR THE MAXIMUM DERIVATIVF TH 
 METHOP. IT MUST NOT EXC 
 
 PW A VFCTOR O c LENGTH N**2*-20 ( 
 GENERATED BY MATSFT AND 
 
 YLSV A VFCTOR OF LENGTH NL . 
 
 DY A VECTOR "IF LENGTH M, OUT 
 
 ERSV A VECTOR OF LENGTH NY. 
 
 Fl A VECTOR OF LENGTH N, OUT 
 
 EQN.VAR VECTORS USED BY MATSET. 
 
 HE FOLLOWING MEANINGS: 
 T STEP. THF ROUTINE 
 SELF, SCALES THF OER IV- 
 I ) AND THEN PERFORMS 
 N UNTIL T > TFND. 
 
 CALLS SHOULD BE MADE 
 I. 
 E LAST STEP, INTEGRATING 
 
 F CURRENT ORDER OF 
 
 AT SHOULD BE USED IN THE 
 EED 6. 
 REAL*4) , 
 USFD BY MATINV.MATMUL. 
 
 PUT OF DIFFUN. 
 
 PUT OF MATMUL. 
 
 C**** 
 
 C* 
 
 C**** 
 
 c 
 
 c* 
 
 ***************************************************************** 
 DIMENSION T(l),G(ll,Y(7,l),YL(l),SAVE(7,l),YMAX(l) 
 DIMENSION ERROR ( 1 ) , PW ( 1 ) , YLSV( 1 ) ,DY ( 1 1 , ERSV ( 1 ) 
 DIMENSION Fill I , EON ( 1),VAP(1),A(7),PEPTST(6,3) 
 DATA PERTST /4 . ,9. , 16 . 0, 25 .0 , 36 . , 49 . , 9. , 16 . 0, 
 
 + 25.0,36.0,49.0,64.0,1.0,1.0,0.25, 
 
 ♦ 2.7889E-2,1.70 569E-3,6.83929E-5/ 
 DATA MF /?./, KZILCH / 0/ 
 
 ********************************************* 
 
 THIS COMMUNICATES THE VALUE OF N TO MATMUL AND MATIN3. 
 
 a**********************.**********.********** 
 
 COMMON N99 
 
 FORMAT </I5,!4,Dll.2,Dl0.2,6Dl5.5/(31X,6D15.5)) 
 
 FORMAT (31X.6D15.5) 
 
 FORMAT ('IMF =«,I2,« N =»iI3i' NL = , ,I3,' FPS =', 
 
 ♦ D9.2,' TFND =',D9.2,' H =',D9.2//I 
 FORMAT (• NS NW H» , 8X , • T( 1 ) • , 8X , 
 
 ♦ «Y< It*) AND YL(*)« //) 
 
 NY = N-NL 
 
 LENGTH IN WORDS OF Y, YL (FOR COPY) : 
 
 LCCPYY = NY*14 
 
 LCOPYL = NL*2 
 
 N99=N 
 
 IF (MAXDEP .GT. 61 MAXDER = 6 
 
 HINV=2.5 
 
 CALL DIFFUN ( T , G, F I , Y , YL ,HI NV ) 
 
 CALL MATSET ( 0, DY , FPS ,EON ,G ,HI NV , 
 
 ♦ PW.F1 ,T, VAR,Y,YL) 
 
 0, M,MF,N,NY, 1, 
 
 
 CALL 
 
 
 IF (JST 
 
 
 WRITE « 
 
 
 WRITE ( 
 
 r *********** * 
 
 C* 
 
 ON THE 
 
 c* 
 
 THE VFC 
 
 r ** ** * ....... 
 
 
 NS = 
 
 
 NW = 
 
 
 DO 20 J 
 
 
 YMAXI 
 
 20 
 
 Y(2,J 
 
 
 NO * 1 
 
 
 BR = 1. 
 
 
 ASSIGN 
 
 
 GO TO ? 
 
 r ************ 
 
 C* 
 
 SET THF 
 
 c* 
 
 TYPF. 1 
 
 c* 
 
 FUP IS 
 
 c* 
 
 DECRFAS 
 
 r ************ 
 
 MATIN1 (PW) 
 APT .GT. 01 GO TO 60 
 6,3) MF.N.NL, EPS, TEND, H 
 6,4) 
 
 ************************************** 
 FIRST CALL SET THE ORDER TO ONE. INITIALIZE 
 TOP YMAX, AND SCALE THE FIRST DERIVATIVES OF Y BY H. 
 ************************************** 
 
 = 1,NY 
 J) = DMAXK l.ODO.DABSIYI l.JI ) ) 
 ) = Y(2,J)*H 
 
 
 
 100 TO IRET 
 
 21 
 
 ************************************** 
 
 COFPFICIENTS THAT DETFRMINF THE ORDER AND THE METHOD 
 EISA TEST FOP FRRORS OF THE CURRENT ORDER NO. 
 TO TEST FOR INCREASING THE ORDER, EDWN FOR 
 INC. THE ORDER. 
 ************************************** 
 
 ^FAS 
 
 72 
 
 HFftS 
 
 073 
 
 DFAS 
 
 074 
 
 OCAS 
 
 075 
 
 DPAS 
 
 76 
 
 PFAS 
 
 077 
 
 T=AS 
 
 07R 
 
 HFAS 
 
 079 
 
 PFAS 
 
 OBO 
 
 PPAS 
 
 0R1 
 
 DFAS 
 
 0B2 
 
 P^AS 
 
 0R3 
 
 DFAS 
 
 0R4 
 
 DFAS 
 
 085 
 
 nFAS 
 
 0B6 
 
 DFAS 
 
 087 
 
 DFAS 
 
 088 
 
 DFAS 
 
 089 
 
 DFAS 
 
 090 
 
 DFAS 
 
 091 
 
 DFAS 
 
 092 
 
 **DFAS 
 
 093 
 
 DFAS 
 
 094 
 
 DFAS 
 
 095 
 
 DFAS 
 
 096 
 
 DFAS 
 
 097 
 
 DFAS 
 
 098 
 
 DFAS 
 
 099 
 
 DFAS 
 
 100 
 
 DFAS 
 
 101 
 
 DFAS 
 
 102 
 
 DFAS 
 
 103 
 
 DFAS 
 
 104 
 
 DFAS 
 
 105 
 
 DFAS 
 
 106 
 
 DFAS 
 
 107 
 
 DFAS 
 
 108 
 
 DFAS 
 
 109 
 
 DFAS 
 
 110 
 
 DFAS 
 
 ill 
 
 DFAS 
 
 112 
 
 DFAS 
 
 113 
 
 DFAS 
 
 114 
 
 DFAS 
 
 115 
 
 DFAS 
 
 116 
 
 DFAS 
 
 117 
 
 DFAS 
 
 118 
 
 DFAS 
 
 119 
 
 DFAS 
 
 120 
 
 DFAS 
 
 121 
 
 DFAS 
 
 122 
 
 DFAS 
 
 123 
 
 DFAS 
 
 124 
 
 DFAS 
 
 125 
 
 DFAS 
 
 126 
 
 DFAS 
 
 127 
 
 DFAS 
 
 128 
 
 DFAS 
 
 129 
 
 DFAS 
 
 130 
 
 DFAS 
 
 131 
 
 DFAS 
 
 132 
 
 DFAS 
 
 133 
 
 DFAS 
 
 134 
 
 DFAS 
 
 135 
 
 DFAS 
 
 136 
 
 DFAS 
 
 137 
 
 DFAS 
 
 138 
 
 DFAS 
 
 139 
 
 DFAS 
 
 140 
 
 DFAS 
 
 141 
 
 DFAS 
 
 142 
 
 DFAS 
 
 143 
 
 DFAS 
 
 144 
 
 DFAS 
 
 145 
 
33 
 
 170 CONTINUE 
 221 
 
 226 A(2) = -2.45D0 
 
 0<=AS 146 
 
 GOTO (221, 222, 223, 2 24,225,226), NO D p4 S 1*7 
 
 A(2J = -l.ODO DF? - S 1*« 
 
 r,o to 2?o ° c ^ ^ 9 
 
 222 A(2) = -1. 500 j>"S |5° 
 
 A(3) = -C.5D0 ° FA ^ 1S1 
 
 GP TO 230 nPiS 152 
 
 22* A(2) = -1.83333333333333333 - Ffi<? l 5 ^ 
 
 A(3) = -l.ODO nPAS 154 
 
 A<4) = -0.1666666666666666"' nFAS 155 
 
 GO TO 230 nCAS l5 * 
 
 224 A(2) = -2.08333333333333233 DFAS 157 
 A<3) = -1.45833333333333333 DF AS 1 5B 
 A<4) = -0.4166666666666667 DFAS 159 
 A(5> = -0.04166666666666667 DFAS 160 
 GO Tn 230 nFAS l61 
 
 225 A(2) = -2.8333333333333333 °FAS l62 
 A(3) = -1.875D0 DFAS I 63 
 A(4) = -0.7083333333333333 °FAS I 6 * 
 A(5) = -0.125D0 DFA< > 165 
 A(6) = -0.0083333333333333333 DFAS I 66 
 GO TO 230 DFAS 167 
 
 DFAS 168 
 
 At3) = -2.2555555555555555 nFA 5 169 
 
 A(4) = -1.2083333333333333 DFAS I 70 
 
 A(5) = -0.2430555555555555 OFAS 171 
 
 A(6) = -0.02916666666666667 DFAS 172 
 
 A(7) = -0.0013888888888888889 DFAS 173 
 
 230 K = NO+1 DFAS 17 * 
 
 ID0U8 = NO nFAS I 75 
 
 FN01 = 0.5/FLDATINQ) DFAS l76 
 
 ENQ2 = 0.5/FLDAT(K) °FAS I 77 
 
 FN03 = 0.5/FL3AT(NQ ♦ 21 DFAS I 78 
 
 PFPSH = EPS**2 DFAS 17q 
 
 E = PFRTSTINQ, 1) *PEPSH DFAS 180 
 
 EUP = PERTST(NQ,2J*PEPSH DFAS 181 
 
 FDWN= PERTST(N0,3I«PEPSH DFAS 182 
 
 BNO = ( EPS*FNQ3 )**2 0FAS I 83 
 
 IWEVAL = 1 OFAS 184 
 
 GO TO TRET, (100,250,630,751) OFAS 185 
 
 £***********************•*********•»********•****** DFAS 186 
 
 C* IF THE CALLER HAS CHANGED H, SCALE THE Y VARIABLES. DFAS 187 
 
 C* HNEW IS THF STEP SIZE USED ON THE LAST CALL. DFAS 188 
 
 £•***********•**********•***********••************* DFAS 189 
 
 60 IF (H .EO. HNEW) GO TO 100 OFAS 190 
 
 R = H/HNEW DFAS 191 
 
 ASSIGN 100 TO IRET OFAS 192 
 
 GO TO 800 DFAS 193 
 
 C************************************************** DFAS 194 
 
 C* BEFORE EXIT, JSTART IS SET TO THE ORDER OF THE METHOD, DFAS 195 
 
 C* AND THE CURRENT VALUE OF H IS SAVFD. DFAS 196 
 
 £***************•************•***•***************•* DFAS 197 
 
 70 KFLAG = -2 DFAS 198 
 
 80 JSTART = NO DFAS 199 
 
 HNEW = H DFAS 200 
 
 RETURN DFAS 201 
 
 r I************************************************ .. DFAS 202 
 
 C* PRINT DATA RELEVANT TO THE STEP, AND TAKE ANOTHER STEP DFAS 203 
 
 C* IF T < TEND. DFAS 204 
 
 £***«*****•****«************•*»********•********•** DFAS 205 
 
 90 NS = NS*1 DFAS 206 
 
 WRITE (6,11 NS,NW,H,T (1) ,(Y( 1, I ), 1=1, NY) DFAS 207 
 
 IF (NL .GT. 0) WRITE (6,2) < YL I I ) , I = 1 , NL ) OFAS 208 
 
 CALL ANSWER(Y,F1,T) DFAS 209 
 
 IF CKNTSPKO) .GT. 2) GO TO 70 DFAS 210 
 
 IF (KFLAG .LT. 0) GO TO 80 DFAS 211 
 
 IF (T(l) .GE. TEND) GC TO 80 DFAS 212 
 
 JSTART = 1 DFAS 213 
 
 C ********* ***************************************** DFAS 214 
 
 C* BEGIN BY SAVING INFORMATION FOR POSSIBLE RESTARTS. DFAS 215 
 
 e*****************************+******************** DFAS 216 
 
 100 CALL COPYZ (SAVE,Y,LCOPYY) DFAS 217 
 
 CALL COPYZ (YLSV,YL,LCOPYL) DFAS 218 
 
3k 
 
 RACUM = 1.0 DFAS 219 
 
 KFLAG = 1 DFAS 220 
 
 HOLD = H DFAS 221 
 
 NOOLD = MO DFAS 22? 
 
 TOLD = T(l) ^FAS 223 
 
 K7ILCH = 1 DFAS 224 
 
 ;*************************************** ********************************DFAS 225 
 
 C* THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY *DPAS 226 
 
 C* MULTIPLYING THE SAVED INFORMATION BY THE PASCAL TRIANGLE *DFAS 227 
 
 C* MATRIX. *DFAS 228 
 
 C* *********************************************** ***********************OFAS 229 
 
 250 T(l) = TdH-H DFAS 230 
 
 CALL S2(T,G) DFAS 231 
 
 CALL MATSET ( , DY , EPS , EON ,G,HI NV , 0, M,MF , N, N Y ,2 , DFAS 232 
 
 + PW,F1,T, VAR,Y,YL J DFAS 233 
 
 CALL MATIN2 <PW) DFAS 234 
 
 HINV = l.DO/H DFAS 235 
 
 DP ?60 J = 2fK OFAS 236 
 
 J3 = K*J-1 DFAS 237 
 
 DO 260 Jl = J f K DFAS 238 
 
 J2 = J3-J1 DFAS 239 
 
 DO 260 I = 1,NY DFAS 240 
 
 260 Y(J2,II = Y(J2,I) ♦ YIJ2*l,I) P p A'S 241 
 
 C******************* ********************************** ******************pFAS 242 
 
 C* UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. CONVFRGENCF IS TFSTED *DFAS 243 
 
 C* BY REQUIRING CHANGES TO BE LESS THAN BND WHICH IS DEPENDENT ON *DFAS 244 
 
 C* THE ERROR TEST CONSTANT. *DFAS 245 
 
 C* THE SUM OF THE CCPRECTIONS IS ACCUMULATED IN THE ARRAY *DFAS 246 
 
 C* ERROR(I). IT IS EQUAL TO THE K-TH DERIVATIVE OF Y MULTIPLIED *DFAS 247 
 
 C* BY H**K/(FACTORI AL (K-1)*A(K) ), AND IS THERE = ORE PROPORTIONAL *DFAS 248 
 
 C* TO THE ACTUAL ERRORS TO THE LOWEST POWER OF H PRESENT. (H«*K) *DFAS 249 
 £**** ******** ******* **************************************************** op as 250 
 
 DO 270 I = 1,NY DFAS 251 
 
 270 ERROR(I) = 0.0 DFAS 252 
 
 DO 430 L = 1,3 DFAS 253 
 
 CALL DIFFUN d , G, DY ,Y , YL, HI NV I DFAS 254 
 
 C ***** ******** ********************************************************** (if AS 255 
 
 C* IF THERE HAS BEEN A CHANGE OF ORDER OR THERE HAS BEFN TROUBLE *OFAS 256 
 
 C* WITH CONVERGENCE, PW IS RE-EVALUATED PRIOR TO STARTING THE DFAS 257 
 
 C* CORRECTOR ITERATION IN THE CASE OF STIFF METHODS. IWEVAL IS *DFAS 258 
 
 C* THEN SET TO -1 AS AN INDICATOR THAT IT HAS BEEN DONE. *DFAS 259 
 
 C ********* ******** *********** ******* ************************** **********qfaS 260 
 
 IF (IWEVAL .LT. 1 .AND. L .GT. I) GO TO 280 DFAS 261 
 
 CALL MATSET ( A ( 2) , CY , EPS , FQN, G,H I NV , 0, M , MF , N.NY , 3 , DFAS 262 
 
 ♦ PW,F1,T,VAR,Y,YL> DFAS 263 
 
 CALL MATIN3 IPW) DFAS 264 
 
 KFLAG * 1 DFAS 265 
 
 IWEVAL « -1 DFAS 266 
 
 NW ■ NW*1 DFAS 267 
 
 280 CALL MATMUL «PW,DY,F1) OFAS 268 
 
 IF (NL .LE. 0) GO TO 300 DFAS 269 
 
 DO 290 I - l.NL DFAS 270 
 
 290 YL(I) - YL(II - FlU+NY) DFAS 271 
 
 300 CONTINUE DFAS 272 
 
 DEL - 0.000 DFAS 273 
 
 DO 420 I - I, NY DFAS 274 
 
 Yd, I) » Yd, I) - F1III OFAS 275 
 
 Y(2,I» - Y(2,I» ♦ A<2»*F1(I) DFAS 276 
 
 ERROR(I) - ERROR(I) ♦ F1(I» DFAS 277 
 
 DEL - DEL ♦ (Fill )/YMAX< I) )**2 DFAS 278 
 
 420 CONTINUE DFAS 279 
 
 IFIL.GE.2) BR - DM AX 1 1 .9*BR , DEL/DEL 1 I DFAS 280 
 
 DELI ■ DEL DFAS 281 
 
 IF(DMIN1IDEL,BR*0EL*2.0).LE.BND| GO TO 490 DFAS 282 
 
 430 CONTINUE OFAS 283 
 
 C************* ******************************** **************************DFAS 284 
 
 C* THE CORRECTOR ITERATION FAILEO TO CONVERGE IN 3 TRIES. VARIOUS *DFAS 285 
 
 C* POSSIBILITIES ARE CHFCKED FOR. IF H IS ALREADY HMIN AND *DFAS 286 
 
 C* THIS IS EITHER ADAMS METHOD OR THE STIFF METHOD IN WHICH THE *0FAS 287 
 
 C* MATRIX PW HAS ALREADY BEEN RE-EVALUATED, A NO CONVERGENCE EXIT *DFAS 288 
 
 C* IS TAKEN. OTHERWISE THE MATRIX PW IS RE-EVALUATED AND/OR THE *DFAS 289 
 
 C* STEP IS RFDUCED TO TRY AND GET CONVERGENCE. *DFAS 290 
 
 C************* **************************** ******************************[)FAS 291 
 
35 
 
 440 T(l) = TOLD nc 4S 292 
 
 IF (IWEVAL) 445,455,450 OFAS 293 
 
 445 IF (H .LF. HMI N* 1 . 000000 1 ) GO TO 460 OFAS 2°4 
 
 450 RACUM = PACUM*0.25 OFAS 295 
 
 455 CONTINUE ^ FA S 29f 
 
 GO TO 750 OFAS 277 
 
 460 KFLAG = -3 OFtS 29R 
 
 470 CALL COPYZ ( Y, SAVF , LC CPYY ) ~>FAS 2 C 9 
 
 CALL COPYZ (YL.YLSV.LCOPYL) OFAS 300 
 
 H = HOLD DFAS 301 
 
 NO = NQOLD DFAS 302 
 
 GO TO 90 nF:A S 303 
 
 r *************************************** *********** OFAS 304 
 
 C* THF CORRECTOR CONVFRGEC ANO NOW THF ERROR TEST IS MADE. OFAS 305 
 
 C ************************************************** OFAS 306 
 
 490 P = 0.0 OFAS 307 
 
 00 500 I = 1,M1 OFAS 308 
 
 500 0=0* (ERROR (I )/YMAX< I ) )**2 OFAS 309 
 
 IWEVAL = OFAS 310 
 
 IF (D .GT. E) GO TO 540 OFAS 311 
 
 C ************************************************** D C AS 312 
 
 C* THE ERROR TEST IS OKAY, SO THE STEP IS ACCEPTED. IF IDOUR OFAS 313 
 
 C* NOW BECOMES NEGATIVE, A TEST IS MADE TO SEE IF THE STFP OFAS 314 
 
 C* SIZE CAN BE INCREASED AT THIS ORDFP OP ONE HIGHER OR OFAS 315 
 
 C* LOWER. THE CHANGE IS MADE ONLY IF THE STFP CAN BE IN- OFAS 316 
 
 C* CREASED RY AT LEAST 10*. IDOUB IS SET TO NO TO PREVENT DFAS 317 
 
 C* FURTHER TFSTING FOR A WHILE. IF NO CHANGE IS MADE, IDOUR DFAS 318 
 
 C* IS SFT TO 9. OFAS 319 
 
 C*************************************+************ DFAS 320 
 
 IF (K .LT. 3) GO TO 5 20 DFAS 321 
 
 DO 510 J = 3,K DFAS 322 
 
 DO 510 I = 1,NY DFAS 323 
 
 510 Y(J,I) = Y(J,I) ♦ A(J)*ERRnR(I) DFAS 324 
 
 520 KFLAG = >1 DFAS 325 
 
 IDOUB = IDOUB-1 DFAS 326 
 
 IF (IDOUB) 550,525,700 DFAS 327 
 
 525 DO 530 I = 1.M1 OFAS 328 
 
 530 ERSVdl = ERROR(I) DFAS 329 
 
 GO TO 700 OFAS 330 
 
 C************************************************** DFAS 331 
 
 C* THE ERROR TEST FAILED. IF JSTART = 0, THE DERIVATIVES DFAS 332 
 
 C* IN THE SAVE ARRAY ARE UPDATED. TESTS ARE THEN MADE TO DFAS 333 
 
 C* FIX THE STEP SIZE AND PERHAPS REDUCE THE ORDER. AFTER OFAS 334 
 
 C* RESTORING AND SCALING THE Y VARIABLES, THE STEP IS OFAS 335 
 
 C* RFTRIED. DFAS 336 
 
 C************************************************** DFAS 337 
 
 540 IF (JSTART .GT. OJ GO TO 548 DFAS 338 
 
 DO 544 I = l.NY DFAS 339 
 
 544 SAVE(2,l) = Y(2,I) DFAS 340 
 
 548 KFLAG « KFLAG - 2 DFAS 341 
 
 IF (H .LE. HMINI GO TO 740 DFAS 342 
 
 T(l) » TOLD OFAS 343 
 
 IF (KFLAG .LE. -5) GO TO 720 OFAS 344 
 
 550 PR2 * (D/E)**EN02*1.2 DFAS 345 
 
 L « DFAS 346 
 
 IF (NO .LE. 1) GO TO 570 OFAS 347 
 
 D - DFAS 348 
 
 DO 560 J ■ 1,M1 DFAS 349 
 
 560 D ■ D*( Y(K,J)/YMAX(J>)**2 DFAS 350 
 
 PR1 - (D/E0WN)«*EN01*1.3 DFAS 351 
 
 IF (PR1 .GE. PR2) GC TO 570 DFAS 352 
 
 PR2 * PR1 DFAS 353 
 
 L - -1 OFAS 354 
 
 570 IF (KFLAG .LT. .OR. NO . GE . MAXDERI GO TO 590 DFAS 355 
 
 D = OFAS 356 
 
 DO 580 J * 1,M1 DFAS 357 
 
 580 D - DM (ERROR(J)-ERSV( J) )/YMAX(J) )**2 DFAS 358 
 
 PR1 = (D/EUPI**EN03*1.4 DFAS 359 
 
 IF (PRl .GE. PR2) GC TO 590 DFAS 360 
 
 PR2 - PRl DFAS 361 
 
 L * 1 DFAS 362 
 
 590 R = 1.0/AMAXHPR2, l.E-5) DFAS 363 
 
 IF (KFLAG .LT. .OR. R .GE. 1.1) GO TO 600 DFAS 364 
 
 IDOUB - 9 DFAS 365 
 
 GO TO 700 OFAS 366 
 
36 
 
 600 NEWO = NO+L DFAS 367 
 
 K = NFWQ+1 DFAS 368 
 
 IF (NFWO .LE. NQI GO TO 620 DFAS 369 
 
 Rl = A(NFWQ)/FLDAT(NEV<Q) DFAS 370 
 
 DO 610 J = 1,NY DFAS 371 
 
 610 Y(K,J> = ERR0R(J)*P1 DFAS 37? 
 
 620 C"NTINUF DFAS 373 
 
 C***#*********#*****************************+****** DFAS 374 
 
 C* IF THE STEP WAS OKAY, SCALE THE Y VARIABLES IN ACCORDANCE DFAS 375 
 
 r* WITH THE NEW VALUE OF H. IF KFLAG < 0, HOWEVER, USE THE DFAS 376 
 
 C* SAVED VALUES (IN SAVE AND YLSV). IN EITHER CASE, IF THE DFAS 377 
 
 C* ORDER HAS CHANGED IT IS NECESSARY TO FIX CERTAIN PARAMETFRS DFAS 378 
 
 C* BY CALLING THE PROGRAM SEGMENT AT LINE 170. DFAS 379 
 
 r **** «•****»#** *****.*,*****************.******<.****:.!<**. DFAS 380 
 
 IDOUB = NO DFAS 381 
 
 IF (NEWO .EQ. NO) GO TO 630 DFAS 382 
 
 NO = NFWO D C AS 383 
 
 ASSIGN 630 TO IPET DFAS 384 
 
 GO TO 170 DFAS 385 
 
 630 IF (KFLAG .GT. 0) GO TO 670 DFAS 386 
 
 RACUM = RACUM*R DFAS 387 
 GO TO 750 ' DFAS 388 
 
 670 P = DMAX1(DMIN1(HMAX/H,R),HMIN/H) DFAS 389 
 
 H = H*R DFAS 390 
 
 IWEVAL = 1 DFAS 391 
 
 ASSIGN 700 TO IRET DFAS 392 
 
 GO TO 800 DFAS 393 
 
 700 DO 710 I = 1,M1 DFAS 394 
 
 710 YMAX(I) = DMAXK YMAXU ) ,DARS(Y( 1, I ) )) DFAS 395 
 
 GO TO 90 DFAS 396 
 
 C ********************************************* ***** DFAS 397 
 
 C* THE ERROR TEST HAS NOW FAILED THREE TIMES, SO THE DFAS 398 
 
 C* DERIVATIVES ARE IN BAD SHAPE. RETUPN TO FIRST-ORDER DFAS 399 
 
 C* METHOD AND TRY AGAIN. OF COURSE, IF NQ = 1 ALREAOY, DFAS 400 
 
 C* THEN THERE IS NO HOPE AND WE EXIT WITH KFLAG = -4. DFAS 401 
 
 C************************************************** DFAS 402 
 
 720 IF (NO .EO. I) GO TO 735 DFAS 403 
 
 NO = 1 DFAS 404 
 
 IDOUB = 1 DFAS 405 
 
 ASSIGN 751 TO IPET DFAS 406 
 
 GO TO 170 DFAS 407 
 
 735 NOOLD * 1 DFAS 408 
 
 KFLAG « -4 DFAS 409 
 
 GO TO 470 DFAS 410 
 
 740 KFLAG = -1 DFAS 411 
 
 GO TO 90 DFAS 412 
 
 r *»********•**•******»**•••*•*•**•.»•**»..».* .... »* DFAS 413 
 
 C* THIS SECTION RESTORES THE SAVED VALUES OF Y AND YL , SCALING DFAS 414 
 
 C* THE Y DERIVATIVES AS NECESSARY, AND THEN RETURNS TO THE DFAS 415 
 
 C* PREDICTER LOOP. DFAS 416 
 
 r ****+**********♦*♦******* *•*•*•«*•******* ********* DFAS 417 
 
 750 H ■ HOLD*RACUM DFAS 418 
 H ■ DMAX1(HMIN,0MIN1(H,HMAX) ) DFAS 419 
 
 751 RACUM - H/HOLO DFAS 420 
 Rl « 1.0 DFAS 421 
 DO 760 J ■ 2,K DFAS 422 
 
 Rl « R1*RACUM DFAS 423 
 
 00 760 I ■ 1,NY DFAS 424 
 
 760 Y(J,IJ ■ SAVE(J,II*Rl DFAS 425 
 
 DO 770 I - liNV DFAS 426 
 
 770 Yd, I) - SAVE(l.I) DFAS 427 
 
 CALL COPYZ (YL.YLSV.LCOPYll DFAS 428 
 
 IWEVAL - 1 DFAS 429 
 
 GO TO 250 DFAS 430 
 
 r I*.**.*.****.***.***.**..*.*.*.,,*.*.*****. ****... DFAS 431 
 
 C* THIS SECTION SCALES THE Y DERIVATIVES BY R. DFAS 432 
 
 r *********************************** *************** DFAS 433 
 
 800 Rl « 1.0 DFAS 434 
 
 DO 810 J - 2,K DFAS 435 
 
 PI * R1*R DFAS 436 
 
 DO 810 I » 1,NY OFAS 437 
 
 810 Y(J,I) « Y(J,I)*R1 DFAS 438 
 
 GO TO IRET, (100,700» DFAS 439 
 
37 
 
 C* 
 
 r ****** 
 
 ENTRY RFDSUB 
 
 ************ *****, 
 
 ♦•♦•♦A************************************** 
 
 PNTFR HERE TO RESET VALUES IN PREP FOP AN AUTO-RESTAPT . 
 I******************************************* 
 
 ENTRY RFDSUB 
 
 IF INTERRUPT OCCURPFD IN MATIN1, Nn RFSTORATIDN NECESSARY: 
 
 IF (KZILCH .EO. 0) RETURN 
 
 "0 910 J = 1,NY 
 
 Y( I, J) = SAVEIl, J> 
 910 Y(2,J) = SAVEI2, JJ/HOLD 
 
 CALL CPPYZ (YL,YLSV,LCOPYL) 
 T( 1) = TOLD 
 RETURN 
 ' END 
 
 DFAS 
 
 440 
 
 n^AS 
 
 441 
 
 OFAS 
 
 442 
 
 DFAS 
 
 443 
 
 PFAS 
 
 444 
 
 OFAS 
 
 445 
 
 DFAS 
 
 446 
 
 OFAS 
 
 447 
 
 DFAS 
 
 448 
 
 DFAS 
 
 449 
 
 DFAS 
 
 450 
 
 DFAS 
 
 451 
 
 DFAS 
 
 452 
 
BLIOGRAPHIC DATA 
 IEET 
 
 1. Report No. 
 
 UIUCDCS-R-T3-5T5 
 
 3. Recipient's Accession No. 
 
 'Title and Subtitle 
 
 DOCUMENTATION FOR DFASUB — A Program for the Solution 
 
 5. Report Date 
 
 of Simultaneous Implicit Differential and Nonlinear Equatior 6 
 
 July 1973 
 
 Author(s) 
 
 R. L. Brown, C. W. Gear 
 
 8. Performing Organization Rept. 
 No. 
 
 Performing Organization Name and Address 
 
 Department of Computer Science 
 
 University of Illinois at Urban a- Champaign 
 
 Urban a, Illinois 6l801 
 
 10. Project/Task/Work Unit No. 
 
 11. Contract /Grant No. 
 
 US AEC AT(ll-l)lU69 
 
 Sponsoring Organization Name and Address 
 
 US AEC Chicago Operations Office 
 9800 South Cass Avenue 
 Argonne, Illinois 
 
 13. Type of Report & Period 
 Covered 
 
 Technical 
 
 14. 
 
 Supplementary Notes 
 
 . Abstracts 
 
 This report surveys the theory of an efficient method for solving 
 .f(y_s y_» "t) + P(t)v = to get approximations to v_(t) and v(t) at 
 discrete time points t. > t given y_(t ) andy_'(t ). The 
 organization and use of a program, DFASUB, which uses this 
 theory is described. 
 
 Key Words and Document Analysis. 17o. Descriptors 
 
 ordinary differential equations 
 
 nonlinear equations 
 
 linear equations 
 
 discrete variable method 
 
 simultaneous implicit differential and nonlinear equations 
 
 b. Identifiers/Open-Ended Terms 
 
 e. COSATI Field/Group 
 
 '■•Availability Statement 
 
 unlimited 
 
 ,R M NTIS-38 (10-70) 
 
 19. Security Class (This 
 Report) 
 
 UNCLASSIFIED 
 
 20. Security Class (This 
 
 Page 
 UNCLASSIFIED 
 
 21. No. of Pages 
 
 k2 
 
 22. Price 
 
 USCOMM-DC 40329-P7 1 
 
orm AEC-427 
 
 (6/68) 
 AECM 3201 
 
 U.S. ATOMIC ENERGY COMMISSION 
 
 UNIVERSITY-TYPE CONTRACTOR'S RECOMMENDATION FOR 
 
 DISPOSITION OF SCIENTIFIC AND TECHNICAL DOCUMENT 
 
 ( See Instructions on Reverse Side ) 
 
 AEC REPORT NO. 
 
 C00-1U69-0225 
 
 2. TITLE 
 
 DOCUMENTATION FOR DFASUB— A Program for the Solution 
 of Simultaneous Implicit Differential and Nonlinear 
 Equations 
 
 i. TYPE OF DOCUMENT (Check one): 
 
 P a. Scientific and technical report 
 
 2] b. Conference paper not to be published in a journal: 
 
 Title of conference 
 
 Date of conference 
 
 Exact location of conference _ 
 
 Sponsoring organization 
 
 fjc. Other (Specify) 
 
 I RECOMMENDED ANNOUNCEMENT AND DISTRIBUTION (Check one): 
 
 [3 a. AEC's normal announcement and distribution procedures may be followed. 
 
 ~ J b. Make available only within AEC and to AEC contractors and other U.S. Government agencies and their contractors. 
 
 ^] c. Make no announcement or distribution 
 
 . REASON FOR RECOMMENDED RESTRICTIONS: 
 
 SUBMITTED BY: NAME AND POSITION (Please print or type) 
 
 C. W. Gear 
 
 Professor and Principal Investigator 
 
 Organization 
 
 Department of Computer Science 
 
 University of Illinois at Urbana-Champaign 
 
 Urbana, Illinois 6l801 
 
 Signature 
 
 ^frj^tofo^. 
 
 Date 
 
 July 1973 
 
 FOR AEC USE ONLY 
 
 AEC CONTRACT ADMINISTRATOR'S COMMENTS, IF ANY, ON ABOVE ANNOUNCEMENT AND DISTRIBUTION 
 RECOMMENDATION: 
 
 PATENT CLEARANCE: 
 
 _l a. AEC patent clearance has been granted by responsible AEC patent group. 
 _| b. Report has been sent to responsible AEC patent group for clearance. 
 U c. Patent clearance not required. 
 
A* 
 
 \ 
 
 \ 
 
 d>