ililli HHHBHBb BhHHBB ^H8_ HHrajHI ■ Bra 3ss&sy IBiBQHBQKQB SB SKI LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 5IO. 84 too. 800 - 8oe> Digitized by the Internet Archive in 2013 http://archive.org/details/blendedlinearmul800skee UIUCDCS-R-76-800 C00-2383-0030 / BLENDED LINEAR MULTISTEP METHODS by Robert D. Skeel and Antony K. Kong June 1976 UIUCDCS-R-76-800 C00-2383-0030 BLENDED LINEAR MULTISTEP METHODS by Robert D. Skeel and Antony K. Kong June 1976 DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN URBANA, ILLINOIS 61801 This work was supported in part by the Energy Research and Development Administration under contract US ERDA AT(ll-l) 2383. ABSTRACT The accuracy of linear multistep formulas suitable for stiff differential systems is limited. Greater accuracy can be attained by including higher derivatives in the formula, but this is not practical for all problems. However it is possible to duplicate the absolute stability region for any given m-derivative multistep formula by taking a linear combination of m multistep formulas. These blended formulas are similar to Lambert and Sigurdsson's linear multistep formulas with variable matrix coefficients, but the approach is different. Imple- mentation details and numerical results are presented for a variable- order, variable-step blend of the Adams-Moulton and the backward differentiation formulas. 910. 1. INTRODUCTION Currently the most popular codes for solving systems of stiff ordinary differential equations are based on the backward differentiation formulas, for example Gear [9], Bray ton et al [ 1 ], Hindmarsh [12], and Byrne and Hindmarsh [3]. The biggest drawback [7, p. 33] of these formulas is the poor stability properties of the higher order formulas when the Jacobian matrix for the problem has eigenvalues close to the imaginary axis. This difficulty can be overcome by restricting the order of the formulas; however these lower order formulas are often fairly inefficient because of their limited accuracy. There are conceivably three situations in which more accurate stiffly stable formulas would be helpful: (i) for nonstiff problems. The casual user of a program from a subroutine library should not be reauired to know whether or not his problem is stiff. And there is presently no auick and reliable way for automatically determining whether or not a problem is stiff. Therefore it may be reasonable to provide a stiff formula as the default. (I f efficiency were very imnortant, the user should not be considered as "casual".) (ii) for stiff problems with moderate to high accuracy requirements, (iii) for stiff problems whose Jacobian matrix has large eigenvalues near the imaginary axis. These are the problems for which the backward differentiation formulas are often inadeouate. One can get more accurate stiffly stable formulas for a system &-f(t.y> of s differential equations by including derivatives of f in the formula: k m . k f . n z o. y M . = z h J z b •• f . ' 1=0 n n_1 J-l 1-0 J1 n " n where f^ = f ^ n -1 ' y n-i ^ Here f denotes the total derivative of f defined recursively by f<°> = f , f U+l) = 9_ f (£) + f 3_ f U) at 3y Formulas of this type have been proposed by Enright [5] and Rrown [2]. The presence of higher derivatives makes it necessary to reauire the user to provide routines for their evaluation or better still to have these routines generated symbolically from f(t, y). Unfortunately there are differential equations for which symbolic differentiation is not possible, and in these cases the user may be unwilling (see [18]) or unable to provide analytic Jacobians. Thus the biggest drawback of these multi derivative formulas is their limited applicability. Stiff problems require implicit formulas, and hence systems of nonlinear equations must be solved. This is usually accomplished by some kind of modified Newton iteration, which means that an approximation to the Jacobian is available during the integration of a stiff system. It may be advantageous to make further use of the Jacobian in our method. We cannot use it to improve the order of accuracy of the method because it may not be exact, but we can use it to enhance the stability properties. This idea has been pursued by Lambert [15] and by Lambert and Sigurdsson [16], who introduce the class of linear multistep formulas with variable matrix coefficients. The formulas presented in this paper are of this type, but the approach is different. The basic idea is to take a linear mul tiderivati ve multistep formula and convert it into a variable coefficient linear multistep formula so that the region of absolute stability and the order of accuracy is unaffected. There is no reason to suspect that this transformation would either improve or degrade the performance of the formula in the case of smooth nroblems, although it should be beneficial for problems with nonexistent higher derivatives . As an example, the (k + 1 )st order Enri ght 'formula [5] is transformed into a blend of the (k + l)st order Adams-Moul ton formula and the k-th order backward differentiation formula: {EnrightF (k+1) } -> {Af1F (k+l) } + k y* h J{BDF (k) }. Mere J is a Jacobian approximation and y* is a negative constant defined K in [11, d. 195] and [10, o. 11 3]. Py adjusting the coe f ficient k y*, it is possible to widen the stability region without sacrificing accuracv. f Blended formulas modified in this way have been put into the CACM alaorithm DIFSUB (Gear [9]) and compared with the BDF's (see §6). The implementation of these formulas required the derivation of the flordsieck representation (§3) and the development of local error estimators (§5). There is a difficulty associated with thp use of either mul ti deri vati ve or variable matrix coefficient multistep formulas: the use of modified Newton iterations to solve the implicit equations requires "inverting" polynomials in the Jacobian matrix. This would be undesirable if sparse matrix techniques are being used, because matrix mul tinl i cations can considerably reduce sparseness. Two ways of avoiding matrix multipli- cations have been suggested in the literature for the case m = ?., "illoughby t A term suggested by C. W. Gear [6, p. 102] suggests doing the LU decomposition in conplex arithmetic. Enright [6] suggests choosing the parameters of the formula so that the quadratic in h J becomes a perfect square; the computation reouired is one LU decomposition (in real arithmetic) and Uio backsolvos. A third alternative is suggested in §4, which also requires one decomoosi tion and two backsolves. These last two methods of avoiding matrix multipli- cations also avoid the multiplication of a vector by a matrix, which is (k) suggested by the appearance of J{BDF VIW } in the AMF-BDF blend. 2. DERIVATION OF BLENDED FORMULAS We discuss two ways of converting a linear n-deri vative multisteo formula into a linear combination of m linear multisteo formulas: (i) mixed order approach) (ii) pure order approach. For each approach an example is presented followed by a brief discussion of the general case. Consider Enright's third order formula: 2h , h f h 2 fl y -y -i =— T"f + o-f -i - r— f • 7 n ■'n-l 3 n 3 n-1 6 n The associated linear operator, , 2h ',^ . h ..,/. ^ h 2 L n y(t) = -y(t) + y(t - h) + ^y'(t) + ^y'(t - h) - £-y"(t) , satisfies L, y(t) = 0(h ). The main idea is to snlit L. into a sum L. + h -rr Lr where L, and L,' are 1-derivative operators (that is, the second derivative of y(t) is not present). Me choose L, so that lI y(t) = 0(h 4 ), which then forces L 2 y(t) = 0(h 3 ). Th^ conditions on L, imply that lJ y(t) = -y(t) + y(t - h) + h 3 y'(t) + h b ] y'(t - h) + ... + h e k y*(t - k h). where the e's are chosen so that L v(t) = 0(h ). There is a unique choice which minimizes the stepn umber k: Ljy(t) = -y(t) + y(t - h) + f|y'(t) + ^|y'(t - h) - ^-v'(t - 2h), which corresponds to the third order Adams-foul ton formula. This choice of L, forces h Lj; y(t) = - g-{- |y(t) + 2y(t - h) - l-y(t - 2h) + h y'(t)} . The expression in braces corresponds to the second order backward differen- tiation formula. Having obtained our decomposition L. + h tt L. , we 1 2 construct the blended formula by replacing L, and L, with their corresponding formulas and by replacing -rr wi th the Jacobian matrix J = f v (t«, y«): In an actual implementation J would be updated whenever f was reevaluated. Two important properties of this blended formula should be noted. First, the blended formula is identical to Enright's formula for the test equation y' = Ay, and hence the region of absolute stability R = {hx: y -> as n •> « for y* = Ay} is the same for both formulas. Second, the normalized truncation error is for some t and x between t - 2h and t, and so the formula is of order three. Moreover, the terms in braces tend to cancel each other, particularly for the equation y 1 = Jy for which the truncation error would be [1 - hJ]-l { ^! y (4)( t ) + 0(h 5 )}. Here 1 denotes the s x s identity matrix. This convention of using scalars to represent scalar multiples of the s * s identity matrix is used throughout the paper. Having presented an example, let us briefly look at the general case where k m . k f .v L. y(t) = - i a. y(t - ih) + z h J i g.. y u, {t - ih) . n i=0 1 j=l i=0 J1 Assume the order p 5 m - 2, and let q = max{k, p - 1}. We write 1 . .1 . l d ,Z , , , m-1 d . m L h - L h + h dt L h + ••• + h ^n" L h where Ljj y(t) = -a 3 Q y(t) - ... - aj| y(t - qh) + h Bfj y'(t) + ... + h s;J y'(t - qh) "I O ml This decomposition is uniquely specified by choosing L, L, , ..., L ~ (in that order) so that and ii = , p - j < i $ q il y(t) = o(h p+2 ~ j ) . 10 m Having obtained L, L, ..., L. , we define the mixed order blended formula to be m . , q . ' . Z (hJ) J_l E (-a] y + h ^ f •) = . (2.1) j=l i=0 n n_1 1 n_1 In the second approach all the formulas in the blend are of the same order. Again we illustrate by means of Enright's third order formula. 1 A 1 O The operator L. is split into a sum L, + h -rr L where L, and L^ are 1 -derivative operators. Except this time the decomposition is uniquely specified by requiring L? to correspond to a third order formula of minimum stepnurnber: L h y(t) = - B { ir y(t) + 3y(t " h) ■ i y(t - 2h) + g-^* - 3h > + hy, (t)}. The expression in braces corresponds to the third order backward differentiation formula. This forces Lj y(t) - - y(t) + y(t - h) + If hy'(t) + |hy'(t - h) - §-y'(t - 2h) + i^y'(t - 3h) , which corresponds to a linear combination of the third order Adams -Bash forth 2 (3) 13 (3) and Adams-Moulton formulas, namely y=- ABF V ' + jr AMP '. The pure order blended formula is defined to be {- y + y n + ^hf +|hf ,4f o + i4f o> -'n J n-^ 36 n 6 n-1 4 n-2 18 n-3 -^(4^-y + 3y ,-|y o + i-y -, + hf } = 0. 6 6 J n •'n-l 2 J n-2 3 n-3 n Again this has the same absolute stability region as the third order Enright formula. Note, however, that this is a 3-step formula whereas the mixed order approach yields only a 2-step formula. In the general case where L, is a p-th order m-derivative k-step operator, we write i .J il d,2 , . m-1 d , m L h " L h + h dt L h + • • • + h TmT L h dt where L J h y(t) = - c^ y(t) - ... - J q y(t - qh) + h 0fj y'(t) + ... + h 6^ y'(t - qh) and q = max{k, p}. The decomposition is unique if L , L ",..., L are chosen (in that order) so that a^ = , p < i * q , and L J h y(t) = 0(h p+1 ) . TO m Having obtained L, L, , ..., L, we define the associated pure order blended formula as in (2.1 ). The tnright formulas are attractive because of their ideal stability behavior near hx = and near hx = °°, and so it would be of interest to examine the corresponding blended formulas in the general case. Enright's (k + l)-st order formula k-1 ? y n "Vl =h ^/li Vi +h 6 20 f n corresponds to the mixed order blended formula given by * k where {AMF^ k+1 h + k y* hJ{BDF^ k h AMF^ K u : y-y 1 =hZ3-f-, J n J n-\ . q l n-i ■ BDF^ M : Z a. y . = h f . i=0 ' n_1 n This relationship between the three formulas is easily verified once we know that 6 k = (-l) k y£ and a k = (-l) k_1 /k. Note that the order of the AT4F^ k+1 ^ - BDF^ blend is not affected by changing the coefficient k y.*j and hence we consider the formula {AMF (k+1) } - y^ hJ{BDF (k) } (k) (k) where y is a free parameter. In particular y could be chosen to maximize the angle a for which the formula is A(a)-stable. (A formula is A(a)-s table if e^jery numerical solution of y' = xy tends to zero for |Arg(-x)| < a.) Values for y numerically determined in this way appear (k) (k) in Table I in. the column y v ' = Yq D 1. The largest angle a for which ea ch 10 formula is A(a)-stable is also given. A(a)-stable formulas up to order 12 were obtained. The second, third, and fourth order formulas are A-stable for a range of values of y » y » and y . The actual choice of these three parameters is delayed until §4. (k) Table I. Values for y v ' and Corresponding a-angles a for a for Stepnumber k Order k + 1 - k Y k (k) Y opt (k) l * y v ' = - k Y k (k) .. (k) Y " Y opt 1 2 .5 [ 0,+-) 90° 90° 2 3 .1666667 [.1250000,+°°) 90° 90° 3 4 .125 [.1218908, .6837917] 90° 90° 4 5 .1055556 .1284997 86.5° 89.4° 5 6 .09375 .1087264 79.6° 87.0° 6 7 .08561508 .09625961 72.5° 82.9° 7 8 .07957176 .08754865 61.5° 77.4° 8 9 .07485229 .08105623 40.2° 70.2° 9 10 .07103299 .07599874 not A(0)-s table 60.7° 10 11 .06785850 .07192936 not A(0)-stable 47.6° 11 12 .06516462 .06857227 not A(0)-stable 28.7° The second last column was obtained from measurements of Figure 3 in [4, p. 21]. The angle a is only one of a number of parameters which have been proposed for measuring the extent of the stability region R. But it is 11 probably the best such measure, especially for methods with automatic stepsize selection. When such methods are applied to the equation y' = Ay + g(t) with A complex, the integration starts out with hA near the origin, and as the integration proceeds, the stepsize h increases and the value hA moves away from the origin. However if hA approaches the boundary of R, this would be detected by the error estimator, and any further move- ment of hA would be prevented. The implication of this is that not all of R is "used", but only that portion which can be reached from the origin by rays lying entirely inside R. The useful part of R can be described mathematically as the set {y e R: t p e R for < t ^ 1}, and for stiffly stable formulas this region is much like a wedge, which is best described by the angle a. This parameter also serves as a good indicator of the problem class for which the formula is suitable, namely those problems for which the "large" eigenvalues A of the Jacobian satisfy | Arg(-x) | < a. Ideally we would like to have o = 5- so that the formula would be suitable for all realistic problems. It may be helpful to present the original motivation for the (k+1) (k) AMF V - BDF V blend. We begin by considering a scalar equation y 1 = f(t, y). If the stepsize h were determined only by accuracy considerations, then the size of -h f would be a good measure of stiffness. If -h f were small, then the AMP + '' would be a good k-step formula to use, and if -h f were large then the BDF V ; would be a good choice. And so it would seem that a weighted sum {AMF (k+1) } + y (k) (-h f ){BDF (k) } 12 would give us the best of both formulas. For small -h f this blend is y nearly the AMP k+1 \ and for large -h f it is nearly the BDP k ^. We now extend this idea to a system. Consider a change of coordinates which decouples the system at (t, y) = (t , y ). Such a coordinate transformation can be accomplished by a matrix T which diagonalizes f (t , y ); that is, f" f v (t , Y n )T= A where a is diagonal. This is almost always possible. Letting z = T~ y, the system becomes z' = g(t, z) = T" 1 f(t, Tz). Since this system is locally decoupled, we can use -h a 11 as a measure of the stiffness of the i-th differential equation. Each equation of the system can be integrated by a formula which is suitable for that equation: Transforming this back to the original coordinate system gives K + Vl + h J 6 J W + T {k) (-h V ( "j 'i y n-.i + h f n> " °- which is the blended formula. The remainder of the paper is restricted to a discussion of the AMF^ k+1 ^ - BDF^ blend with Y ^ = y^i- 'opt 13 3. THE NORDSIECK REPRESENTATION In the Nordsieck representation we work with an approximation ^ = [y n > h y^, ..., h k y^/k!] T to the scaled derivatives [y(t n ), h y'(t ), . .., h k y^ k ^(t n )/k!] T of the solution. A linear Nordsieck formula can be thought of as consisting of a prediction ^ = P Vl followed by a correction *n = En + *■ A n where a is chosen so that y 1 = f(t , v.,); that is, a is defined implicitly n •'n n n n v J by A n = ju^h f(t , p n + £ n a J - h p'} Here p^ = [p n> h p^, . . . , h k p^/k!] 1 , £ = [£ Q , a } £ k ] T , and P is the Pascal triangle matrix given by P id = ( ] ) , « i, j $ k. Nordsieck formulas are closely related to multistep formulas. For every linear k-step formula of order >, k there is a unique (k + l)-value linear Nordsieck formula such that the values y , y -,,..., y ■ generated by the Nordsieck formula from v . satisfy the difference equation of the multistep formula. Let the multistep formula be p(E) y , = c(E) h f . x . J n-k n-k where E is the shift operator, p(0 = Z a. £ , and aU) = E B 5 k_1 . 1=0 ] i=0 ^ The corrector vector _£ of the "equivalent" Nordsieck formula satisfies 14 p(c) = u - D k+1 &f(5 I - P)' 1 L (3-D and oU) = (E - D k+1 ej(e I - P)" 1 * (3.2) where e, = [0, 1 , 0, . . . , 0] and e« = [1, 0, 0, ...,0]. It can be shown from these equations that l- = a Q and z Q = 6 . (3.3) Proofs of these relationships will appear in [19]. (k) (k) Let a_ v ' and tr ' denote the corrector vector corresponding to the k-step AMP k+1 ) and BDP k ' respectively. It can be shown [19] that a (k) = . k ; ] (.dj.iw) and b (k) u .lr k+1 1 • m Q(1)k D j k! L j+1 J ' J UU ' K ' where the brackets denote Stirling numbers of the first kind (see e.g. Knuth [13, p. 66]). Values of a^^ and b/ k ^ are given in Appendix A for k = 1(1)11. Consider the Nordsieck formula with £ = a • y h J b , where we have suppressed the superscripts (k). Applying (3.1) and (3.2) gives p(0 = p AM U) - y h J p BD U) and aU) = a AM (?) - y h J a BD (s). Therefore the values computed by this Nordsieck formula satisfy P (E) y n _ k -ha (E) f p _ k - y h J{p (E) y p _ k -ha (E) f p _ k > = , (k+1) (k) which is the difference equation for an AMF V ' - BDF X ' blend. 15 In an actual implementation J might be reevaluated every few steps; in other words, where i = a - y h J b. And a blended Nordsieck formula with a changing -n — ' n — 3 3 J is not equivalent to a blended multistep formula with a changing J. For example, the blended 3-value Nordsieck formula generates values y , y, , y _p which satisfy P V2 + e^(I -!„*[> P Vl hj, n-l + §0^ hy n ' 16 which is a system of two equations relating y. y , , y „ hy', h y 1 . h y;_ 2 , and h 2 y^/2. yields equation (3.5). 2 2 h y' os and h y" /2. Eliminating h y" /2 and substituting f for y 1 n-d n-£ n-£ n n 17 4. THE CORRECTOR ITERATION The implicit equation satisfied by a is h p 1 + a, A - h f(t, p + & a) = where for convenience we suppress dependence on n. The Newton iteration for this equation is A (0) = °» A (m + 1) = A (m) + ^1 - h ( f y )( m ) ^ ]_1{h f (m) " (h P' + *1 A (m) )} where for a function the subscript (m) denotes evaluation at (t, p + £ Q A/ v). This iteration is the same for both the multistep form and the Nordsieck form of the formula. Convergence of this iteration is seldom a problem, and so a considerable saving in time results if a recent Jacobian approximation J is used in place of ( f y )( m ) : A (m+1) = A (m) *CV hJ V" 1{h f (m) " (h P' + £ 1 A (m) )} ' < 4 '« From i = a -y h J b and from (3.3) we have that £, - h J i Q = (1 - y a h J) - h J(b - y h J) BD AM where a = a n and 3 = $ Q . This means that (4.1) involves the LU decomposi- tion of a quadratic polynomial in h J. But forming matrix products can be expensive, and this is particularly true if sparse matrix techniques are being used because matrix multiplications reduce sparseness. One way of avoiding matrix multiplications is suggested by Willoughby [6 , p. 102], and it is applicable as long as the roots of 1 - (y a + e)z + y z are complex conjugates. This technique calls for an LU decomposition in complex arithmetic, which has the disadvantage of being four times slower and requiring twice as much storage compared to real arithmetic. 18 A second possibility is to choose the parameter y so that 2 2 1 - (y a + 3)2 + y z is a perfect square (1 - c z) . Then multiplication -2 by (1 - c h J) can be accomplished by decomposing 1 - c h J and performing two backsolves. Enright [6] derives a second set of formulas based on this idea. This approach, however, reduces the number of free parameters, and so one order of accuracy must be sacrificed. There is a third way of avoiding matrix multiplications, and that 2 2 is to approximate 1 - (y a + b)z + y z by a perfect square (1 - c z) so that instead of (4.1) we use | . Clearly we must choose c > 0, and so the expression inside the vertical bars is analytic for Re y $ 0. By the maximum modulus theorem the maximum is attained for Re y = 0, and so max | (t «+B - 2c)1 M - (c 2 - T ) M 2 | -°° b - 2c) 2 + (c 2 -tfJl \_\ + C w J Differentiation with respect to u yields extrema at J = 0, 2 iii and 2 OJ = (y a + 3 - 2c) c 2 (y a + 6 - 2c) 2 - 2(c 2 - y) 2 if 2(c 2 - Y ) 2 < c 2 (y a + 6 - 2c) 2 ; therefore V - ^ if 2(c 2 - Y ) 2 * c 2 ( Y a + 6 - 2c) 2 , Mc)f - ( max { (c 2 - v) 2 v« (y a + 3 - 2c)' 4[c 2 (y a + 6 - 2c) 2 - (c 2 - y) 2 ] } otherwise The value of c which minimizes 41(c) was determined numerically for k = 1(1)11. These values appear in Table II. 21 Table II. Values for c (k) Stepnumber k (k) *(c< k )) 1 .2928932 2 .3374973 .1184662 3 .3335427 .1161818 4 .3427329 .1139836 5 .3169058 .0995555 6 .2992971 .0894459 7 .2862392 .0819169 8 .2760327 .0760633 9 .2677630 .0713660 10 .2608834 .0675036 11 .2550426 .0642651 From Table II we see that even in the worst case, x pure imaginary and k = 2, the asymptotic error coefficient is only .118, which means that each iteration is good for an additional 0.93 digits. 22 5. ERROR ESTIMATION AND CONTROL The purpose of estimating the local error is to assist in the selection of stepsize and order in such a way that (i) the least possible work is done for the accuracy attained. (ii) the global error behaves in a reasonable manner as a function of the local error tolerance. Hence the estimation of the local error is not in itself of much importance The theory of error estimation is not well developed. Usually an attempt is made to devise a local error estimate which is asymptotically equal to the local error under the assumption that the order and stepsize do not vary. This approach must be used cautiously. For example the (strongly A-stable) formula y = y , + .501 h f ., + .499 h f has the v 3 J ' J r\ -'n-l n+1 n asymptotically correct error estimate .001 h(f ., - f ). However, this will often give a \iery poor estimate, simply because the leading term in the asymptotic expansion of the local error is not the dominant source of error unless h is \ievy small. A more realistic approach would be to make use of a good bound on the truncation error. This approach would yield the estimate .250001 hi If - f ,1 for our example. It is worth noting 1 ' n n- 1 ' ' r that the highly rated code DE of Shampine and Gordon [17] uses error estimates that are not asymptotically correct; for example, the estimate .5 h(f - f , ) is used for the trapezoid formula. v n n-1 K Our implementation of the blended formulas uses an error bound as the basis for the error estimate. In order to derive the error bound, we first split the AMF^ k+1 ^ into a sum of the AMF^ and the ABF^: AMF< k+1 > =^AMF< k > + ^ULA ABF ( k ) Y k-1 Y k-1 23 where y. is given in Henrici [11, p. 193] and Gear [10, p. 108]. This equality follows from the fact that e Q = y, for the AMP '. From Henrici [11, equation (5-22)] we have that y. - y*-. = y* , and hence Y k-i Y k-i If we consider the AMF-BDF blend applied to a scalar equation, we get that the (normalized) truncation error is h w (1 . Y h x) -igJi y (t) . IkJl y (W)( T ., . rU, * f(t, y) '""cH «- c x h if desirable then / W <- 1 - cH x f (t, y) ^decompose W solve W x d = h x f - hy' for d do_ / solve W x d = d for d d^ Y x h>< (d- d)/cH y *- y + a n x d + d hy' ^ hy' + d + b 1 x d A^A + d;A«-A + d rf d + d small enough then go to testerror deal with iteration convergence failure testerror: rf (a + A)/(k +1) too large then go to tryagain for J + 2(1 )k do hV J "Vj! «■ hV AJ1 + a. * A + b. x a In both cases one function evaluation is required for each iteration. However the corrector overhead for the AMF-BDF blend is twice as great as for the BDF. The overhead is defined to be any computation other than function evaluations. A significant portion of the overhead of a stiff system solver often consists of LU decompositions. We would not expect there to be a difference between the two formulas in the number of decompositions if we were to use the same type of iteration in both cases. For the BDF we use A,„.,x = A, _, + [£, - H J £ n ] _1 {...} V." (m+1) ~ "(m) 1 '0- 26 where H is the size of the step when J was last evaluated. The same type of iteration for the blended formula would be A (m+1) = A (m) + [(1 " y « h J) - H J(g - Y h J)]" 1 !...} However what we actually use is Vl) =A (m) + [1 -cH J]- 2 (...}. This second iteration might be expected to fail more often than the first, but it is difficult to analytically compare the two iterations even for the equation y' = Ay, and therefore we have to rely mainly on empirical evidence. If the preceding discussion is generalized to blends of m formulas, it can be seen that the correction iteration overhead is proportional to m. Hence blends of m formulas are competitive only if they require significantly fewer function evaluations to achieve a given accuracy. How many fewer depends on the ratio of the cost of a function evaluation to the cost of a backsolve. For problems with very simple right hand sides f(t, y) , the cost of a function evaluation might be roughly equal to the cost of a backsolve and in this case it would not pay to use a blend of m formulas unless the number of function evaluations compared to that of the BDF's is smaller by a factor of 2/(m +1). Here we are assuming that the number of decompositions is independent of m, which is fairly reasonable since it is primarily the heuristics used by the program that determine the number of decompositions. We have compared the efficiency of the blended formulas to that of the BDF's on several test problems. In order to make the comparison fair, the blended formulas were implanted into a well-known BDF code (Gear's DIFSUB [9]) in such a way that only the formula-dependent part of the code 27 was changed. And therefore the code was not tuned for the blended formulas. There was one modification made to the codes and that was to replace the call to MATINV by calls to DECOMP and SOLVE [8]. A listing of the source ft code appears in Appendix B. In order to facilitate testing, only problems with known analytic solutions were used. The stepsize h was initially set to the local error tolerance e, and on subsequent steps it was set to min{h, t f - t} where h is the stepsize recommended by the method and [0, t f ] is the interval of integration. The parameter h was set to t f - t and h . was (mistakenly) 3 max t mm 1 -13 set to j u t f where u is the unit roundoff error, 16 , of the IBM 360. The weight vector w was initialized to all ones. The Jacobian f (t, y) was evaluated numerically.. The relative (global) error was defined by i i / \ max / J , y n " y ^n\ 2xl/2 - ' OsnsN { . , * i ; ; i=l where w = max{l, y n , y, , ..., y }. Tables III through VII give the results of comparing DIFSUB with blended DIFSUB for five different types of problems. Each problem was -2 -3 -10 integrated for nine different error tolerances, e = 10 ,10 , ..., 10 For each integration the following statistics were collected: max order - the order of the highest order formula selected by the code. steps - the total number of steps taken. func evals - the total number of function calls (including those needed to approximate the Jacobian). back solves - the total number of calls to the subroutine SOLVE, which solves a linear system that is in factored form. LU decomps - the total number of calls to the subroutine DECOMP, which performs an LU factorization on a matrix. 28 accurate digits - the negative of the base 10 logarithm of the relative (global) error, which was defined previously. time - machine time in centiseconds on the IBM 360 model 75 at the University of Illinois at Urbana-Champaign. These timings are not completely reliable due to multiprogramming. It should be stressed that the intent is to compare two classes of formulas and not two methods or two programs. Also, only one aspect of the relative utility of the two formulas is being compared, namely machine time versus global accuracy. Problem I. The following is a linear problem whose Jacobian matrix has eigenvalues -.1, -50, -120: -.1 -49.9 -50 70 -120 y(o) ■ This was integrated on [0, 15] 29 Table III. Numerical Results for Problem 1 max steps func back LU accurate time £ order evals solves decomps digits DIFSUB ID" 2 3 29 104 76 9 1.9 30 _3 10 5 49 145 108 12 2.9 38 io- 4 4 70 202 171 10 3.9 51 10" 5 6 106 286 240 15 5.0 75 10" 6 5 193 508 453 18 5.7 128 IO" 7 6 176 474 422 17 6.7 128 10" 8 6 204 551 505 15 7.5 143 IO" 9 6 270 771 719 17 8.4 185 IO" 10 6 353 1024 969 18 9.3 251 bier ided DIFSUB IO" 2 4 30 101 146 9 3.2 33 10 J 5 43 141 220 10 3.9 45 IO" 4 6 69 195 316 12 4.6 71 IO" 5 7 89 238 396 13 5.4 86 IO" 6 8 120 308 524 15 6.5 no IO" 7 11 139 410 662 26 7.4 155 IO" 8 9 169 443 782 17 8.5 165 IO" 9 11 209 540 952 21 9.3 228 IO" 10 11 234 595 1056 22 10.4 261 An examination of the last two columns reveals that there is a slight though definite improvement in performance when the blended 30 formulas are substituted for the BDF's. The reason that blended DIFSUB does more LU decompositions might be due to the fact that fresh decomposi- tions are performed whenever the order is changed but not when the stepsize is changed. Problem II, The following is problem 12 in Krogh's [14] set of test problems; it is also used by Gear [10, p. 218]: - 800z^ + {z ? y' = u -lOOOz 1 + (z 1 ' .2 . ,_2 where z = Uy and 10z 3 + (z 3 4 / 4 -.001z H + (z* -1 U= I 1 y(o) = -l -l -l l -l l -i This was integrated over [0,1000]. The eigenvalues of the Jacobian matrix are equal to -1002, -802, 8, and -2.001 at t = 0, and approach -1000, -800, -10, and -.001 respectively as t ■*■ ». 31 Table IV. Numerical Results for Problem 2 e max order steps func evals back solves LU decomps accurate digits time DIFSUB 10 " 2 3 64 277 184 23 1.7 65 10 " 3 4 105 374 273 25 2.7 91 ID' 4 5 167 496 387 27 3.6 131 lO" 5 6 250 733 608 31 4.1 206 10" 6 6 284 778 677 25 5.4 229 lO" 7 6 380 1076 915 40 6.4 318 10" 8 6 467 1311 1178 33 7.0 416 lO" 9 6 590 1641 1516 31 8.2 586 lO" 10 6 768 2083 1954 32 9.1 699 bier ided DIFSUB in" 2 4 61 276 358 24 2.7 75 lO" 3 6 98 333 488 22 3.7 105 lO" 4 7 146 449 688 26 4.6 160 10 D 7 200 594 970 27 5.4 231 10" 6 9 262 750 1242 32 6.4 314 lO" 7 9 325 900 1502 37 7.3 419 10 -8 9 381 1062 1818 38 8.4 521 lO" 9 11 474 1303 2244 45 9.4 734 10" 10 11 558 1548 2678 52 10.4 852 Once again there is a slight though definite improvement when the blended formulas are substituted for the BDF's. 32 y. = y(o) = Problem ill. This is problem B5 in the set of stiff problems considered by Enright et al [7]: -.1 This was integrated over the interval [0, 20]. Because the Jacobian matrix has eigenvalues -10 ± ilOO, -4, -1, -.5, -.1, this problem can be troublesome for the BDF's. - 10 100 -100 -10 -4 -1 -. 33 Table V. Numerical Results for Problem 3 £ max order steps func evals back solves LU decomps accurate digits time DIFSUB io- 2 (4) (1001) (3002) (2959) (7) (1.1) "too much work' ' ; integration s topped at t = 8.3 blende d DIFSUB ID" 2 6 125 493 756 19 2.9 215 10" 3 6 194 691 1152 19 3.7 348 IO" 4 7 277 922 1590 21 4.5 517 10 3 8 368 1196 2114 23 5.4 729 IO' 6 9 468 1494 2686 25 6.4 1040 IO" 7 11 592 1831 3288 31 7.3 1546 IO" 8 11 721 2178 3958 33 8.5 1772 IO" 9 12 895 2644 4878 34 9.4 2369 IO" 10 (12) (1001) (2917) (5580) (21) (10.3) "too much work' ' ; integ ration s topped at t = 5.1 There is a remarkable difference between the abilities of the two types of formulas to solve this problem. Problem IV. The following is problem 13 in Krogh's set of problems; it is also used by Enright et al [7, problem E4]: 34 lOz 1 + 10z 2 + \ (z 1 ) 2 - \ (z 2 ) 2 1 2 12 -lOz 1 + lOz + z 1 z -lOOOz 3 + (z 3 ) 2 L -01z 4 + (z 4 ) 2 where z = Uy and -1 1 1 1 -1 1 1 1 -1 1 1 1 -1 y(o) = o -2 -1 -1 u = i 1 This was integrated over [0,1000]. The eigenvalues of the Jacobian matrix are equal to 8 ± lOi, -1002, and -2.01 at t = and approach -10 ± lOi, -100, and -.01 as t -> ». 35 Table VI. Numerical Results for Problem 4 E max order steps func evals back solves LU decomps accurate digits time DIFSIJB io- 2 3 79 388 263 31 1.3 78 10- 3 4 137 555 402 38 2.0 116 IO" 4 5 289 836 715 30 3.3 249 IO" 5 6 268 802 677 31 4.2 229 IO" 6 (6) (1001) (2777) (2652) (31) (5.0) "too much work' " ; integration s topped at t = 150 blende d DIFSUB IO" 2 4 79 347 484 26 2.6 95 IO" 3 5 127 463 708 27 3.5 140 IO" 4 7 170 558 890 28 4.7 191 -5 10 ° 7 244 750 1258 30 5.4 266 IO" 6 7 325 952 1646 32 6.4 376 IO" 7 9 399 1165 2008 40 7.4 536 IO" 8 11 483 1367 2364 46 8.3 724 IO" 9 (11) (1001) (5511) (6708) (539) (9.2) 1389 "too much work 1 1 ; integ ration s topped at t = 977 For this problem the blended formulas performed significantly better. Problem v. The following problem (problem 8 of Krogh's set of test problems) tests the integrator's ability to solve nonstiff problems: 36 y ■ 1,, K2 , 2v2v-3/2 -y ((y ) + (y ) ) 2,, K2 / 2v2 N -3/2 -y {(y ) + (y ) ) This was integrated over [0, 20]. y(o) = 37 Table VII. Numerical Results for Problem 5 e max order steps func evals back solves LU decomps accurate digits time DIFSUB ID" 2 4 50 323 210 28 -0.4 73 ID" 3 5 81 345 264 20 0.4 91 IO" 4 5 95 453 324 32 0.8 115 ID" 5 5 140 424 403 5 1.6 130 ID" 6 6 161 487 462 6 3.1 158 ID" 7 6 233 681 656 6 4.2 224 ID" 8 6 294 865 840 6 5.1 283 _Q 10 3 6 401 1174 1149 6 6.1 383 io- 10 6 588 1662 1637 6 7.3 534 ble nded DIFSUB in' 2 5 59 278 410 18 1.4 100 IO" 3 7 63 212 366 7 2.8 100 IO" 4 8 77 247 436 7 3.7 113 10" 5 8 99 307 556 7 4.5 145 IO" 6 10 106 335 596 9 5.1 165 IO" 7 10 128 394 714 9 6.0 228 IO" 8 10 157 470 866 9 6.9 273 IO" 9 11 170 505 928 10 8.2 299 IO" 10 11 200 579 1076 10 9.5 361 Again there is a significant difference in performance. 38 7. IMPLICIT DIFFERENTIAL EQUATIONS One of the advantages of the Nordsieck representation is that a predicted value p' for y' is available, and hence implicit differential equations F(t, y, y') =0 can be treated without additional complication (cf. Gear [10]). We simply chose A to satisfy F a< 2 > a< 3 > «W a< 5 > a< 7 > 1 2 5 12 3 8 251 720 95 288 19087 60480 1 1 1 1 1 1 1 2 3 4 11 12 25 24 137 120 1 6 1 3 35 72 5 8 1 24 5 48 1 120 17 96 1 40 1 720 42 a'< 8 > a< 9 > ,(10) a( ]1 > ,(12) 5257 17280 1070017 3628800 25713 89600 26842253 95800320 4777223 17418240 1 1 1 1 1 49 40 363 280 761 560 7129 5040 7381 5040 203 270 469 540 29531 30240 6515 6048 177133 151200 49 192 967 2880 267 640 4523 9072 84095 145152 7 144 7 90 1069 9600 19 128 341693 1814400 7 1440 23 2160 3 160 3013 103680 8591 207360 1 5040 1 1260 13 6720 5 1344 7513 1209600 1 1 8960 29 96768 121 40320 193536 1 1 / 72576 11 362880 272160 1 11 3628800 7257600 1 39916800 43 (1) (2) (3) (4) (5) (6) 1 1 1 1 1 1 1 3 2 11 6 25 12 137 60 49 20 1 2 1 35 24 15 8 203 90 1 6 5 T2 17 24 49 48 1 1 35 24 8 1 120 144 7 240 1 720 1 363 11 1 761 44 ,(9) 1 7129 140 280 2520 469 29531 6515 180 10080 2016 967 267 4523 720 160 2268 7 1069 95 18 1920 128 23 360 9 80 3013 17280 1 13 5 180 960 192 1 1 29 5040 1120 12096 1 1 40320 8064 1 362880 (10) 7381 2520 177133 50400 84095 36288 341693 362880 8591 34560 7513 172800 121 24192 11 30240 n 725760 1 3628800 (ID 83711 27720 190553 50400 341747 129600 139381 120960 242537 725760 1903 28800 10831 1209600 n 13440 1 20736 1_ 604800 1_ 39916800 45 APPENDIX B. SOURCE CODE FOR BLENDED DIFSUB IMPLICIT REAL*8( A-H.O-Z) CCMMOK /UmTGER^L IPRPIF .NFfN.NC TNV.NSQLVE CCMMON /INMAIN/ ICUIT CALL UNDERZI 'OFF') £E*D ( 5.800. ENC=2C ) I F FOB . IPO i£R *RITE(6.lOCC) IPROB 1000 FOPMAT(/////'l», t FRRnP Xl«nO .0.00.3*0.00 / 2CDO,1G00.00.15.DO,15.DC20.D0.2C.D0.10CO.D0.3*O.DO 1 .00. 15*0 .DO* - 1.00*^1 .00.-1. 00 ♦- 1^00^16 *0^D0» 1 COO. 19*0. DO. 2 .DO. 1 .00 .2. CO .17*0. DO. 1 .QO.C.OO.C.nC.l .00.16*0.00. **» .nn.i A*n T nn. CATA CATA TF YO 60*0.00 / ! 0.D0.-2.O0.-1 .00.-1 .00 .16*0. DO. CATA NOIM / 1.4,1,3.4.6.4.3*0 / CATA FCURU / Z334 COOC OOOOOOCO / CATA MAXIMUM / 1CCC • NFCN = *AXNC = NCSTEP - NSCLVE = KOINV = ERFMAX = O.CO T MXERR = COO NECN = NDIM(IPROE) __ _ T = TO(IPRCE) ITIME = CC 10 I = 1 .NECK Y( 1 . [ ) = Y0< I . IPRCB) YMAX(I) = DMAX1 { I .DC.CABSi Y<1 .1 ) ) ) F.MIN = F0URU*QMAX1(0AFS< T) .CABS F = TCL JSTART = IF { T.GE.( TF( IPRCB) - HMIN) ) GO TO hMAX = TF(IPROB) - T F = DMAX1 ( FM IN.DMIN1 (HMAX.H) ) NCSTEP = NOSTEP -k I IF ( NCSTEP. LE.MAXNU* ) GO TC 35 IQUIT = 1 GO TO 50 35 CALL CIFSUB(NEQN.T. Y .SA VE.H .HM1 N»FWAX.TC4_^*V*1AX. ERROR. 1 K FLAG. JSTART. 11, Pte, IP .OMEGA .OMEGA 2) CALL CALERRIT^JY-^VtiAX.GLEEBR ) IF ( GLBEWR.LE. ERRMAX ) GO TO 40 ERRMAX - GLBERfi T MXERR = T JSTART = 1 IF ( KFLAG.EQ.l ) 10 30 70 40 GO TO 30 5-0 *AITE(6.9Q9) TQL . NOST EP . NFCN . NO INV . NSQLVE . M AXNQ . EflCM AX i T MX EBP. 1 ITIME.KFLAG.T 999 FORMAT ( »0 • .3X.D9.3.2X.IS.2X.XC 3X.I4.SX.I4.7X.12. 1 8X.D11.5.7X.D11.S.5X.I6.4X.I2.4X.010.4) RETURN 70 IftRITE (6.999 ) TOL . NOSTEP . NFCN . NO INV. NSOL VE .MAXNQ .ERRMAX. T MXERR. 1 I T I * E RETURN END ____ _. 47 SUBROUTINE CECCMP ( P* • h • M « I P ) T REAL*8< A-H.Q-Z) RE^L+4 PW CIMENSICN PI(M)«IP(»i) CCMMON /INTGER/ I PROE .NFCN .NO I NV .NSOLVE NCINV = NOINV + 1 IP(N) = 1 Cd 6C.K.s__l^jj IF ( K.EC.N ) GC TC 5C KP1 = K ♦ 1 - KM1M = (K - 1 ) *M IPIVCT = K CO 10 I = KP1.N IF ( ARSCPJUC I»KM1 W) ) .fiT . ARM P*< I P I VHT+KMl M) ) ) IPLVOT s._l_ 10 CONTINUE IP(K) = IFIVOT IF ( IPIVCT. NE.K ) IP(N) = - IP(N) IPKM = IPIVCT ♦ KM1* KKM = K ♦ KMIM TEMP = PW(IPKN) ' PW< IFKM) = PX( KKM) PW(KKM) = TEMP IF ( TEMP.EQ.O.DO ) GC TO 50 CO 20 I = KPl.h IKM = I ♦ KMIM 20 PWIIKM) = - PW < IKK|/TFMP CO 40 J = KP1,N J VIM = (J - 1 )*M IPJM = IPIVOT ♦ JM1M KJM = K ♦ JM1M . TEMP - FW( IPJM) Pl( IPJNLl. 5„EmiKJM > Pto(KJM) = TEKP IF ( TEMP. EQ. CDC X GO TH,AO^ CO 30 I = KPl.N IJM = I ♦ JM1M 30 PW(UM) = Ptt(IJM) + P*( I+KM1M) *TEMP 40 CONTINUE 5C IF ( PM(KKM) .EC.O.C ) IP(N) = 60 CCNTINUE RETURN ENC 48 SUBROUTINE SOLVE (PM . N • M« IP , B ) U4F1_ICLT RFAI *fl- Z ) CIMENSION Yd3,2G) ,DY(2C) CCMMON /INTGEay I PROB .J-tFCN. NO XJNV .NSOLME NFCN = NFCN ♦ 1 GO TO (10. 2C, 30. 40, 50. 6C. 70. aC. 90. 100). I PROB 10 DY(1 ) = - Y( 1 ,1 i £E IJJRN 2C ZTEMP = .5DC*(Y<1,1) 4 Y(l,2) ♦ Yd. 3) ♦ Yd, 4)) 21 = ZTEMP - Yd.lJ 22 = 2TEMP - Y( 1 .2) 23 = 2TEMP - Yd. 3) 24 = 2TEMP - Yd.4) 21 = Z1*(Z1 - LOCO. DC ) 22 = 22*(22 - 8CC.DC) 23 = 23*(Z3 ♦ 10. DO) Z4 = 24*(Z4 - .0C1D0) 2TEMP = .5DCMZ1 ♦ Z2 ♦ 23 + Z4 | CY(1) - ZTEMP - 21 CY(2) = ZTEMP - 22 CY<3) = ZTEMP - 23 CY(4) = ZTEMP - Z4 RETURN 30 CY<1) = - 2C0.D0*Y(1. 1 J ♦ 2000. DQ - (1991. DO ♦ 1 199.DC*T)*DE>P<-T) RETURN 40 CY(1) = - .1D0*Y(1,1) - 49.9C0*Y(1.2 ) CY(2) = - 5C.DC*Y< 1 ,2) CY(3) = 70.CO*Y(1,2) - 12O.D0*Y< 1.3) RETURN 50 CV(1 ) = Y( 1 .3) CY(2) = Yd. 4) CY(4) = ( Y( 1, 1 ) *Y( 1 . 1 ) ♦ Y( 1 ,2)*Y(1 .2))**1 .5 CY(3) = - Y( 1 .1 )/DY(4) CY<4) = - Y(l ,2)/DY(4 ) RETURN 60 CY(1) = - lC.D0*Yd.l) ♦ 100.00*YC 1 .2) CY(2) = - 1 CO^J-a-YXl-d . - l0.nO»YC1.2) CY(3) = - 4 ,D0*Y( 1.3) CY(4) = - Y(l. 41 CY(5) = - ,5DC*Y< 1,5) CY(6) = - . lD0*Yd.6) __. RETURN 7C ZTEMP s >5DC*(Y(1.1) ± Yd. 2) + Yd. 3) ♦ Yd. 4)) 21 = 2TEMP - Y( 1 ,1 ) 22 = ZTEMP - Y(i .2) 23 = ZTEMP - Y( 1 ,3) 24 = ZTEMP - Y(l. 41 HI = .5D0*(Z1*ZI - Z2*Z2| _2 _-_21*-22 Tl = Zl Zi = 1C.DC*T1 ♦ 10.OC+22 ♦ _-l 22 = - 1C.D0*T1 ♦ 10.D0*Z2 ♦ *2 23 = 23*123 - 1000.DC1 Z4 = Z4* )**2 ((Y(l*3) - ZTEMP ♦ Z2)/YMAX(3) )**2 ( ( Y ( 1 , 4j - ZTEMP-*.. DSQRT(GLBERP) 74)/YMAX(4) >♦*? GLEERR RETURN 30 GLEERR 1 RETURN 40 Tl = DEXP(- 5Q.CC*TJ DAES((Y(1,1) - 10. CO ♦ (1C.00 4 T)*DEXP(-T) 1C.D0*DEXP(- 2CC.DC*!) )/YMAX(l)l GLeERR = ! ♦ 5C GLBERR RETURN CCS T = DCCSdl SIN T * DSIM T) ((Y(l,l) - DEXP(- .1C0*T) - T1J/YMAX(1))**2 ( )/YMAX(3 ) )**2 DSQRT(GLSERR) CCS SIN SIN CCS T J/YMAXC 1 J J**2 T)/YMAX(2) )**2 T1/YMAXC3) 1**2 T)/YMAX(4) )**2 CLEERR = ( (Y ( 1.1 ) 1 ♦ (< Y( 1 .2) 2 «■ { (Y( 1 ,31 3 + ( I »*2 1 ♦ ( >**2 2 + ( < Y C 1 ,3) - ZTEMP ♦ Z3)/YMAX(3) )«*2 3 ♦ ( C* SET THE COEFFICIENTS THAT DETERMINE THE ORDER ANC THE METHOD C* TYFE. CHECK FCR EXCESSIVE CRDER. — THE LA ST TW O ST A T E M E NTS OF C* THIS SECTICN SET I fcEVAL .GT.O IF P* IS TO BE RE-EVALUATED 55 C* BECAUSE OF THE ORDER CHANGE. AND THEN REPEAT THE INTEGRATION * C* STEP IF IT HAS NOT VET BE EN D O N E (I RET = 1) OR S KI P TO A FINAL - * C* SCALING BEFORE EXIT IF IT FAS EEEN COMPLETED (IRET = 2). * C***** 4 **************************** ****** ************ ******************* 170 IF (VF.EQ.O) GO TO 18C IF (NC.GT.ll) GC TO 190 GO TO ( 221, 222. 222.224, 225. 226, 227,228. 229,230. 231) .NO IBO IF (NC.CT.12J G£ TO ISO GO TO ( 201 .202.203. 204,205,206.207,208,209.2 10.21 1,2121 ,NQ 190 KFLAG = -2 RETURN C******* *************************************************** ************* C* TFE FOLLOWING COEFFICIENTS SHOULD BE DEFINED TO THE MAXIMUM * C* ACCURACY PERMITTEC EY T KE KAChlhE. THEY ARF. list THE J3RDER USED.. * C* * C* -1 * C* -1/2.-1/2 * C* -5/12. -3/4. -1/6 * C* -3/8.-11/12.-1/3.-1/24 * C* -251/720. -25/24. -35/72 ,-5/48 .-1/ 120 * C* -95/288. -137/12C.-5/8, -17/96. -1/40.-1/720 * C* - 1 9C e 7/6 C4 80. -49/40 • -203/2 70 • -4 5/ 192. -7/ 1 44, -7/ 1 440 * -1/504 * C* * C* -1 * C* -2/3.-1/3 • * C* -6/1 1 .-6/11.-1/1 1 * C* -12/25. -7/1C. -1/5,- 1/50 * C* -12C/274. -225/274. -e5/274, -15/274, -1/274 * C* -ieC/441. -58/62, -15/36, -25/252, -3/252. -1/1 764 * C** ************************************************************ ********* 20 1 A ( 1 ) = - .500 CO TO 2 25 202 fi ( 1 ) = - .4 166666666666667D0 A(3) a - .500 GO TO 235 203 A( 1 ) = - .375C0 A( 3) = - .7500 A (4) = - .166666666666666 700 CO TO 235 204 A(l) = - .34861 1 1 11 1 1 1111100 A(3) = - .916666666666666700 A(4) = - .333333333322333300 A(5> = - .04166666666666667DO GC TO 235 . 205 A(l) = - .329861111111111100 A(2) = - 1.C41666666666667DC fl{4) = - .486111111111111100 A(5) = - . 1 C4 1666666666667U0 A(6) = - .0C8332233332333333DC GO TC 235 _ _ .. ._. 206 A(l) = - .2 155919312169212D0 A(3) - - 1 . 14 1666666666667D0 A(4) = - .62500 A(5) = - . 1770832233333333D0- — A(6) = - .C25D0 a(7) = - < oci38££8aaeeaeaaa9co GO TO 235 207 Ml) = - .304224537037027000 A(3) - - 1.225D0 56 4(4) = - .751851651851851900 JU-5-1 -=- - .2552083231223133D0 Ait) = - .0486111111111111100 A(7) = - . C4 66 1 1 1 1 1 1 11 1111 ICO Ai8) = - .000 19841269e41269£400 GC TO 235 208 4(1) = - .2948680004409171D0 4<3).J= - 1 .?9642fi5714?84flfinr. 4(4) = - .868518518518518500 AC 5) = - . 32 5 7£3££a£££8£aaDJI 4(6) = - .0777777777777777700 4(7) = - .0106481481481481500 4(6) = - .000793650752650793700 4(91 = - .QO0024fni5fl7.ini5f 7?no GC TO 235 209 4(1) - - .2869754464265714QD 4(3) = - 1.3589285714285700 4(4) = - .9 76 55 423280 4233X10- 4(5) = - .417187500 4(6) = - T l 1 l 3S A l fiA^fi^ft7nn 4(7) = - .01875DC 4(8) = - .00193452380952380900 4(9) = - .00011 160714285714300 4(10) = - .C0C0C275573192239£59DO C-C TC 235 210 4(1) = - .? 80 189 5 64 4 3 936 700 Ai3) = - 1.4144841269841300 4(4) = - 1.C772156084656100 4(5) = - .498567C1940C35300 4(6) = - • 148437500 4(7) = - .029060570987654300 4(8) = - ,0 03 720 2 3809 52 381 DC 4(9) = - .00029968584656084700 4(10) = - .00001377865961199290-0- 4(11) = - .CCC00C275573192239€6D0 GC TC 235 211 4(1) = - .27426554003159900 ■Jg.^ _1 .464484 1269PA1 imp = - 1.1715145502645500 = - . 5793581 9QQ3S2 7 300 .18832 26 61 55202 eDO .0 414303626543 21000 .0062 11 144 79 £94 1800 .0 006 25 2066 7 9894 18000 4(10) = - .OOC040417401528512600 A(ll) = - .COCO 01 51 56525-573 192200 4(12) - - .COC000025C521C836544170C GC TO 235 212 4(1) = - .269028£467736492D0 ^ - -4~. -53*39 186 7 ?4 38 6 700 = - 1.2602711640211600 = - .659234 ie.2096265DQ- .2 3045 80 0264550300 .05569 72461 C52 322QC .0094 39 48412 69 84130 .001 1 192 7496 692 12 20 .C0C09C939 1534 3915 34DC .COCA 4 8225 30 £6 4 1 9 1 5 3 00 .COCOOC150 3126503 1265CDO 57 A ( 1 3 ) = - .C0C000002C676756S878681O0 GC TO 235 E( 2) = - l.DO GC TO 201 E<2) = - 1 .500 E( 3) = - .SCO CC TO 202 E(2) = - I. 833333 3 3 3 3.33 333D0 E( 3) = - 1 .00 E(4) = - .1666666666666667D0 GC TO 203 E(2) = - 2.C82333233223333Q0 E(3) = - 1 .4583333333233330C E(4) = - .4 16666666d6£6i66ZQC-- E<5) = - .04166666666666667D0 GC TC 204 8(2) - - 2.282333232223333DC E<3) = - l.e75D0 6(4) = - .7C8333333323333300 ECS) = - • 12 500 6(6) = - .00833333333232333300 GC TO 205 E(2) = - 2.45D0 £(3) = - 2.25S555555E555560C E(4) = - 1.G20832233323333D0 E(5) = - .24305S5555555556DC E(6| = - .C291666666666666700 E(7) - - .oci3eaEeeee£888889oo C-C TC 2C6 E(2) = - 2.5928571428571400 E(3) = - 2.6055S5555SS5556DC E(4) = - 1 .3430S555555S556QC - 6(5) = - .388688888888888900 E(6) = - .oe3eee£e8ee£ee889DO £(7) = - .OC5S5555555S5555560C E(e) = - .CCO 1984126S64126S800 GO TO 2C7 E(2) = - 2.71785714-26571400. „ _ EC3) = - 2.9296626984127C00 E(4) = - 1.66875CC E(S) = - .556770833332333300 E(6) = - .1 12500 E(7) = - .0 135416666666666700 E(6) - - .00C8926571 4.2B5X14 3 C E(<5) = - .000024601587201567300 C-C TC 2C8 E(2) = - 2.828968253 = -«2Afl 58 21 7592^9 2600 E<7) = - .043478009259259300 E<8) = - .00500 16E34a91EO4A0il E(9) = - .000363756613756614C0 E(10) = - .00C 01 £156525573 192200 E<11> = - .00000027557319223965900 24J) 231 E(2) = - 3.CI 967734487735D0 = - 3.7808134920634900 = - 2.6369367283950600 = - 1.15 229CC13 J 2i-l£XDQ_ .3 34 1834 766 3 1393 00 . Q 66.026 3f RP.t fflPflfl90Q . CC895 4 1 9973544 S74D0 .0CC818452JSaC9S23aiD0 E(10) = - .000048225306641975300 E(ll) = - .00000 165343915343915011^ E(12) = - .C3C0000250521083854417D0 CO TO £11 235 K = NC+1 ID CUB = K _ MTYP = (4 - MF)/2 ENG2 = .5/FL0AT(NC ♦ 11 ENC3 = .5/FL0AT**2 ED*N =(PERT£T < NG . MTYP «3 J *PEPShl±A2 IF (EDfcN.EQ.O) GC TO 780 ENC = FPR*FNQ.1/riF I nAT 240 IHEVAL = MF GO TO ( 250 , 68C ) .IRET C* THIS SECTION COMPUTES ThE EBEOICIEQ VALUES -EUL EFFECT IV ELX C* MULTIPLYING THE SAVEO INFORMATION BY THE PASCAL TRIANGLE C* MATRIX. Q* ************************************************************ ********* 250 T = T ♦ H CO 260 J = 2.K DO 260 J I = J.K J2 s K - Jl + J - 1 -— CO- 26.0 L =- 1 .N 260 Y(J2,I) = Y(J2.I) ♦ Y(J2+1.I> C**********4****t*********4*******4**** ********************** ********** C* LP 10 3 CORRECTOR ITERATICNS ARE TAKEN. CONVERGENCE IS TESTED C* EY REQUIRING CHANGES TO EE-LESS THAN END WHICH IVDEPENOENT ON C* THE ERROR TEST CONSTANT. C* THE SU M OF— THE- CO RR E CT I CNS I S A CCU MULA T ED IN THE A RR A Y C* ERRCR(l). IT IS EQUAL TC THE K-TH DERIVATIVE CF V MULTIPLIED C* BY H**K/(FACTGfil AL(K-1 )»A >) NT = NT - 420 CONTINUE IF CNT.LE.0) GC TC 490 430 CONTINUE C + + ** + *++ + + ***m + ++ + * + + + m+*i f ++4* + +4++** + ++4 + * + t + + + + + 4 + + + + + + ++ + m + **++m + C* THE CORRECTOR ITERATION FAILED TC CONVERGE IN 3 TRIES- VARIOUS C* PCSSIEILITIES ARE OECKE0 FOR. IF H IS ALREADY HMIN AND C* THIS IS EITHER ACAMS METHOD OR THE STIFF METHOD IN WHICH THE C* MATRIX PW HAS ALREADY BEEN RE-EVALUATED, A NO CONVERGENCE EXIT C* IS TAKEN. CThERWlSE THE VA TRI X Pin IS RE-EVALUATEn AND/OR THF C* STEP IS RECUCEC TO TRY ANO GET CONVERGENCE. c**»*» ♦♦»♦»♦♦♦♦»♦ »♦♦♦♦» ♦»»♦♦♦ ♦♦♦♦♦♦ ♦♦♦♦♦*♦♦♦♦»♦♦♦ ♦»»♦♦♦♦ »♦♦♦♦♦»♦♦»♦♦♦ »J 440 T = T - H IF ( ( h,LE. (h« IN* 1 .00CC1 ) ) .ANC . I { IWEVAL - M T YP ) .LT .- 1 ) ) GO TO 460 IF ( (MF.EQ.OI .OR. ( I WEVAL.NE .0)) RACUM = RACUM*0.25D0 IteEVAL = f*F IRET1 = 2 GC TC 750 46C KFLAG = -3 470 CO 480 I = 1,N CO 480 J - l.K 480 Y(J.I) = SAVEIJ.IX H = FOLD NQ = NGCLD JSTAPT = NG RETURN C* THE CORRECTOR CONVERfiFC ANC (TNTRni TS PASSED TO STATEMENT 520 C* IF THE ERROR TEST IS O.K., AND TC 540 OTHERWISE. C* IF THE STEP IS O.K. IT IS ACCEP T ED. I F-IODUfl H AS BEEN REDUCED _. C* TO CNE, A TEST IS MADE TC SEE IF THE STEP CAN BE INCREASED C* AT THE CURRENT ORDER OR BY COXKG-TQ- ONE HIGHER CR ONE LOWER-- C* SUCH A CHANGE IS ONLY MADE IF THE STEP CAN BE INCREASED BY AT -C-W-LEAST -X^l^ IF NU CH A N GE IS P O SS IBLE IOOUB I S S E T TO 1 T O C* PREVENT FUTHER TESTING FOR 10 STEPS. C* IF A CHANGE IS POSSIBLE. IT IS WACS AND IOQUB-J-S SET TO C* NQ ♦ 1 TO PREVENT FURTHER TESING FOR THAT NUMEER OF STEPS. C* IF THE ERROR WAS TOC LARGE* THE OPTIMUM STEP SIZE FOR THTS OR C* LCWER ORDER IS COMFLTED. ANC THE STEP RETRIED. IF IT SHOULD XA-EAIL TWT-CE- MQRF IT IS AN INDICATION THAT THE DERIVATIVES THAT C* HAVE ACCUMULATED IN THE Y ARRAY HAVE ERRORS OF THE WRONG ORDER C* SC THE FIRST CEfllVATIVES ARE RECOMPUTED AND THE CROFR IS ^SEI . 1 C* TC 1. c*» ♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦»j.»it»»»»Aj|.jt ♦»♦♦ ♦♦♦♦♦♦ ♦» »♦ »» »♦♦♦ »»^t »♦♦»♦♦»♦♦ »»jp 450 C at 0.0 ERROR(I) = CMEGA(I) ♦ CMEGA2(I) 500 D = D + (ERROR CI ) / YMAXX I^.U<*2 1 IHEVAL = 61 IF (C.GT.E) GO TO 54C IF (K.LT.3) dO TC 520 C** ******* 4 4************************************* *********************** C* COMPLETE THE CORRECTION CF THE HIGHER OROER DERIVATIVES AFTER A * C* SLCCESFUL STEP. * C**** * 4***4 4************************************************************ 00 510 J = 3.K CO 510 I - 1«* 510 Y(Jtl) = Y(J,I) 4 A < J)*CKEGA< I > ♦ B ( J ) *GMEGA2 ( I ) 520 KFLAG = +1 FNEW = H CXXXXXXXXXXPLEASE REMO VEX X X XX X X XX XX XXXXXXXXXX XXXXXXXXXXXXXX XXXXXXXXXXXXX MAXNG = MAXC ( MAXNC. K ) CXXXXX>XXXXOOES NCT BELONG HERE XXJLX XX X X X X X X XX XXXXJLXXX-XXXXXXXXXXXXXXXXXX X IF (ICCUB.LE.l) GO TO 550 IDCUE = ICCLE - 1 IF (ID0U8.GT.1) GO TO 7CC CO 530 I = 1 .N 530 £AVE(16.I) = EPR0fi(l> GC TO 700 Q********* ******************************************************* ******* C* RECOCL THE FAILURE FLAG COUNT TO CHECK FOR MULTIPLE FAILURES. * C* RESTORE T TO ITS ORIGINAL VALUE AND TRY AGAIN UNLESS THERE HAVE * C* THREE FAILURES. IN THAT CASE THE CERIVATIVES ARE ASSUMED TO HAVE * C* ACCUMULATED ERRORS SO A RESTART FRGM THE CURRENT VALUES OF Y I S * C*TRIED. _....._.. ._.____ * C + ************* ************************************************ ********* 540 KFLAG = KFLAG - 2 IF (h.LE. (HMN*1 .000C1) ) GO TC 740 T = TOLD IF (KFLAG. LE. -5) GO TC 72C C***44 4 4444 4*4***4 44* 44 4* 44 4* ********* *** ********* ********************** C* PR1. PR2. AND PR3 KkILL CONTAIN THE AMOLNTS BY WHICH THE STEP SIZE * C* COULD EE DIVIDED AT ORDER ONE LCteER. AT THIS ORDER. AND AT ORDER * C* ONE HIGHER RESPECTIVELY. * C*4*4 4* ************************************************ ***************** 550 FR2 = (C/E)**ENQ2*1 .2 FR3 = 1.E42C IF ( (NQ.GE.WAXDER). OR. (KFLAG. LE.-l) ) GO TO 57C C = .0 CO 560 I = 1 .N 560 C = C + ((ERROR(I) - SA VE ( I 6 • I ) J • YMAX( HI **Z PR3 = (D/EUP )**ENCJ*1 .4 57C PHI = 1.E4-2C - - IF (NG.LE.l ) GO TO 59C C = CO CO 5E0 I = l.N 580 C = D + ( f ( K. I )/YMAX{ I ) ) 44Z PR1 = (C/£D*N)**ENQ1*1 .3 590 CONTINUE __ - . _.. IF (PR2.LE.PR3) CO TC 650 IF (PR3.LT.FR1) GC TC 660 600 R = 1 .C/AMAX1 (PR 1 , 1 .E-4) NEteQ = NQ - 1 610 IDCUB = 10 IF ( ( KFLAG. EQ.l ) «AND^XR^LT-^-Ll . 1) ) ) GO J-O-700- IF (NE*Q.LE.NC) GC TC 630 C** ******** 4*** ************************************************** ******* C* COMPUTE ONE ACCITIONAL SCALED DERIVATIVE IF ORDER IS INCREASED. * 62 C****** ****************************************************************, -DC &2.Q La 1*M 620 YCNE*Q*1.I> = EfiRCP(I)*A(K)/DFLOAT(K) 630 K = NEWQ ♦ 1 IF (KFLAG. EC. 1) GC TO 670 KACUM = RACUM*R IRET1 = 3 _£D_IO 750 I 640 IF (NE*Q.EQ.NQ) GO TG 250 KG = NEWQ : GC TC 170 650 IF (PR2.GT.PR1) GO TO 600 NEteQ = NQ fi = 1.0/AMAXl(ER2»l^£2=AJ GC TC 610 660 R = 1 .0/AMAX1 (PR3.1 .E-4) NEteQ = NQ ♦ 1 GO TC 6 10 670 IRET = 2 fi = DMINKR .HMAX/r.AflF (H> ) F = F*R hNE* = H r IF (NO.EQ.NEtoQ) GO TO 68C NG = NEHQ GC TO 170 680 Rl = 1.0 CC 690 J = 2.K £1 = R1*R CO 650 I = l.N 690 Y(J,I) = Y(J.I)*fil IDOUB = K 700 CC 710 I = l^fel 710 YMAX(I) = DMAX1 (YMAX( I ) ,DAES( Y< 1 . I ) ) ) JSTART = NQ RETURN 720 IF (NC.EQ.l) GC TC 780 CALL DIFFUN ( T . Y • SA VE (N2 • 1 ) ) fi ..a HyFOLQ- CC 73C I = l.N Yll.IJ = SAVEC1.I* Ml = Nl 4 I SAVE (2, I) = hCLD*SAV£ = SAVE(2.I)*R Mi _^=_4 KFLAG = I GO TO 170 7AC KFLAG - -1 F.NE* = H JSTART = NC TURN C***** ***************************************************************** i C* THIS SECTION SCALES ALL VARIABLES CONNECTED WITH H AND R E TURNS J C* TC THE ENTERING SECTION. ' C****************************** ************** ************* *** ** ******** ) 750 RACUM = DMAX1 (CAfcSCHM IN/HOLD) .RACUM) -RACUM = DMlKUfiACLM,CABS(HMAX/HCLO) > Rl = 1.0 CC 760 J = 2.K Rl = R1*RACUM 780 63 DO 760 I = l«K 760 Y(J.I) = SAVEC J.I l*&± h - HCLD*RACUM DC 77C I = l«h 770 Y( 1 .1 > = SAVEC 1 »I ) ICCU6 = K CZ TC ( 130 KFLAG = -4 GC TC 470 £NC 25C . 640 ),IPET1 64 REFERENCES [ 1] BRAYTON, R. K. , GUSTAVSON, F. G., and HACHTEL, G. D. A new efficient algorithm for solving differential-algebraic systems using implicit backward differentiation formulas. Proa. IEEE 60 (1972), 98-108. [ 2] BROWN, R. L. Multi -derivative numerical methods for the solution of stiff ordinary differential equations. UIUCDCS-R-74-672, U. of Illinois, Urbana, IL, August 1974. (Also Ph.D. Th., Dep. of Computer Science, U. of Illinois, 1974.) [ 3] BYRNE, G. D., and HINDMARSH, A. C. A polyalgori thm for the numerical solution of ordinary differential equations. ACM Trans. Math. Software I, 1 (March 1975), 71-96. [ 4] ENRIGHT, W. H. Studies in the numerical solution of stiff differential equations. Tech. Rep. 46, Dep. of Computer Science, U. of Toronto, Toronto, Ont. , Canada, Oct. 1972. (Also Ph.D. Th., Dep. of Computer Science, U. of Toronto, 1972.) [ 5] ENRIGHT, W. H. Second derivative multistep methods for stiff ordinary differential equations. SIAM J. burner. Anal. II, 2 (April 1974), 321-331. [ 6] ENRIGHT, W. H. Optimal second derivative methods for stiff systems. In Stiff Differential Systems, R. A. Willoughby, Ed., Plenum Press, New York, 1974, 95-111. [ 7] ENRIGHT, W. H., HULL, T. E., and LINDBERG, B. Comparing numerical methods for stiff systems of o.d.e.'s. BIT 15, 1 (1975), 10-48. [ 8] FORSYTHE, G. E., and MOLER, C. B. Computer Solution of Linear Algebraic Systems. Prentice-Hal 1 , Englewood Cliffs, N.J., 1967. [ 9] GEAR, C. W. Algorithm 407: DIFSUB for solution of ordinary differential equations, Comm. ACM 14, 3 (March 1971), 176-179. [10] GEAR, C. W. Numerical Initial Value Problems in Ordinary Differential Equations. Prentice-Hall, Englewood Cliffs, N.J., 1971. [11] HENRICI, P. Disarete Variable Methods in Ordinary Differential Equations. Wiley, New York, 1962. [12] HINDMARSH, A. C. GEAR: Ordinary differential equation solver. UCID-3001, Rev. 3, Lawrence Livermore Laboratory, U. of California, Livermore, 1974. [13] KNUTH, D. E. The Art of Computer Programming, Vol. 1: Fundamental Algorithms. Addison-Wesley, Reading, Mass., 1968. 65 [14] KROGH, F. T. On testing a subroutine for the numerical integration of ordinary differential aquations. J. ACM 20, 4 (October 1973), 545-562. [15] LAMBERT, J. D. Linear multistep methods with mildly varying coefficients. Math. Comp. 24, (1970), p. 81-94. [16] LAMBERT, J. D., and SIGURDSSON, S. T. Multistep methods with variable matrix coefficients. SIAM J. Numer. Anal. 9, 4 (December 1972), 715-733. [17] SHAMPINE, L. F., and GORDON, M. K. Computer Solution of Ordinary Differential Equations: Initial Value Problems. Freeman, San Francisco, Calif., 1975. [18] SHAMPINE, L. F., and GORDON, M. K. Typical problems for stiff differential equations. ACM SIGNUM Newsletter 10, 2 & 3 (November, 1975), 41. [19] SKEEL, R. D. Equivalent forms of multistep formulas. In preparation, Dep. of Computer Science Tech. Rep., U. of Illinois, Urbana, IL. ormAEC-427 U.S. ATOMIC ENERGY COMMISSION plSfSoi UNIVERSITY-TYPE CONTRACTOR'S RECOMMENDATION FOR DISPOSITION OF SCIENTIFIC AND TECHNICAL DOCUMENT I See Instructions on Reverte Side ) I. AEC REPORT NO. C00-2383-0030 2. TITLE BLENDED LINEAR MULTISTEP METHODS 3. TYPE OF DOCUMENT (Check one): K] a. Scientific and technical report I | b. Conference paper not to be published in a journal: Title of conference Date of conference Exact location of conference Sponsoring organization □ c. Other (Specify) 4. RECOMMENDED ANNOUNCEMENT AND DISTRIBUTION (Check one): IY1 a. AEC's normal announcement and distribution procedures may be followed. ~2 b. Make available only within AEC and to AEC contractors and other U.S. Government agencies and their contractors. 3 c. Make no announcement or distribution. 5. REASON FOR RECOMMENDED RESTRICTIONS: 3. SUBMITTED BY: NAME AND POSITION (Please print or type) C. W. Gear Professor and Principal Investigator Organization University of Illinois at Urbana-Champaign Department of Computer Science Urbana, IL 61801 Si90atUre ^C& ^ '<2~- Date April 1976 FOR AEC USE ONLY r. AEC CONTRACT ADMINISTRATOR'S COMMENTS, IF ANY, ON ABOVE ANNOUNCEMENT AND DISTRIBUTION RECOMMENDATION: I. PATENT CLEARANCE: a. AEC patent clearance has been granted by responsible AEC patent group. □ b. Report has been sent to responsible AEC patent group for clearance. 1 I c. Patent clearance not required. IBLIOGRAPHIC DATA HEET 1. Report No. UIUCDCS-R-76-800 Title and Subtitle BLENDED LINEAR MULTISTEP METHODS 3. Recipient's Accession No. 5. Report Date May 1976 Author(s) Robert D. Skeel and Antonv K. Kong 8. Performing Organization Rept. No ' UIUCDCS-R-76-800 Performing Organization Name and Address Department of Computer Science University of Illinois at Urbana-Champaign Urbana, IL 61801 10. Project/Task/Work Unit No. 11. Contract /Grant No. US ERDA AT(ll-l) 2383 {. Sponsoring Organization Name and Address US ERDA Chicago Operations Office 9800 South Cass Avenue Argonne, IL 60439 13. Type of Report & Period Covered 14. Supplementary Notes . Abstracts The accuracy of linea systems is limited. Greater tives in the formula, but thi possible to duplicate the abs multistep formula by taking a blended formulas are similar with variable matrix coeffici detials and numerical results the backward differentiation r multistep formulas suitable for stiff differential accuracy can be attained by including higher deriva- s is not practical for all problems. However it is olute stability region for any given m-derivative linear combination of m-multistep formulas. These to Lambert and Sigurdsson's linear multistep formulas ents, but the approach is different. Implementation are presented for a blend of the Adams-Moul ton and formulas. Key Words and Document Analysis. 17a. Descriptors ordinary differential equations stiff equations multistep methods Adams method absolute stability linear multistep methods b. Identifiers/Open-Ended Terms e. COSATI Field/Group Availability Statement Unlimited RM NTIS-35 ( 10-70) 19. Security Class (This Report) UNCLASSIFIED 20. Security Class (This Page UNCLASSIFIED 21. No. of Pages 70 22. Price USCOMM-DC 40329-P7! JUNE1 REC'O