UNIVERSITY OF 
 
 ILLINOIS LIBRARY 
 
 AT URBANA-CHAMPAIGN 
 
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. 
 To renew call Telephone Center, 333-8400 
 
 UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN 
 
 JAW 2 5 "1935 
 
 JAN 2 ]9<J7 
 
 DEC 1 1 1996 
 
 L161— O-1096 
 

 
 
 
Oa . /o/</ 
 if. 1 
 
 maJ-ij 
 
 UIUCDCS-R-80-1014 
 
 UILU-ENG 80 1711 
 COO-2383-0067 
 
 March 1980 
 
 UNIFIED MODIFIED DIVIDED DIFFERENCE 
 IMPLEMENTATION OF ADAMS AND BDF FORMULAS 
 
 by 
 
 C. W. Gear 
 
Digitized by the Internet Archive 
 in 2013 
 
 http://archive.org/details/unifiedmodifiedd1014gear 
 
UIUCDCS-R-80-1014 
 
 UNIFIED MODIFIED DIVIDED DIFFERENCE 
 IMPLEMENTATION OF ADAMS AND BDF FORMULAS 
 
 by 
 
 C. W. Gear 
 
 March 1980 
 
 DEPARTMENT OF COMPUTER SCIENCE 
 
 UNIVERSITY OF ILLINOIS AT URB ANA- CHAMPAIGN 
 
 URBANA, ILLINOIS 61801 
 
 Supported in part by the U.S. Department of Energy, Grant ENERGY/EY-76-S-02-2383 
 
ABSTRACT 
 
 The representation of both full variable-step and fixed-leading- 
 coefficient implementations of BDF formulas using the modified divided 
 differences normally associated with Adams' formulas is developed. 
 Simple recursion relationships for the coefficients are derived. The 
 technique permits a simple, compact implementation of variable-step and 
 fixed-leading-coefficient forms of blended formulas. 
 
- 2 - 
 
 J.. INTRODUCTION 
 
 A careful implementation of modified divided differences (MDDs) 
 leads to an efficient implementation of variable-order, variable-step 
 multistep codes (see Krogh [6] and Shampine and Gordon [8]). These 
 codes can be more efficient than use of the Nordsieck vector (Gear [2] 
 
 and Hindmarsh [4] ) because the Nordsieck vector approach uses order of 
 
 2 
 
 k /2 additions and stores in the predictor phase plus an additional k 
 
 multiplications and additions and stores in the corrector phase. (These 
 numbers are per component of a system. We will ignore other operations 
 that are independent of the number of components.) In contrast, MDDs use 
 2k additions and k multiplications in the predictor phase followed by k 
 additions and stores in the corrector phase. To this must be added k - 
 s multiplications and stores if the step has been constant for only s 
 steps. By keeping the step constant whenever possible, MDDs have 
 considerably less work per component for any but the smallest values of 
 k, and even if the step changes frequently there is a saving for modest 
 to large k. Although there is some additional work in coefficient 
 calculation for MDD-based codes making them perhaps less attractive for 
 a very small number of equations, those problems do not normally use 
 much computer time anyway. (Also, if the step is constant for a while, 
 the coefficients become constants, saving the overhead.) 
 
 One drawback of MDDs is that a substantially different scheme must 
 be used for Adams' formulas (for nonstiff equations) to BDF formulas 
 
- 3 - 
 
 (for stiff equations). When Adams' formulas are used, the MDDs of y' = 
 dy/dt are stored, while when BDF formulas are used, the MDDs of y are 
 stored. This means that when a switch from a nonstiff to a stiff method 
 is wanted (or vice-versa) a significant disruption is necessary in the 
 processing. Nordsieck-based codes, on the otherhand, simply change some 
 coefficients and continue. This makes a "continuous transition" from 
 nonstiff to stiff, as done in the blended methods of Skeel and Kong [9] , 
 a reasonable addition. 
 
 This paper gives a technique for implementing a variable-step form 
 of the BDF formulas while storing the MDDs of y' (instead of y). The 
 terra "variable-step form" is used in a technical sense here to mean that 
 the technique is exactly equivalent to using a BDF formula based on 
 unequal intervals. Formulas based on unequal intervals have been shown 
 to have generally better stability properties than those based on 
 constant steps with interpolation for step changing (see Gear and Tu 
 [3]). The latter method is the simplest to use with the Nordsieck 
 scheme although the former have been used (Byrne and Hindmarsh [1]). 
 The former method is more suited to MDD-based codes . 
 
 Although this technique is more expensive than storing the MDDs of 
 y (by k multiplications and additions), it allows a simple transition 
 from Adams to BDF formulas and thus permits a simple implementation of 
 the blended methods. The additional cost is not too significant in BDF 
 methods because the corrector normally requires the solution of implicit 
 
- 4 - 
 
 equations, and this is the most time-consuming part of the process. 
 
 The paper also describes the very simple change necessary to get a 
 fixed-leading-coefficient form of the method, and the application of 
 this form to a "Blended" version. 
 
 The second section of this paper describes the basis of the 
 technique. This is illustrated by a simple fixed-step second order BDF 
 method. The third section discusses the calculation of the BDF 
 coefficients in a variable-step method and develops difference equations 
 which are simple to evaluate. The final section discusses fixed- 
 leading-coefficients and blended methods. 
 
- 5 - 
 
 2. POLYNOMIAL INTERPOLATION AND UPDATING 
 
 This section first describes the basic MDD approach to an Adams' 
 formula-based code in general terms. (Details can be found in [6] and 
 [8].) From this viewpoint it is then clear how a BDF formula-based code 
 can be implemented storing the same quantities as the Adams' scheme. 
 Finally this Is illustrated with second-order examples. 
 
 Suppose that the solution has already been calculated at a set of k 
 
 points t , t ...... t ,,, and we have kept the values 
 
 r n n-1 n-k+1 r 
 
 y y', y' ,,..., y' , , . . In the MDD-based method we actually store the 
 J n, J n y n-l ^n-k+1 J 
 
 1* 
 modified divided differences V J y' , j *» 0, . . . , k-1, which, for constant 
 
 steps, are equal to the backward differences of order j. The first 
 
 process in calculating the value at t , , ■ t + h is prediction, and can 
 r n+1 n 
 
 be written as 
 
 k-1 
 J-0 
 
 (1) 
 
 'n+1 * 7 n y n 
 
 k_1 4* 
 P^i = y + h V Y 4 V J y' 
 r n+l 'n A L n In n y n 
 
 j-0 J 
 
 (2) 
 
 where the y^ are the Adams-Bashforth coefficients (which are functions 
 jn 
 
 of the stepsize ratios). It is important to note that prediction 
 
 consists of nothing more than calculating the value and derivative at 
 
 t . of the k-th degree polynomial which matches y at t and y' at 
 
 t ., < I < k-1. The second process is to calculate a correction 
 n-i r 
 
- 6 - 
 
 term. This depends on which type of method, namely PE, PEC, PECE,..., 
 
 or solving the corrector by Iteration, is used. Whatever method is 
 
 used, we are computing the value y^, which is such that ?' +1 - f(y_.i) 
 
 for the differential equation y' » f(y), where ?',, is defined as the 
 
 derivative of the polynomial passing through 
 
 y .i» y y' t y' 1 1 • • • » y' • The corrector order is R+l and normally 
 y n+l •'n, J n ^n-1 J _. 
 
 n-K 
 
 R » k. If the corrector is iterated to convergence, $',-. ■ f(y^_i) anc * 
 
 y',, = y'., • For the other methods we have 
 y n+l y n+l 
 
 PE! K+i - "*i s ° 1*1 - *Vl 
 
 1*1 • ^.H-l' 
 
 PEC! Kn ' 1*1 - £( >Vi> 
 
 PECE: r a+l ■ Up^) 
 
 Finally we have to update the saved Information so that we have the 
 
 representation of y .,, y' , , , y' , . . . , y' , . „ . for the next step if the 
 
 n+1 n+1 J n y n-k+2-j 
 
 order is the same (j = 0), increased by one (j = 1) or decreased by - j . 
 
 The values of y . , and y' +1 are defined above. The remaining MDDs can 
 
 be updated by observing that y' . - p' . is the difference between y' . 
 
 and the derivative at t . of the k-th degree polynomial whose 
 
 derivatives are y' at t . , < I < k. Hence this is exactly the k-th 
 n-i n-i J 
 
 modified divided difference of y' , -1 < 1 < k, V , , y' . From this the 
 
 •'n-j' J ' n+1 J 
 
 iteration 
 
 1 i* 1+1 
 V-J » V J + V J ^J j - k-1,..., 
 n+1 n n+1 j » » 
 
- 7 - 
 
 can be used, where the MDD r M is the 1-th order MDD at t . , based on 
 
 n+1 n+1 
 
 h ■ t . - t . It now has to be multiplied by a step-ratio dependent 
 
 1* 
 scale to get V J , the MDD at t , , based on the next stepslze 
 n+1 n+1 r 
 
 t _ - t .. Details are given In Krogh or Shamplne and Gordon op_. cit . 
 
 After the update in Eq. (3), the MDD table contains one additional 
 
 entry, V and represents a polynomial through y' , -1 < j < k. This 
 
 is appropriate if the order is to be increased by one. If not, we must 
 
 discard y' , . , y' , ., ■ . • as necessary — that is to say we must update 
 
 the MDD table so it represents a polynomial of minimum degree passing 
 
 through the saved values. For Adams' method the saved values are the 
 
 most recent derivatives so all that Is necessary is to discard the 
 
 highest order MDDs of y' . . 
 
 n+1 
 
 Now let us examine the BDF formulas using the same representation. 
 
 Suppose that we have already computed y,y .,..., y ,, and are 
 
 representing this information by saving y and the MDDs of the 
 
 derivatives at t , < j < k, of the k-th degree polynomial passing 
 
 through y , < 1 < k. The prediction process again consists of 
 
 calculating the value and derivative at t . of the polynomial. The 
 
 predicted values are given by the same formulas (1) and (2) as the 
 
 Adams' code. The correction process is also identical (usually 
 
 consisting of solving the corrector equation starting with the predicted 
 
 value and using a quasi-Newton iteration). However, note that we now 
 
 get the value of y ,, such that ?',,, the derivative at t (1 of the 
 y n+l ■'n+l n+1 
 
 polynomial passing through y . , -1 < 1 < k, approximately satisfies 
 
- 8 - 
 
 f , ■ f(y .,). The main difference comes in the updating process. Now 
 n+ln+1 
 
 we want to calculate the MDDs V y' which are the MDDs of the 
 
 n+1 n+1 
 
 derivatives of the polynomial passing through y , -1 < i < k-l+j where 
 j is +1 if the order is increased by one, if the order is unchanged, 
 and negative if the order is to be reduced. 
 
 How can we do this? Let us first consider the case of increased 
 order. At the start of the step the saved values represent the unique 
 
 k-th degree polynomial passing through y , < i < k. For such a 
 
 k* 
 polynomial, V y = 0. At the end of the step the saved values must 
 
 represent the unique (k+l)-st degree polynomial passing through 
 
 y , -1 < i < k. Hence, the polynomial will have been modified by the 
 
 addition of a (k+l)-st degree polynomial which is zero at 
 
 t . , < i < k and non-zero at t . , (for example, 
 n-i , n+1 r 
 
 k 
 
 P ,(t)» IT (t-t ) or any scalar multiple of it). Hence, the new 
 n,k 1=Q n-l 
 
 MDDs at t are 
 n 
 
 (4) 
 
 V J 0' * V J + aV J P' , 
 n y n+l n n n,k 
 
 1* 
 
 for some scalar a, where by V J 9',, we mean the MDDs at t of the 
 
 n y n+l n 
 
 derivative of the polynomial through y ., > y > y 1 »»««» y t « The 
 
 coefficient a must be chosen so that 
 
 P n+1 + aP n,k (t n+l ) ~ ^n+1 (5) 
 
 and 
 
- 9 - 
 
 »; + i + ^.k'^i' a £( w 
 
 Finally, we confute V^_ ?' from V^ ?' . using Eq. (3). 
 
 This gives the MDDs when the order is raised by one. What happens 
 if it is to stay the same or is lowered? We will analyze this by 
 considering how to solve the following problem: Given the MDDs of the 
 derivative of the k-th degree polynomial through y , < i < k, how 
 can we find the MDDs of the derivative of the (k-l)-th degree polynomial 
 through y , < i < k? If we can do this, we can apply the process 
 repeatedly to lower the order as many times as needed. 
 
 The modification needed consists of the addition of a k-th degree 
 
 polynomial which is zero att ., < i < k to the original polynomial 
 
 n— l 
 
 to reduce its degree to k-1. Thus we add a multiple of 
 
 k-1 
 
 P , . - n (t - t .). 
 n,k-l i=Q n-i 
 
 The modification to the MDDs are 
 
 1* -t* i* 
 V J y' <- V J y' + aV J P' , , 
 nn nn n n,k-l 
 
 (7) 
 
 k-1* 
 such that the new value of V y' is zero. Hence, 
 
 n J n 
 
 yJ* y - + v j * y' - V j * P' -—2- 2— , j - 0,..., k-1 
 
 n n n y n n n,k-l yk-1* _, J 
 
 n,k-l (8) 
 
 Thus, the integration process for the BDF formula is as given below. In 
 
- 10 - 
 
 this process we will write 8 for V P' , <J> , for / 6 , , and 
 
 n,k n n,k n,k . . n,k 
 
 -1 1 ^ 
 
 <> . for P . (t , .) (<+> J , is the i-th order MDD of P' , at t , , based on 
 n, k n,k n+l n,R n,K n+l 
 
 h - t . . - t ). 
 n+l n 
 
 Single Step 
 
 (a) Predict using Eqs. (1) and (2). 
 
 (b) Correct by "solving" 
 
 p' + ot<}> « f (p + a<}> ) 
 
 n+l n,k r n+l n,k 
 
 for a . 
 (c) Update MDDs by 
 
 (9) 
 
 i* i* i 
 v f\ , - V y' + a6 x , < i < k 
 n n+l n n n,k 
 
 k* k 
 n n+l n,k 
 
 (d) Advance MDDs to t . , by 
 
 n+l 
 
 V k v' =» V k * 9' 
 n+l y n+l n y n+l 
 
 i i* i+1 
 
 v \i y j.i ■ v 9'j.i + v It y\i » i - k-i,..., o 
 
 n+l y n+l n y n+l n+l y n+l 
 
 i* 
 
 (e) Modify the divided differences to get V y'+i • 
 
 (f) Lower order if necessary by 
 
 - „k* , /n k 
 n n+l 7 n+l n+l,k 
 
 i* i* i 
 
 n+l y n+l n+l y n+l n n+l,k 
 
 This order lowering can be repeated for k-1, k-2,..., as 
 
 required. 
 
 (10) 
 
 (ID 
 
 (12) 
 (13) 
 
- 11 - 
 
 In fact, (c) and the last application of (f) from the previous step can 
 be combined to reduce arithmetic. This is described in the next section 
 on coefficient calculation. 
 
 Example . 2nd order BDF method with constant stepsize. 
 
 and 
 
 This starts with the values y , , y which are represented by y 
 
 'n-l'n r J J n 
 
 y' = (y - y ,)/h. P , and its MDDs are given by 
 J n J n J n-l n,l G J 
 
 P - (t - t ) (t - t + h) 
 
 n, i n n 
 
 e° - h , 8* - a , 
 
 n,l n, l 
 
 4> , ■ 3h , <{> , = 2h 
 n,l n,l 
 
 The predictor is 
 
 P^i - y + hy' ( = 2y - y . ) 
 r n+l ; n J n J n y n-l 
 
 p 0+ i - "„ 
 
 The corrector solves 
 
 p; +[ + 3ho = «p + 2h 2 a) 
 
 or 
 
 4 1 2h mt , 
 
 y n+l '3 y n" IVl + T f ( W 
 
 where 
 
- 12 - 
 
 y n+l - P n+1 y n+l " 2y n + y n-l 
 . - 
 
 2ti 21i 
 
 This is precisely the 2nd order BDF method. At the increased order we 
 update the MDDs to get 
 
 followed by 
 
 n y n+l y n 2h 
 
 _ y n-t-l " y n-l 
 2h 
 
 and 
 
 v l* a y n+l ' 2y n + y n-l 
 n y n+l h 
 
 v l , y n+l ~ 2y n + y n-l 
 n+1 y n+l ' h 
 
 3y n+l " 4y n + y n-l 
 
 7° y 
 n+1 y n+l 2h 
 
 i* 
 At fixed stepsizes there is no further modification to get V . , and we 
 
 can lower the order by computing 
 
 ,0* , „0 n l , Q /Q 1 
 
 7 
 
 n+1 y n+l n+1 y n+l n+1 y n+l n+l,l' n+1,1 
 
 y n+l ' y n 
 
- 13 - 
 
 and 
 
 Note that the BDF corrector has one higher order than the predictor-. 
 This is different to the BDF formulas implemented in DIFSUB [2] where 
 the predictor and corrector have the same order. 
 
- 14 - 
 
 3. COEFFICIENT CALCULATION 
 
 Let P . (t) - IT (t - t ). Define Q , (t) - P' (t). It is easy 
 
 n, k in n n, k n, k. 
 
 to see that 
 
 «n.k (t) " (t " 'n-k' "n.k-l (t) + P n.k-l (t) (M) 
 
 Let the i-th order divided difference of Q , at t be w ' . We have 
 
 n,k m n,k 
 
 ^l - Q t (t ) 
 
 n,k ^n,k m ^^ 
 
 M»»J +1 - [J*'} •\r m }' t ]/(t - t . ,). 
 n,k n,k n,k m m-l-i .. , s 
 
 Ultimately we want to compute the modified divided differences of Q at 
 
 t . First we will examine the regular divided differences w ' . We 
 n ° n,k 
 
 would like to generate these simply so we look for a recursion on k 
 
 using Eq. (14). Noting that P , , (t ) - for n-k+1 < m < n we have 
 
 n , k— 1 m 
 
 ""Ik - < e . - Vk' Vk-i'V 
 
 for the same range of m. Now we claim that 
 
 w m »j . (t - 1 .) *"•£ . + w^*- 1 
 
 n,k m n-k n,k-l n,k-l 
 
 - (t . - t . )w m »* . + W*'*" 1 
 m-i n-k n,k-l n,k-l 
 
 (17) 
 
 if n-k+i < m < n and 1 < i. 
 Proof 
 
 For 1 ■ 1 we have by definition 
 
- 15 - 
 
 W* 
 
 ,1 (t m - t n-k ) Q n.k-l (t m ) " ( Vl " Vk) Sufe^fcl? 
 
 n,k t - t . 
 
 m m-1 
 
 U m " n-k' t - t , 
 
 m m-1 
 
 + Vw'Vi' 
 
 or 
 
 ,. , q n.k-l (t „> " ^.k-l'Vl* 
 
 U m-1 " n-k' t - t . 
 
 m m-1 
 
 + Vk-lV 
 In either case the result follows from Eq. (15). 
 
 Now consider 1 > 1 and use induction. Substituting Eq. (17) into 
 
 the definition of w * we get 
 
 n,k ° 
 
 4 W*'*" 1 - W*-?' 1 - 1 
 
 «*» i , n » k H*J£ 
 
 n,k t - t . 
 
 ra m-i 
 
 ct - 1 . ) *■»£-} - (t , - 1 .) w^?^- 1 
 
 m n-k n,k-l m-i n-k n.k-1 
 
 m m-i 
 
 Q.E.D, 
 from which the result follows directly. 
 
 The modified divided differences are given by 
 
- 16 - 
 
 n,k nri-1 m nrfl m+l-i n,k 
 For our purposes it will be convenient to introduce a scale of P , and 
 
 Dy K 
 
 its divided differences so that Q , (t ) » P' , (t ) - 1. This is 
 
 Ti,k m n,k n 
 
 achieved by dividing by IT (t - t .). Hence, we define 
 
 1-1 n Brl 
 
 m,i _ (t mfl "* ^ * " (t ari-l " t mfl-i ) „m,i 
 
 n,k " (t - t )...(t - t ) n,k 
 n n— I n n— k 
 
 The 6 coefficients defined earlier are then given by 9 ■ ' . Eq. 
 
 LI y K U. j K 
 
 (17) now yields the recursion 
 
 Q i _ (t n-I " fc n-k } 6 n.k-l + (t n+l " t n+l-i ) 6 nTk-l 
 
 n,k " t - t , 
 
 n n-k 
 
 (18) 
 
 1 < i<k 
 Because of the scaling, 6 , - 1. We will use the notation of Krogh 
 (op » cit «), namely, 
 
 V n+1) " Sh-i - Vi 
 
 6 (n+l) - 1 
 
 P ± (n+1) - B^^n+1) C 1 . 1 («H-l)/5 1 _ 1 (n) 
 
 from which we find that 
 
 "n.k " V n+1 > "S.'k 
 
 It is easy to show that W 11 ', - k+1. (It is the k-th divided difference 
 
 n,k 
 
 of the derivative of a monic (k+l)-st degree polynomial.) Hence, we know 
 
 Ok i 
 
 9 , and 9 , k - 0, 1 so we can calculate 9 using the 
 
 n, k n, K n, Ic 
 
- 17 - 
 
 recursion Eq. (18) rewritten as 
 
 < IC.(n+l) - C.(n+1)] 6* + ? (n+1) 6*~) 
 e 1 , J 1 n.1-1 i-1 n.1-1 
 
 a, J Cj(n+1) - C (n+1)' 
 
 (19) 
 
 for i - 1,..., j-1; j - 2,..., k . Note that if the steps izes have been 
 
 constant for s steps (that is, 
 
 'nfl " fc n ' fc n - C n-1 — - ^2-8 " t nfl-s ) » then 6 j (n+1) " X and 
 
 e ■ (j+D/(j+l-i) for j < s. Thus we can start with j - s - 1 and 
 n » J 
 
 calculate 6 explicitly for i - 1,..., j. If s > k, the previously 
 
 calculated values of 9 , , can be used. We also need 
 
 n-l,k 
 
 *n,k ' j Q 9 n,k " J-0.....k 
 
 (20) 
 
 and 
 
 ♦;!k - v* +i > e k (n+i) 
 
 (21) 
 
 in the corrector step 
 
 Finally, note that the action represented by Eq. (13) can be 
 
 delayed until the start of the next step when Eqs. (1), (2), and (10) 
 
 1* 
 must be executed. Suppose that V J y' are the MDDs left by Eq. (11) in 
 
 the previous step update to the new stepslze. Eqs. (1) and (2) can be 
 
 written as 
 
- 18 - 
 
 *n+l /; n n 'n n-1 n,k 
 
 J-o 
 
 J* -' _ a ,.° 
 
 (22) 
 
 *- rt n ; n n-1 n,k n,k 
 
 J-o 
 
 k_1 i* 1 
 
 p^i - y + h I Y 4 (V J y' - 6 . 9 J ) 
 r n+l 'n .^ jn n 'n n-1 n,k 
 
 k-1 k-1 
 
 » y + h [ y. V J y' - h6 , £ y 4 J 
 
 n j-o Jn n n n " 1 j-o Jn n,k 
 
 The order lowering represented by Eq. (13) can now be handled by 
 changing the first member of Eq. (10) to 
 
 V 1 * ?' . - V 1 * y' + (a - 6 , ) 6 1 
 n 7 n+l n •'n n-1 n,k 
 
 (23) 
 
 (24) 
 
 Alternatively, the 6 terms can be omitted from Eqs. (22) and (23) if 
 the summations are taken from to k. This corresponds to using a 
 predictor of one order higher — which is possible if the order did not 
 increase in the previous step. 
 
 Examination of the equations given above reveals that several 
 divisions can be saved by introducing another normalization of so that 
 9 ■ l. Our program computes using Eq. (18) and then makes this 
 
 11) Iv 
 
 final normalization before computing $ and <J> 
 
- 19 - 
 
 4.. FIXED LEADING COEFFICIENT FORMULAS AND BLENDED METHOD 
 
 The corrector equation which must be solved at each step for stiff 
 equations is 
 
 p' + a<J> « f (p + a<j> ) 
 
 (25) 
 
 where we have dropped the numerous subscripts. For stiff equations f 
 
 is large so a quasi-Newton iteration is often used to solve Eq. (25). 
 
 This means that the matrix [I - J<(> /$ ] must be decomposed, where 
 
 J « f . Unfortunately, in variable-step methods <j> /<J> is only constant 
 
 if the step and order have been fixed for at least k steps. This is in 
 
 contrast to interpolatory methods, usually implemented with Nordsieck 
 
 vectors, in which <j) /<}> is a function of the current stepsize and order 
 
 only. In the latter methods, it is not necessary to redecompose 
 
 [I - J<J> /<J> ] unless the stepsize or order changes. In [5], Jackson and 
 
 Sacks-Davis discuss fixed-leading-coefficient methods. These methods 
 
 have the advantage, along with interpolatory methods, that the ratio 
 
 <t> /<)> does not change unless the stepsize or order changes. There are 
 
 two equivalent ways of viewing these methods. In a variable-step 
 
 framework, we want to use the k-th order formula 
 
 k+1 
 
 n ° n i»i n* 1 a -1 
 
 (25) 
 
 where P. is the coefficient for the fixed-step, k-th order, k-step BDF 
 
 method. This determines the a uniquely unless all stepsizes are 
 
- 20 - 
 
 equal, in which case we get uniqueness (and continuity) by choosing 
 
 a , . , ■ 0. Another view of the method is to start with the unique k-th 
 n,k+l 
 
 degree predictor polynomial through y , 1 < i < k+1, interpolate 
 
 values of y" . to a set of k equally spaced points at 
 
 t . - (i-l)h, 1 < i < k, and use the fixed-step, k-th order BDF formula 
 
 k 
 
 y - B ft hy' + £ a. f . 
 y n ; n L . i 'n-i 
 
 (26) 
 It is simple to show that these are equivalent by considering the 
 expansion of J into linear sums of the y to convert from Eq. (26) 
 
 to Eq. (25). This is a unique mapping and since Eqs. (25) and (26) are 
 uniquely determined by the k-th order criterion; the results must be 
 identical. 
 
 This fixed-leading-coefficient form can be implemented using the 
 MDD scheme as follows: 
 
 a. Predict using Eqs. (1) and (2) with sums to k. 
 
 b. Solve the corrector equation 
 
 y' - p' + a$° - f(p + a*" 1 ) - f(y) 
 where 4> and <J> are the fixed-step coefficients, 
 c Update using Eqs. (24) and (11). 
 
 Blended methods [9] can be implemented by noting the difference 
 
- 21 - 
 
 between Adams and BDF methods arise only in the definition of 8 and <J>. 
 
 B B 
 Let us write 9 and <J> for the coefficients introduced above. The 
 
 equivalent Adams coefficients are 
 
 eA, i - 6 <u » < i < k 
 n,k ik 
 
 where 6., is the Kronecker delta, 
 ik 
 
 n,k 
 
 ♦ i - hy . 
 n,k n,k 
 
 k* 
 where y , is the Adams-Bashforth coefficient of 7 y' . (As before, we 
 n,k n ■'n 
 
 can have a fixed-leading coefficient form by choosing y , ■ y.» the 
 fixed-step-coefficient.) In blended methods we form the linear 
 combination 
 
 A - rJB =» 
 
 where A ■ represents the Adams method, B - represents the BDF 
 method, J is the (approximate) Jacobian^ and r is a real scalar. For 
 constant J we can save the MDD representation of [I - rJ] (P. - rJP R ) 
 where P and P n are the polynomials through the saved values for Adams 
 and BDF methods, respectively. To accomplish this, we simply solve 
 
 p' + <!> A ' a - r<> B ' ja - f(p + .J,*'" 1 a - r* 8 '" 1 ! a) 
 
 and then update the MDD's using Eq. (10) modified to 
 
- 22 - 
 
 V 1 * y' . v 1 * y' + a6 A 'J - rJa9 B »J 
 n 'n+1 n J n n,k n,k 
 
 As before, this can be modified to combine steps to get the replacement 
 
 for Eq. (24) by replacing a by (a - 6 ) In Eq. (27). With the 
 
 It 
 
 normalization of 6 , ■ 1 (for Adams and BDF), 6 , is precisely the a 
 
 n,k n-1 r J 
 
 from the previous step. 
 
- 23 - 
 
 REFERENCES 
 
 [1] Byrne, G.D. and A.C. Hindmarsh, A Poly algorithm for the Numerical 
 Solution of Ordinary Differential Equations, ACM Trans . 
 Mathematical Software 1 (1), March 1975, 71-96. 
 
 [2] Gear, C.W. , DIFSUB for Solution of Ordinary Differential Equations 
 [D2], Communications ACM 14 (3), 1971, 185-190. 
 
 [3] Gear, C.W. and K.W. Tu, The Effect of Variable Mesh Size of the 
 Stability of Multistep Methods, SIAM J. Numer . Anal . 11 (5), 
 October 1974, 1025-1043. 
 
 [4] Hindmarsh, A., GEAR: Ordinary Differential Equation Solver, Report 
 UCID-30001, Rev. 3, Lawrence Livermore Laboratory, CA, 1974. 
 
 [5] Jackson, K.R. and R. Sacks-Davis, An Alternative Implementation of 
 Variable Stepslze Multistep Formulas for Stiff ODEs, University of 
 Toronto, Dept. Comp. Science Report 121, July 1978. 
 
 [6] Krogh, F.T., Changing Stepsize In the Integration of Differential 
 Equations using MDDs, in Lecture Notes in Mathematics 362 , 
 Springer-Verlag, NY, 1974, 22-71. 
 
 [7] Sacks-Davis, R., Fixed Leading Coefficient Implementation of Second 
 Derivative Formulas for Stiff ODEs, to appear ACM TOMS. 
 
 [8] Shampine, L.F. and M. Gordon, Computer Solution of Ordinary 
 
 Differential Equations : the Initial Value Problem, W.H. Freeman: 
 San Francisco. 
 
 [9] Skeel, R.D. and A.K. Kong, Blended Linear Multistep Methods, ACM 
 Trans . Mathematical Software , _3 (4), December 1977, 326-345. 
 
FormAEC-427 U.S. ATOMIC ENERGY COMMISSION 
 
 appm^Jqi UNIVERSITY-TYPE CONTRACTOR'S RECOMMENDATION FOR 
 
 DISPOSITION OF SCIENTIF C AND TECHNICAL DOCUMENT 
 
 ( See Instructions on Reverse Side ) 
 
 1. AEC REPORT NO. 
 
 COO-2383-0067 
 
 2. TITLE 
 
 UNIFIED MODIFIED DIVIDED DIFFERENCE IMPLEMENTATION 
 OF ADAMS AND BDF FORMULAS 
 
 3. TYPE OF DOCUMENT (Check one): 
 
 a. Scientific and technical report 
 
 1 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): 
 
 [xj 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. 
 I | 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 U-C 
 Urbana, IL 61801 
 
 2 
 
 Signature 
 
 ±^=. — ,^ % — V — ^s 
 
 Date 
 
 March 1980 
 
 FOR AEC USE ONLY 
 
 7. AEC CONTRACT ADMINISTRATOR'S COMMENTS, IF ANY, ON ABOVE ANNOUNCEMENT AND DISTRIBUTION 
 RECOMMENDATION: 
 
 8. PATENT CLEARANCE: 
 
 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. 
 
 1 I c. Patent clearance not required. 
 
JIBLIOGRAPHIC DATA 
 >HEET 
 
 1. Report No. 
 
 UIUCDCS-R-80-1014 
 
 2. 
 
 3. Recipient's Accession No. 
 
 . Title and Subtitle 
 
 UNIFIED MODIFIED DIVIDED DIFFERENCE IMPLEMENTATION 
 
 5. Report Date 
 
 March 1980 
 
 OF ADAMS AND BDF FORMULAS 
 
 6. 
 
 . Author(s) 
 
 C. W. Gear 
 
 8. Performing Organization Rept. 
 
 No. R-80-1014 
 
 . Performing Organization Name and Address 
 
 Department of Computer Science 
 
 10. Project/Task/Work Unit No. 
 
 University of Illinois U-C 
 Urbana, IL 61801 
 
 11. Contract /Grant No. 
 
 EY-76-S-02-2383 
 
 2. Sponsoring Organization Name and Address 
 
 Department of Energy 
 9800 S. Cass Avenue 
 
 13. Type of Report & Period 
 Covered 
 
 technical 
 
 Argonne, IL 
 
 14. 
 
 5. Supplementary Notes 
 
 6. Abstracts 
 
 The representation of both full variable-step and fixed-leading-coefficient 
 implementations of BDF formulas using the modified divided differences 
 normally associated with Adams' formulas is developed. Simple recursion 
 relationships for the coefficients are derived. The technique permits a 
 simple, compact implementation of variable-step and fixed-leading- 
 coefficient forms of blended formulas. 
 
 7. Key Words and Document Analysis. 17a. Descriptors 
 
 integration 
 
 differential equations 
 multistep methods 
 blended formulas 
 
 7b. Identifiers/Open-Ended Terms 
 7c. COSATI Field/Group 
 
 8. Availability Statement 
 
 unlimited 
 
 19. Security Class (This 
 Report) 
 
 UNCLASSIFIED 
 
 21. No. of Pages 
 27 
 
 20. Security Class (This 
 Page 
 
 UNCLASSIF 
 
 22. Price 
 
 ORM NTIS-35 ( 10-70) 
 
 
 
 
 
 USCOMM-DC 40329-P7 1 
 
FE&il 
 
 n\