Hi ffi 
 
 ■MB 
 
 VBV 
 
 Hr 
 mm 
 
 M 
 
 ■ mm 
 
 Hfl 
 
 BM 
 
 H 
 
 BNiEl 
 
 m maSBSSt 
 
 WKttSBSSsM 
 
 11 
 
 WSSSSM 
 
 MIM 
 
IIUCDCS-R-76-839 
 
 >%^Ci 
 
 COO-2383-0035 
 
 X 
 
 \tJ3f 
 
 Future Developments in Stiff Integration Techniques - 
 Stability of Methods that do not use an Exact Jacobian 
 
 by 
 C. W. Gear 
 
 lecember 1976 
 
UIUCDCS-R-76-839 
 
 Future Developments in Stiff Integration Techniques - 
 Stability of Methods that do not use an Exact Jacobian* 
 
 by 
 
 C. W. Gear 
 
 December 1976 
 
 Department of Computer Science 
 
 University of Illinois at Urbana-Champaign 
 
 Urbana, Illinois 61801 
 
 *Written while on leave at Yale University, Department of Computer Science. 
 
 Tbis work supported in part by the US ERDA under grant US ERDA/EY-76-S-02-283 and 
 
 US AFOSR under grant 75-2854. 
 
Digitized by the Internet Archive 
 in 2013 
 
 http://archive.org/details/futuredevelopmen839gear 
 
1. Abstract 
 
 Although there are many open theoretical questions for stiff (and 
 non-stiff) variable order, variable step methods, it seems unlikely that 
 there will be any major computational advances as far as integration methods 
 for the general non-linear problem. The major drawback of many of the 
 methods that have been proposed in the last few years is that they involve 
 too many matrix operations, or operations on larger matrices than in some 
 of the more straightforward techniques such as the backward differentiation 
 formulas. However, advances in techniques for handling large, sparse 
 matrices may affect the choice of integration method for very large problems, 
 This paper examines possible improvements that might occur through the use 
 of very inaccurate Jacobians. Such Jacobians might be chosen so that the LU 
 decomposition is particularly simple. Some methods are shown to remain 
 stable in spite of very large errors in the Jacobian. These results are 
 preliminary, but suggest that further speed increases in very large problems 
 are possible. 
 
 2. Introduction 
 
 Although a large number of A-stable methods are known (for example, 
 implicit Runge-Kutta - Butcher [3], multiderivative - Enright [6], Genin [9], 
 Brown [2], extrapolation - Lindberg [16], and other clever techniques - 
 Bickart and Picel [l], Liniger and Odeh [17], Rosenbrock [22], Skeel and 
 Kong [24], Sloate and Bickart [25], and Watanabe [26]) the method used in 
 the majority of large codes is based on the backward differentiation formulas 
 which are not A-stable for order greater than two, and which do not have 
 very small error coefficients. The principal reason that this is true, apart 
 
from the existence of fairly robust automatic software incorporating these 
 methods but not the other methods (Hindmarsh [10], Hindmarsh and byrne [ll], 
 Krogh [14], Shampine [23], and Gear [8]), is that most other methods either 
 require exact values of the Jacobians at every step in order to achieve 
 their accuracy (e.g. multiderivative methods and t'le Rosenbrock-like methods), 
 or require that a larger system of non-linear equations be solved at each 
 step, or both. While neither of these requirements presents a particular 
 problem for small systems, small systems do not consume much computer time 
 so there is little incentive to consider other methods if there is a 
 reasonable subroutine in the library, particularly if it calls for no more 
 than a minimum amount of programming on the part of the user. However, in 
 large systems, such as those that arise in electrical network analysis, 
 chemical kinetic problems coupled with fluid flow considerations, 
 lasing, and spatial discretization of a system of partial differential 
 equations to get a system of ordinary differential equations, significant 
 amounts of time are spent in solving non-linear equations by repeated 
 solution- of linear problems in a quasi-Newton iteration. The figure of 
 30% of the total solution time has been observed for systems of less than 
 a thousand ODEs, and this figure increases with n. Considerable savings 
 could be realized if better ways to handle the non-linear systems could be 
 developed. Further, since, for problems which require modest accuracy and 
 for which function evaluations do not cos.t much, one-step methods such as 
 implicit Runge-Kutta methods may be superior, it may be that these methods 
 will turn out to be superior when better matrix methods are developed. 
 The improvement may be particularly significant in problems whose solution 
 displays repeated transient regions (in which the solution varies rapidly) 
 
 2 
 
separated by relatively quiescent periods in which the solution is slowly 
 varying, as can happen in lasing problems. The reason for the potential 
 advantage of one-step methods is that in variable step methods, the step 
 gets large in the quiescent regions, and then a sudden transient looks , 
 very much like a discontinuity that cannot be handled by current multistep 
 methods. 
 
 In this paper we will suggest some approaches that may lead to 
 improved techniques for handling the non-linear systems of equations. This 
 is a very preliminary report of ongoing work (which may encounter insuper- 
 able difficulties and have no value) but it is hoped that the report will 
 indicate some of the problems and possible solutions. 
 
 3. Why is a Jacobian needed? 
 
 Consider the differential equation 
 
 y f = a.[y - u(t)] + u'(t) (1) 
 
 whose solution is 
 
 y = u(t) + c.exp(at) 
 If a << in the time scale in which u(t) is slowly varying, the problem 
 is stiff. In fact, as a tends to minus infinity, the solution tends to 
 u(t) with a boundary layer at t = in the case that y(0) # u(0). We can 
 also write (1) as the singular perturbation problem 
 
 ey ' - eu* = y - u(t) = g(y,t) (2) 
 
 where e = 1/a. We see that for very small e we are really faced with the 
 solution of the equation g(y,t) = at each time step. In the case of a 
 non-linear differential equation we observe that if the Jacobian J = 3f/3y 
 is very "large" (that is, has very negative eigenvalues), small deviations 
 
 3 
 
of y from the value that makes f(y,t) almost zero lead to very large 
 derivatives which will cause the solution to be returned near to the 
 solution of f(y,t) = very rapidly (or the problem is unstable). Conse- 
 quently whatever we do must amount to solving the non-linear problem 
 f(y,t) = approximately. Many methods for the aporoximate solution of 
 non-linear equations make some use of the Jacobian. The following simple 
 predictor-corrector scheme for the problem y 1 = f(y,t) illustrates the use 
 
 of the Jacobian. Let y and y' be approximations to y and its derivative 
 
 n n J 
 
 at t = t . A single step consists of 
 
 Euler predictor p = y + 1 y' 
 
 n+1 n n 
 
 One-step corrector y = y + bhy' + (1 - b)hy' 
 
 n+1 n n+1 n 
 
 (3) 
 
 If b :> .5 this is good for stiff problems. Since the corrector is implicit 
 
 (as are all methods that have good stability properties for stiff problems) , 
 
 we must decide how to solve for y . If we have an approximation A to the 
 
 J n+l rK 
 
 Jacobian 3f/3y, we can use one step of a Newton method starting from the 
 predictor (as is done in most current codes) to get 
 
 ^n +1 " P n+ 1 - [I " bhAl "' [ P„ +1 " *n " !*«W " (1 - h) < ] <4) 
 
 (The dependence of f on t is not explicitly indicated where unimportant.) 
 If f(y) = Ay, this solves the corrector in one step. If A is reasonably 
 close to the Jacobian J of the system, the method may still have good 
 stability properties. (See Klopfenstein and Davis [13] who consider some 
 PECE methods which are stable for large values of J provided that A is 
 close to J in a sense that requires the difference to remain of order 1.) 
 To see what happens when the Jacobian of a system is very large, consider 
 a sequence of systems with Jacobian c J . Replace J with cJ in (4) and take 
 
 4 
 
the limit as c goes to infinity. Noting that f(y,t) will have to be 
 replaced with cf(y,t) we get 
 
 y ■ P - J~ l f(p , t ) (5) 
 
 y n+i n+i vp n+i ' n+i K J 
 
 which is a Newton step for the solution of f(y,t) = 0. This is no surprise, 
 but it has been introduced to motivate consideration of some other tech- 
 niques for the solution of linear and non-linear equations. First let us 
 consider ways of solving the system g(x) =0. We will assume that all of 
 the eigenvalues of g are strictly in the right half-plane. If g is linear 
 and there is some knowledge of the location of these eigenvalues, there are 
 a number of iterative methods that can be used. The recent thesis of 
 Manteuffel [18] discusses a Chebyshev type of iteration that attempts to 
 estimate the location of the critical eigenvalues, and references much of 
 the literature in the general area. However, let us start with the simplest 
 (and usually useless) iterative technique possible, functional iteration. 
 The iteration 
 
 x , = x - g(x ) 
 m+1 m B m 
 
 only converges if we can bound the norm of I - J to be less than one, where 
 
 J is the Jacobian of g (that is, if the eigenvalues of J are strictly inside 
 
 the unit circle centered at 1) . Since that is an unreasonable restriction^! 
 
 we look for a circle of radius R centered at R that contains all of the 
 
 eigenvalues of J. If they are all in the positive half -plane, such clearly 
 
 exists. Then we can use the iteration 
 
 x = x - R" 1 g(x ) (6) 
 
 m+i m m 
 
 It converges, albeit very slowly. We can view equation (6) as a 
 
quasi-Newton method where R is an approximation to the Jacobian. In fact, 
 the rate of convergence is 1 - d/R where d is the smallest distance from 
 the radius R circle centered at R to any eigenvalue of J. However, we are 
 going to examine the application of a single iteration of this type to the 
 corrector of a differential equation because we hope that we can get 
 accuracy with a good predictor, and use the corrector only to achieve 
 stability. Later, we will examine other iterative techniques that make use 
 of a greater knowledge of the Jacobian than a bound on the eigenvalues. 
 
 4. Regions of Absolute Stability 
 
 The region of absolute stability of a method is that set of values 
 of ha such that the numerical solution of y' = ay decays when step size h 
 is used. A method is A-stable if the region of absolute stability includes 
 all of the negative half-plane. Many methods are not A-stable, but usually 
 if the eigenvalues of the Jacobian are inside the region of absolute 
 stability we will get reasonable error behavior. (See Odeh and Liniger [19] 
 and Dahlquist [4], [5].) The region of absolute stability for the corrector 
 in (3) is as shown in Figure 1. If b i .5 the method is A-stable. However, 
 suppose that we try to "solve" the corrector with a single modified 
 functional iteration as given in (4) using -R as an "approximation" to A in 
 analogy with (6). (Strictly speaking, we must use -R times the identity as 
 the approximation.) If we now consider the behavior of the linear scalar 
 test equation y' = ay we get 
 
 y = [1 + ha + [1 - bhR] _1 b(ha) 2 ]y 
 
 or 
 
y . = [1 - bhR]" [1 + h(a - bR) + bh 2 a(a - R)]y 
 n+l n 
 
 2 
 Because h a(a - R) appears in the numerator, while only hR appears in the 
 
 denominator, this cannot be stable for other than a small region of the ha 
 
 plane (of order 1). The problem occurs because the calculation of y' 
 
 introduces a factor a. This happens because we are doing the equivalent of 
 
 a PECE method. Suppose instead we look at a PEC method. If y and d are 
 
 n n 
 
 the numerical approximations to y(t ) and hy'(t ), we can write a PEC 
 
 n n 
 
 method in the form 
 
 Predictor 
 
 Evaluation 
 
 Corrector 
 
 n+l n n 
 
 d 4., = " hR y 4., + t hR P J., + hf (P 4.,)] 
 
 n+l n+l n+l n+l 
 
 y = y + bd , + (1 - b)d 
 
 n+l n n+l n 
 
 (7) 
 
 Here we have replaced the Jacobian by -R times the identity as before. The 
 behavior of this process for the scalar linear problem y' = ay is described 
 by 
 
 -1 
 
 n+l 
 
 n+l 
 
 n+l 
 
 -b 
 
 - [1 + bhR] 
 
 hR 
 
 -1 
 
 1 - b 
 
 h(a + R) h(a + r) 
 
 
 y n 
 
 
 d 
 n_ 
 
 1 + bh(a + R) 1 - b + bh(a + R) 
 
 ha 
 
 ha + bhR 
 
 
 *-. — 
 
 
 y n 
 
 
 d 
 n 
 
 = Sw 
 
 (8) 
 
 The numerator and denominator of all elements of S are linear in ha and 
 hR, so it is possible that S is power bounded for large values of a (in 
 particular, up to values of order R) . In fact, that is the case if b ^ .5. 
 The region of absolute stability of (7) can be found by plotting the locus 
 
 7 
 
of values of ha for which S has a unit eigenvalue for a fixed value of hi'. 
 It has the form shown in Figure 2 provided that b >■ .5. 
 
 The interesting aspect of this is that the region of absolute 
 stability can be made arbitrarily large by choosing R large enough. This 
 means that equation (7) is a first order (second order if b = .5) method 
 that has a very large region of absolute stability. Unfortunately, very 
 
 little comes for free, and an examination of the error coefficients reveals 
 
 1/2 
 a factor of (hR) in the error growth. To see this, apply equation (7) 
 
 to equation (1) , setting the global error e to be a vector with components 
 
 y - u(t ) and d - hu'(t ). A little computation shows that 
 n n n n 
 
 e , = Se + (1 + bhR) 
 n+1 n 
 
 + (1 + bhR) 
 
 -1 
 
 2b - 1 - bh(R + a) 
 
 _-hR(2b - 1) - h(R + a^ 
 3b - 1 - bh(R + a) 
 
 h 2 u"/2 
 
 -1 
 
 h 3 u m /6 + OCh 4 ) (9) 
 
 -hR(3b - 1) - h(R + a) 
 where S is defined in equation (8). The truncation error, which is the last 
 three terms in (9) , is bounded as R gets large because of the R in the 
 denominator, and S is stable, so the solution of (9), which is 
 
 N-l 
 
 _N-n-l Q N 
 e„ = l b q + b e 
 n 
 n=o 
 
 'N 
 
 (10) 
 
 where q is the truncation error in the n-th step, is bounded. However, the 
 n 
 
 N 
 bound on S involves a factor R. This can be seen by considering the form 
 
 of S for very large hR. It tends to the matrix 
 
 1 1 
 
 1 
 
which is unstable. To find the growth as a function of hR we must look at 
 the matrix more closely. We show in lemma 1 below that the matrix Q which 
 takes S to a diagonal form via a similarity transformation is such that 
 
 -1 
 
 1/2 
 
 jQ -|| MqII ■- 0(hR) 
 can then deduce from (10) that 
 
 If we assume that q is bounded by Kh , we 
 
 i i -i n i i 
 
 ||e N || * ||s N ||MKh 2 + ||S N 
 
 -1 
 
 £ (KNh + lie II) I | Q 
 = 0((Nh 2 + | |e ||) (hR) 
 
 1/2 
 
 (ID 
 
 ,-1 
 
 LEMMA 1 If S is as given in equation (8) , hR is large, and Q S Q is 
 
 — 1 1/2 
 
 diagonal, then ||Q | |q|) = 0(hR) is the smallest bound possible. 
 PROOF The eigenvalues of S are not identically one when hR is finite. 
 When hR is finite, det(zl - S) = (z - l) 2 , but when hR is finite, 
 det(zl - S) is a monic polynomial in z with coefficients that are linear 
 in (1 + bhR) because 
 
 S = 
 
 "l 
 
 - — 1 
 
 1 
 
 + (1 + bhR) 1 
 
 "b 
 
 
 
 1 
 
 
 1 
 
 [ha ha - f] 
 
 = A + ( 1 + bhR)" 1 B 
 where B is a rank one matrix. Furthermore, since the constant term in 
 
 det(zl - S) is det(-S) which is not identically one, we can conclude that 
 
 -1/2 
 the eigenvalues of S are 1 + 0((hR) ) for large hR. (This can be seen 
 
 by expanding the eigenvalues as as power series in an arbitrary power of 
 
 hR.) The columns of Q are scalar multiples of the eigenvectors of S, and 
 
 T 
 it is easy to see that the eigenvectors of S have the form (1 p) , 
 
 l/2 -1 
 
 where p = 0(hR) . From this we see that the rows of Q have the form 
 
(1 1/p) • No change of scale will avoid the term 1/p, so we conclude 
 that ||Q _1 | I ||q|| = 0(p) ■ 0((hR) 1/2 ). QED 
 
 The analysis above shows that if a, the eigenvalue of the method, 
 
 1/2 
 is small compared to R, the error has a component (hR) . However, if a 
 
 is of the order of -R, then a similar analysis shows that the error does 
 
 not depend on such a term. This means that in a system of equations with 
 
 several eigenvalues, we will get large error terms in the components with 
 
 small eigenvalues. This suggests that in equation (7) we should use a 
 
 matrix A in place of R that has large eigenvalues corresponding to the 
 
 large eigenvalues of J, but small (or zero) eigenvalues corresponding to 
 
 reasonable size eigenvalues of J. To put this precisely, we can show that 
 
 if such -an A were simultaneously diagonalizable with J and the eigenvalues 
 
 of the two matrices corresponded in this way, the method would be both 
 
 1/2 
 stable and have an error term that did not include the (hR) term. 
 
 Unfortunately, such matrices are not necessarily "simpler" than 
 
 the original matrix J (by which we mean that they not necessarily easier 
 
 to LU decompose). It appears to be important that A have eigenvectors 
 
 that are close to the eigenvectors of J corresponding to the large 
 
 eigenvalues. (See, for example, Lambert [15].) 
 
 5. Higher Order Methods 
 
 The analysis can be extended to higher order methods. Consider 
 the PEC scheme 
 
 10 
 
p J. = £ 
 
 n+l p 
 
 d . = -hRy , + hRp ^ + hf (p . ) (12) 
 
 n+l 'n+l *n+l Vl 
 
 y , =bd , + I 
 J n+1 n+l c 
 
 where E and Z are linear combinations of past values of y and hy' for 
 p c 
 
 the predictor and corrector respectively. By considering y' = ay, we see 
 that 
 
 y , , - (1 + bhR)" 1 [Z + bh(R + a)Z ] 
 n+l c p 
 
 d ,, = (1 + bhR) [h(R + a)I - hRZ ] 
 n+l p c 
 
 Some straightforward calculation on this shows that if the value of ha for 
 which an eigenvalue of the corrector process alone is y is written as 
 ha = u (y) » and if the value of ha for which an eigenvalue of the conven- 
 tional predictor corrector PEC process is y is written as ha = u (y) » then 
 
 pc 
 
 the value of ha for which an eigenvalue of the process in equation (13) is 
 Y can be written as 
 
 ha = u( Y ) = u ( Y ) + hR[ P ° ( > ° ] (14) 
 
 c 
 
 The region of absolute stability of equation (12) can be found by examining 
 
 its boundary which is included in the set of values of ha given by equation 
 
 i ft 
 
 (14) with y = e . For many pairs of predictors and correctors, this region 
 can be made arbitrarily large by increasing R. 
 
 Unfortunately, a study of the error term similar to that in the 
 
 k/(k+l) 
 previous section shows that a factor (hR) appears in the error growth 
 
 rate, where k is the step number of the method, so the higher order may not 
 
 help too much. 
 
 11 
 
6. Tentative Conclusioas 
 
 The stability of these methods indicates possible approaches. in 
 particular, two avenues seem to be worth exploring. One is to attempt to 
 use the ideas of Chebyshev iteration which can hasten the convergence of 
 iterative methods for linear equations (see Manteuffel [18]), and the other 
 is to consider using different approximations R to the Jacobian J on dif- 
 ferent steps. Such approximations could be chosen to "beat down" the 
 eigenvectors corresponding to the large unwanted components in the Jacobian 
 in the spirit of Lambert [15]. In fact, the alternating direction method 
 is a form of this (see Richtmeyer and Morton [21], p. 176). In the ADI 
 method, an approximation is chosen that has the same eigenvectors as J, and 
 alternately damps the large eigenvectors corresponding to the x and y 
 spatial directions. 
 
 7. Acknowledgments 
 
 This work was supported in part by the US ERDA under grant 
 AT (11-1) 2383, by the US AFOSR under grant 75-2854, and by Yale University 
 Department of Computer Science. The author would like to thank Pam Farr 
 for manuscript preparation, and to acknowledge especially Professor 
 Eisenstat of Yale University for many helpful discussions. 
 
 8. References 
 
 [ 1] Bickart, T., Picel, P., High Order Stiffly Stable Composite Multistep 
 Methods for Numerical Integration of Stiff Differential Equations , 
 1973 Report, EE and Comp . E. Dept., Syracuse University. 
 
 12 
 
[ 2] Brown, R. L . , Multiderivative Methods for the Solution of Stiff 
 
 Ordinary Differential Equations , Report UIUCDCS-R-74-672 , Dept. of 
 Computer Science, University of Illinois at Urbana-Champaign, 1974. 
 
 [ 3] Butcher, J. C, "Implicit Runge-Kutta Processes," Math . Comp . , 
 Vol. 18, pp. 50-64. 
 
 [ 4 Dahlquist, G. , On the Stability and Error Analysis for Stiff 
 
 Non-linear Systems . Part I, Report TRITA-NA-7 508, Dept. of Inf. 
 Processing, Royal Inst, of Technology, Stockholm, 1975. 
 
 [ 5J Dahlquist, G., Error Analysis for a Class of Methods for Stiff 
 
 Non-linear Initial Problems . Report TRITA-NA-7511, Dept. of Inf. 
 Processing, Royal Inst, of Technology, Stockholm, 1975. 
 
 , 6] Enright, W. H., Studies in the Numerical Solution of Stiff Ordinary 
 Differential Equations , Report #46, Dept. of Computer Science, 
 University of Toronto, Canada, October, 1974. 
 
 [ 7] Enright, W. H. , Hull, T. E. , Lindberg, B. , "Comparing Numerical 
 
 Methods for Stiff Ordinary Differential Equations," BIT , Vol. 15, 
 pp. 10-48, 1975. 
 
 [ 8] Gear, C. W. , "DIFSUB for the Solution of Ordinary Differential 
 Equations," Comm . ACM , Vol. 14, pp. 185-190, 1971. 
 
 [ 9] Genin, Y., An Algebraic Approach to A-Stable Linear Multistep- 
 
 Multiderivatjve Integration Formulas , Report R244, MBLE Research 
 Lab, Brussels, Belgium, February, 1974. 
 
 [l0] Hindmarsh, A. C. , GEAR: Ordinary Differential Equation System Solver , 
 Lawrence Livermore Laboratory, Report UCID-30001, Rev. 3, December, 
 1974. Program available from Argonne Code Center, Argonne National 
 Laboratory, Argonne, IL 60439. 
 
 [ll] Hindmarsh, A. C, Byrne, G. D., EPISODE: An Experimental Package 
 for the Integration of Systems of Ordinary Differential Equations , 
 Lawrence Livermore Laboratory, Report UCtD-30112, May, 1975. Program 
 available from Argonne Code Center, Argonne National Laboratory, 
 Argonne, IL 60439. 
 
 [12] Hull, T. E., Enright, W. H. , Fellen, B. M. , Sedgwick, A. E., 
 
 "Comparing Numerical Methods for Ordinary Differential Equations," 
 SIAM J. Numerical Anal . , Vol. 9, pp. 603-637, 1972. 
 
 [13] Klopfenstein, R. W. , Davix, C. B., "PECE Algorithms for the Solution 
 of Stiff Systems of Differential Equations," Math . Comp . , Vol. 25, 
 pp. 457-473, 1971. 
 
 13 
 
[14] Krogh, F. T. , VODQ/SVDQ/nVPQ - Variable Order Integrators for the 
 Numerical Solution of Ordinary Differential Kquations , TU L)oc . 
 CP-2308, NPO-11643, Jet Propulsion Laboratory, Pasadena, CA , May, 
 1969. 
 
 [15] Lambert, J. D. , The Numerical Integration of a Special Class of 
 Stiff Systems , Report //13, Dept. of Mathematics, University of 
 Dundee, Scotland, November, 1975. 
 
 [16] Lindberg, B. , "On Smoothing and Extrapolation for the Trapezoidal 
 Rule," BIT, Vol. 11, pp. 29-52, 1971. 
 
 [17] Liniger, W. , Odeh, F., "A-Stable Accurate Averaging of Multistep 
 Methods for Stiff Differential Equations," IBM Journ. of_ Res. and 
 Dev . , Vol. 16, pp. 335-347, 1972. 
 
 [18] Manteuffel, T. A., An Iterative Scheme for Solving Non-Symmetric 
 Linear Systems with Dynamic Estimation of Parameters , Report 
 UIUCDCS-R-75-758, Dept. of Computer Science, University of Illinois 
 at Urbana-Champaign, October, 1975. 
 
 [19] Odeh, F., Liniger, W. , Non-Linear Fixed-h Stability of Multistep 
 Formulae , IBM Research Report RC 5717, Yorktown Heights, NY, 
 November, 1975. 
 
 [20] Prothero, A., Robinson, A., "On the Stability and Accuracy of One- 
 Step Methods for Solving Stiff Systems of Ordinary Differential 
 Equations," Math . Comp . , Vol. 28, pp. 145-162, 1974. 
 
 [2l] Richtmeyer, R. D., Morton, K. W. , Difference Methods for Initial- 
 Value Problems , Wiley Interscience , New York, 1967. 
 
 [22] Rosenbrock, H. H., "Some General Implicit Processes for the Numerical 
 Solution of Differential Equations," Comp . J_. , Vol. 5, pp. 329-330, 
 1963. 
 
 [23"] Shampine, L. F. , Gordon, M. K., Computer Solution of Ordinary 
 
 Differential Equations , W. H. Freeman and Co., San Francisco, CA, 
 1975. 
 
 [24] Skeel, R. D., Kong, A. K. , Blended Multistep Methods , in preparation, 
 Dept. of Computer Science, University of Illinois at Urbana-Champaign, 
 Urbana, IL 61801, 1976. 
 
 [25] Sloate, H., Bickart, T. , "A-Stable Composite Multistep Methods," 
 J ACM , Vol. 20, pp. 7-26, 1973. 
 
 [26] Watanabe, D. S., "Block Implicit One-Step Methods," to appear, 
 Math . Comp . 
 
 [27] Willoughby, R. A., ed., Stiff Differential Systems , Plenum Press. 
 
 14 
 
ha-plane 
 
 unstable 
 
 stable 
 
 Figure 1 
 
 Figure 2 
 15 
 
 ha-plane 
 
FormAEC-427 U.S. ATOMIC ENERGY COMMISSION 
 
 pXi 6 ?™ UNIVERSITY-TYPE CONTRACTOR'S RECOMMENDATION FOR 
 
 DISPOSITION OF SCIENTIFIC AND TECHNICAL DOCUMENT 
 
 ( S*e Instructions on Rmnrse Sid* ) 
 
 1. AEC REPORT NO. 
 
 COO-2383-0035 
 
 2. TITLE 
 
 Future Developments in Stiff Integration Techniques - 
 Stability of Methods that do not use an Exact Jacobian 
 
 3. TYPE OF DOCUMENT (Check one): 
 
 a. Scientific and technical report 
 
 [~l b. Conference paper not to be published in a journal: 
 
 Title of conference 
 
 Date of conference 
 
 Exact location of conference. 
 
 Sponsoring organization 
 
 □ c. Other (Specify) 
 
 4. RECOMMENDED ANNOUNCEMENT AND DISTRIBUTION (Check one): 
 
 Ixl a. AEC's normal announcement and distribution procedures may be followed. 
 
 I | b. Make available only within AEC and to AEC contractors and other U.S. Government agencies and their contractors. 
 
 f~| c. Make no announcement or distrubution. 
 
 5. REASON FOR RECOMMENDED RESTRICTIONS: 
 
 6. SUBMITTED BY: NAME AND POSITION (Please print or type) 
 
 C. W. Gear 
 
 Professor and Principal Investigator 
 
 Organization 
 
 Department of Computer Science 
 
 University of Illinois at Urbana-Champaign 
 
 Urbana, IL 61801 
 
 Signature ^4^%^ f\ < a 
 
 Date 
 
 December 1976 
 
 FOR AEC USE ONLY 
 
 7. AEC CONTRACT ADMINISTRATOR'S COMMENTS, IF ANY, ON ABOVE ANNOUNCEMENT AND DISTRIBUTION 
 RECOMMENDATION: 
 
 8. PATENT CLEARANCE: 
 
 Cj a. AEC patent clearance has been granted by responsible AEC patent group. 
 LJ b. Report has been sent to responsible AEC patent group for clearance. 
 I) c. Patent clearance not required. 
 
'bibliographic data 
 
 SHEET 
 
 1. Report No. 
 
 UIUCDCS-R-76-839 
 
 4. Title and Subtitle 
 
 Future Developments in Stiff Integration Techniques - 
 Stability of Methods that do not use an Exact Jacobian 
 
 3. Recipient's Accession No. 
 
 5. Report Date 
 
 December 1976 
 
 7. Author(s) 
 
 C. W. G 
 
 ear 
 
 8. Performing Organization Rept. 
 
 No - UIUCDCS-R-76-839 
 
 9. Performing Organization Name and Address 
 
 Department of Computer Science 
 
 University of Illinois at Urbana-Champaign 
 
 Urbana, IL 61801 
 
 10. Project/Task/Work Unit No. 
 
 COO-2383-0035 
 
 11. Contract/Grant No. 
 
 US ERDA/EY-76-S-02-283 
 
 12. Sponsoring Organization Name and Address 
 
 United States Energy Research and Development Administration 
 9800 South Cass Avenue 
 Argonne, IL 60439 
 
 13. Type of Report & Period 
 Covered 
 
 Research 
 
 14. 
 
 15. Supplementary Notes 
 
 16. Abstracts 
 
 Although there are many open theoretical questions for stiff (and non-stiff) variable 
 order, variable step methods, it seems unlikely that there will be any major computa- 
 tional advances as far as integration methods for the general non-linear problem. 
 The major drawback of many of the methods that have been proposed in the last few 
 years is that they involve too many matrix operations, or operations on larger 
 matrices than in some of the more straightforward techniques such as the backward 
 differentiation formulas. However, advances in techniques for handling large, 
 sparse matrices may affect the choice of integration method for very large problems. 
 This paper examines possible improvements that might occur through the use of very 
 inaccurate Jacobians. Such Jacobians might be chosen so that the LU decomposition 
 is particularly simple. Some methods are shown to remain stable in spite of very 
 large errors in the Jacobian. These results are preliminary, but suggest that 
 further speed increases in very large problems are possible. 
 
 17. Key Words and Document Analysis. 17a. Descriptors 
 
 stiff equations 
 
 17b. Identifiers/Open-Ended Terms 
 
 17e. COSATI Field/Group 
 
 18. Availability Statement 
 
 unlimited 
 
 PORM NTIS-35 (10-70) 
 
 19. Security Class (This 
 Report) 
 
 UNCLASSIFIED 
 
 20. Security Class (This 
 
 Page 
 UNCLASSIFIED 
 
 21. No. of Pages 
 19 
 
 22. Price 
 
 USCOMM-DC 40329-P7I 
 
m ?\m