Digitized by the Internet Archive 
 in 2013 
 
 http://archive.org/details/equivalentformso940skee 
 
O///;UIUCDCS-R-78-940 
 WW 
 
 lYl&bL 
 
 UILU-ENG 78 1730 
 
 EQUIVALENT FORMS OF MULTTSTEP FORMULAS 
 
 by 
 
 Robert D. Skeel 
 
 August 1978 
 
UIUCDCS-R-78-940 
 
 EQUIVALENT FORMS OF MULTISTEP FORMULAS 
 by 
 Robert D. Skeel 
 
 August 1978 
 
 DEPARTMENT OF COMPUTER SCIENCE 
 
 222 DIGITAL COMPUTER LABORATORY 
 
 UNIVERSITY OF ILLINOIS AT URBANA- CHAMPAIGN 
 
 URBANA, ILLINOIS 61801 USA 
 
 Research sponsored by the Air Force Office of Scientific Research, Air Force 
 Systems Command, USAF, under Grant No. AFOSR-75-2854. The United States 
 Government is authorized to reproduce and distribute reprints for Governmental 
 purposes not withstanding any copyright notation hereon. 
 
1. Introduction. Multistep methods are the oldest and the most important 
 
 class of numerical methods for solving systems of ordinary differential equations, 
 
 Implementations of these methods have become increasingly sophisticated 
 
 in the last two decades. One paper having a considerable impact on the 
 
 development of these algorithms is that of Nordsieck [1962], which 
 
 introduces a class of formulas closely related to multistep formulas. 
 
 It is the purpose of the present paper to explore thoroughly the 
 
 relationship between multistep and Nordsieck formulas on a uniform mesh. 
 
 The results apply to a limited extent to variable order variable stepsize 
 
 algorithms which discourage order and stepsize changes, because such 
 
 algorithms often produce sequences of values computed with the same 
 
 formula and stepsize. 
 
 Consider the system y' = f(t, y) with exact solution y(t). 
 
 Given starting values y., y! , j = 0(l)k-l, a linear k-step formula on a 
 
 mesh t := t_ + nh, n = 1(1)N, determines approximations {y } to the values 
 n n 
 
 {y(t )} by means of p(E)y . = ha(E)y' . and y' = f(t , y ) where 
 n n-k n-k n n n 
 
 k k 
 
 p(£) := Z a. C J , a(£):= £ 6. ? 3 , 
 j=0 J j=0 J 
 
 2 2 
 and E is the shift operator. It is assumed that ex ^ and a, +3, > 0. 
 
 k k 
 
 For convenience a normalization such as a = 1 or Z$ = 1 is not imposed. 
 If 3 n $ 0, the formula is implicit and requires the solution of a system 
 of nonlinear equations, which is always possible for small enough h if 
 f(t, y) is Lipschitz continuous as a function of y. In §2 a simple 
 reformulation of the general k-step formula is given which requires that 
 only k values be saved between steps. In practice, k+m values might be 
 stored, from which it is possible to construct a predictor formula of 
 order k+m-1 except for pathological cases where p(£) and a(£) have a 
 common factor. 
 
In §3 linear multistep formulas are described in terms of a 
 
 local polynomial approximation p (t) which interpolates the k-hn values 
 
 which are retained after completion of the n-th step. It is from this 
 
 basis free description that equivalent forms of multistep formulas are 
 
 derived. One convenient representation of the polynomial p (t) is in 
 
 n 
 
 terms of the scaled derivatives h p J (t )/j!, j = 0(l)k+m-l. Multistep 
 
 formulas represented in this way belong to the class of formulas 
 
 introduced by Nordsieck. Another is in terms of ordinates p (t .), 
 
 n n-j 
 
 i = 0(l)k -1, and derivatives hp'(t .), j = 0(l)k o -l, where k n + k. = 
 1 n n-j 2 12 
 
 k+m, which yields the modified multistep formulas of Gear [1971a, p. 150] 
 
 Backward differences of ordinates or derivatives is also possible. 
 
 General linear Nordsieck formulas compute approximations 
 
 a =: [y , hy',..., h y /q!l to the vector of scaled derivatives 
 ~n n n n 
 
 [y(t ), hy'(t ),..., h q y^(t )/q!] by the recurrence 
 n n n 
 
 a := Pa . + £6 
 ~n ~n-l ~ n 
 
 where 6 is chosen so that y' = f(t , y ); that is, 6 is implicitly 
 n n n n n 
 
 defined by 
 
 e.Pa , + JL6 = hf(t , e_Pa . + l n 6 ) . 
 ~1 ~n-l In n ~0 ~n-l n 
 
 Here P is the Pascal triangle matrix defined by 
 
 P. :- (.h, 0<i<j<q, 
 
 T T 
 
 e. is the i-th unit vector, and I =: [£., £,,...,£] is a vector of 
 ~j 1 q 
 
 parameters which characterize the formula. It is assumed that £ ^ 0. 
 Note that zero-origin indexing is being used. Nordsieck [1962] derives 
 particular choices for % which optimize stability and accuracy for 
 q=2, 3, 4, 5, 6, and he implements the formula for q = 5 as a variable- 
 step computer program for the Illiac I. By direct calculation he shows 
 
that there is a connection between this formula and the sixth order 
 
 Adams-Moulton formula; specifically, it is shown that the zero-th 
 
 and first component of a , a ,,..., a c satisfy the difference equation 
 r ~n ~n-l ~n-5 
 
 of the Adams-Moulton formula. 
 
 Further understanding of Nordsieck formulas resulted from 
 
 the work of Descloux [1963]. An actual equivalence is shown between 
 
 Nordsieck' s (q+l)-value formula and the (q+l)-st order Adams-Moulton 
 
 formula, in the sense that a linear transformation was constructed 
 
 relating the Nordsieck vector to the vector of values [y , hy',..., 
 
 n n 
 
 T 
 hy' , , ] . An equivalence was also shown for the q-th order Adams-Moulton 
 n-q+1 
 
 formula and a (q+l)-value Nordsieck formula with a slightly different 
 £ vector. This second family of formulas is used in the nonstiff option 
 
 of DIFSUB (Gear [1971b]) and its descendants. Furthermore, Descloux 
 
 T 
 showed that in both cases the predicted value y _:= e,_Pa , satisfies 
 
 n,U ~U ~n-l 
 
 the difference equation of the k-th order Adams-Bashf orth formula. 
 
 General k-step formulas were also considered and equivalent Nordsieck 
 
 formulas were constructed having the more general form 
 
 (1.1) a = Aa _ + £6 
 
 ~n ~n-l ~ n 
 
 where the dimensionality of the arrays might be as high as 2k and where A 
 
 is not necessarily the Pascal triangle matrix. In the appendix of Descloux' s 
 
 report is a theorem stating that the values y and f computed by formulas 
 
 n n 
 
 even more general than (1.1) satisfy the difference equation of an 
 associated linear multistep formula. The proof is not given by Descloux 
 because it "is rather painful and long." However, in a later paper of 
 Osborne [1966] a similar result is stated and a short and elegant proof 
 is given. 
 
In §4 a one-to-one correspondence is established between (i) 
 the class of all (q+l)-value linear Nordsieck formulas and (ii) the 
 class of all linear multistep formulas of formal order at least q 
 and stepnumber at most q. This correspondence is given in closed form 
 by 
 
 ^"VO - det(Sl-P) Jai-P)" 1 *, 
 
 and 
 
 C m_1 a(5) = det(£l-P) eJcCl-P)" 1 ^ . 
 
 Osborne [1966] gives an expression for E, p(£)/(£ - l) similar to that 
 
 given above and proves a correspondence between Nordsieck and multistep 
 
 formulas of order at least q+1 (although it is not clear that I is 
 
 correctly chosen in the proof). Presumably at that time only formulas 
 
 of optimal order were considered to be of practical value. More 
 
 recently Wallace and Gupta [1973] obtain expressions for E, p(£) and 
 
 E, o(E,) which are equivalent to those given above. 
 
 In §5 it is shown that the zero-th and first components of 
 
 a , a ..,..., a . satisfy the difference equation of the corresponding 
 ~n ~n-l ~n-k 
 
 multistep formulas. Also, if p(£) and o(E,) have no common factor, then 
 
 the Nordsieck formula is shown to be equivalent in the sense of Gear 
 
 [1971a, p. 143] to the corresponding multistep formula, which means that 
 
 a is a linear transformation of certain linear combinations of the 
 ~n 
 
 values y . , y 1 ., i = 0(l)k-l. It is proved that there is no equivalence 
 y n-j •'n-j 
 
 to a multistep formula if p(?) and a (£) have a common factor. 
 
 In §6 predictor-corrector Nordsieck formulas are considered. 
 An equivalence to predictor-corrector multistep formulas is shown to 
 be the case for P(EC)* formulas but not for P(EC)*E formulas. 
 
 In §7 recent applications of equivalence results are discussed. 
 
2. Minimum Storage for Multistep Formulas. The computation of each solution 
 
 value y and its derivative y' by a general linear multistep formulas would 
 ^n "n 
 
 seem to require saving 2k previously computed values y ,,y „,..., y . 
 
 n-i n-z n-k 
 
 and y' .. , y' _,..., y' . . For this reason it has been considered desirable 
 n— i n— 2 n-k 
 
 to have formulas with most of either the a or 3 coefficients set to zero 
 (see, for example, the "minimum storage methods" of Dill and Gear [1971]). 
 However, there is an easy way to get by with just k values without re- 
 evaluating the function f. Define 
 
 s J := I (-a. ...y . + hR ...y' .) 
 n . n k-i+1 n-i k-j+1 n-i 
 
 i=0 J 
 
 for j = 0(l)k-l. In terms of the partial sums s , the linear multistep 
 formula becomes 
 
 (2.1) solve ay = hg f(t , y ) + s " for y , 
 U n (J n n n-i n 
 
 set y' = f(t , y ) , 
 n n n 
 
 S r, = " a iy„ + h3 l^ + S „ 1 ' 
 
 n In In n-1 
 
 s = -a n y + h3 y 1 + s , 
 n k-1 n k-1 n n-1 
 
 s = -a. y + h8 y ' . 
 n k n k n 
 
 Thus only the k partial sums s -,,s .,,..., s , need to be saved in order 
 
 n-1 n-1 n-1 
 
 to advance the computation. Of course, this is not the only set of values 
 
 that could be saved; any nonsingular linear transformation of them could 
 
 be used instead. 
 
 In practice the implicit equation (2.1) is approximately solved 
 
 by some iterative process. An initial guess y is obtained with a 
 
 n, 
 
 predictor formula 
 
k* k* 
 
 a*y + I a*y . = h Z 3*y' 
 n,0 j'n-j j'n-j 
 
 and subsequent iterates y are given by an iteration of the form 
 
 n,u 
 
 y , := y -[a n l - h3„J ] _1 [a„y - hg rt f - s k "^] , 
 ^n,y+l 'n,y n,y l (Tn,y n,y n-1 ' 
 
 y = 0(l)M(n)-l , 
 
 where f := f(t ,y ) , 1 is the identity matrix and J is a matrix 
 n,y n n,y n,y 
 
 depending on the type of corrector iteration. For functional iteration 
 
 J := 0, for the Newton iteration J := f (t y ), and for the 
 n,y n,y y n/n,y" 
 
 chord iteration J := f (t ,y .,) . The number of iterations M(n) 
 n,y y n ^n,0 J 
 
 might vary from step to step. After determining y := y w , . a final 
 
 n n,M(n) 
 
 function evaluation may or may not be performed. For a P(EC)* formula 
 
 y' := f ... N , is accepted and for a P(EC)*E formula y' := f = f w/ .. 
 ^n n,M(n)-l r J n n n,M(n) 
 
 For our purposes the essential difference between the two types of 
 
 predictor-corrector schemes is that the values y and y' satisfy the 
 
 n n 
 
 corrector formula for the P(EC)* scheme and they satisfy the differential 
 
 equation for the P(EC)*E scheme. In fact, one could view the P(EC)* 
 
 formula as computing values y , y' that satisfy a slightly different 
 
 n n 
 
 differential equation, which implies that equivalence results for linear 
 formulas apply equally well to P(EC)* formulas. 
 
 How many values must be saved between steps for these predictor- 
 corrector formulas? The answer would seem to be the row rank of the 
 following matrix: 
 
a T 
 
 a 
 
 K-l 
 
 a. 
 
 a* 
 K 
 
 a* 
 K-l 
 
 a. 
 
 a. 
 
 °K 
 
 a 
 
 K 
 
 K 
 
 K-l 
 
 P K 
 
 6* 
 P K-1 
 
 K 
 
 ft* 
 P K 
 
 K 
 
 a* 
 
 a* 
 
 • "a*. 
 
 3* 
 
 RA 
 
 • 3* 
 
 1 
 
 2 
 
 K 
 
 1 
 
 2 * 
 
 K 
 
 where K := max{k, k*}. The row rank of this matrix is at least K and 
 could be as high as k+k*. For reasons of storage economy the predictor 
 should be chosen so that the rank of the matrix is low. The most 
 economical choice would be to use a predictor of the form 
 
 1 L J k-l 
 
 V0 = Vn-l + Vn-l + - +T k-l S n-l 
 
 where the coefficients Y. are chosen to make the formula of order k-l 
 or greater. 
 
 A predictor of order k-l may not be accurate enough for a 
 corrector of order k or greater. For this reason it may be preferable 
 to store more than k values between steps. For a predictor of order at 
 least q where q >_ k-l, one would expect to need an additional m = q-(k-l) 
 values. This can be done by including information from preceding meshpoints 
 Since each step of the computation introduces exactly one new item of 
 information, namely y 1 , the following set of k+m values is suggested: 
 
 n-m 
 
 j = 0(l)k-l , 
 
 hy* . , j = 0(l)m-l . 
 
 From these values one can apply the corrector formula to generate 
 
y .. , 8 ,.,..., s ,.» j " l(l)m. Moreover, if the corrector 
 n-m+j n-mf j n-m+j 
 
 formula is of order q or more, then the polynomial of degree q or less 
 which interpolates the original set of k+m values also interpolates 
 the values generated from them. An alternative set of values from 
 which the others can be recovered (by reversing the forward step 
 procedure) is the following: 
 
 s j , j = 0(l)k-m-l* , 
 
 n-m 
 
 y n-j' j = °*( 1 ) m - 1 ' 
 
 hy' ., j = 0(l)m-l , 
 n-j 
 
 where 0* - 0, 1* = 1 if m < k and 0* = 1, 1* = if m = k+1. This set 
 
 has the important advantage of requiring values from only k+0* meshpoints 
 
 instead of k+m meshpoints, although it is a little more awkward to specify, 
 
 The case m = k+1 can only apply to the trapezoid formula, the Milne 
 
 formula, and other (unstable) k-step formulas of maximal order 2k. For 
 
 this case the set of values y ., i = l(l)k,and hy' ., i = 0(l)k, is used 
 
 n-j J n-j 
 
 to define a polynomial of degree at most 2k. This polynomial interpolates 
 
 y because y is determined from the other values by a corrector formula 
 n n 
 
 which is exact for polynomials of degree at most 2k. 
 
 As an example, a k-th order predictor for the (k+l)-st order 
 Adams-Moulton formula would be constructed from the values 
 
 j 
 
 . _ k-i+i n-l-i n n 
 
 i=0 J 
 
 or equivalently linear combinations of these values: y , and hy' ., 
 
 n J n-2 
 
 j = 0(l)k-l. The predictor so constructed is the k-th order Adams- 
 Bashforth formula. 
 
It is not obvious, however, that such a predictor can always 
 
 be constructed. Clearly the question is of practical interest only for 
 
 implicit formulas (p, a) of order at least q. It is claimed that the 
 
 construction of a suitable predictor is possible if and only if a unique 
 
 polynomial p(t) of degree q or less is uniquely determined from the 
 
 values L^p(t ), j = 0(l)k-l, and p'(t . ), j = 0(l)m-l, where the 
 h n-m n-j 
 
 operators L? are defined by 
 
 i j 
 L^z(t) = Z (-a. .^z(t-ih) + h6, .,.z'(t-ih)) . 
 i=0 k_ J +:L k-j+i 
 
 Clearly, the possibility of constructing p(t) implies the existence of 
 
 a suitable predictor. On the other hand, if p(t ,,) can be determined 
 
 n+1 
 
 by a predictor, then from the corrector formula we can get p'(t . ) . 
 
 This can be repeated to generate p(t l0 ), p(t ,_),... until there are 
 
 n+z n+3 
 
 enough values to construct p(t) by interpolation of ordinate values. 
 Necessary and sufficient conditions for the possible construction of 
 p(t) is given by the theorem that follows. First a precise definition 
 of order is needed. 
 
 Definition. A linear multistep formula (p, a) is of order 
 at least q if 
 
 Mi±f} . io g (i +z) + o (z s +1 ) 
 
 and is of formal order at least q if 
 
 p(l+z) = log(l+z) a(l+z) + 0(z q+1 ) 
 
 (cf. Gear [1971a, p. 119]). 
 
 THEOREM 2.1. Let (p, a) be a linear k-step formula of order 
 
 at least q. Values L^p(t ), j = 0(l)k-l, and p'(t .), j = 0(l)m-l, 
 
 h n-m n-j 
 
 uniquely determine a polynomial p(t) of degree at most q if and only 
 if p(£) and a(^) have no common factors. 
 
10 
 
 Proof. If polynomials are identified with the column vector 
 of their scaled derivatives, then the unique existence of p(t) is 
 
 equivalent to the linear independence of the rows of the matrix 
 
 i n T -k T ^T T T -1 T n l-nK . 
 T := col{A , A ,..., ^ » e , e P ,..., e^V ) where 
 
 ,T ^ r T , Ti_— i— m 
 A. := Z {-a e + 3, ...e.JP 
 ~J ±=0 k-j+i~0 k-j+i~l 
 
 and P is the (q+1) * (q+1) Pascal triangle matrix. The row vectors 
 
 T 
 A , j = 0(l)k-l, have a simple generating function, for 
 
 (-p(S)eJ + a(?)e^)(I-?P" 1 )"V m 
 
 = Z C k " j (-a.e^ + B.e*) Z 6" 1 " 1 
 j=0 J ~ U J i=0 
 
 k " 1 £ T 
 £=0 
 
 T k-£ T - 
 
 where we have used the fact that A P = . Let p(£) = tt(£) p(£) and 
 
 a(£) = tt(D a(0 where tt(^) is the g.c.d. of p(£) and a(^). Then 
 
 — — — T— T— 1— 1 — m 
 
 (p, a) is of order q or greater, and so (- p(^)e + o(E,)e ) (1-^,? ) P 
 
 -T 
 is a generating function for the y.» which are defined in the obvious 
 
 way. The relation 
 
 (-P(C)eJ + cr(DeJ)(I-£P~ 1 rV" ,n 
 
 = Tr(0(-P(0ej + a(?)eJ)(I-?P" 1 )" 1 P -m 
 
 implies that the first k rows of T are linear combinations of the row 
 
 -T 
 vectors y.» J = 0(l)k-l. If p(£) and a(£) have a common factor so that 
 
 k < k, then p(t) is obviously not uniquely determined. Assume that 
 
 p(£) and a(£) have no common factors. Suppose that the rows of T are 
 
 linearly dependent. Then Tv = for some v ^ 0. Hence 
 
11 
 
 (-p(Oel + c(Oe*)(I-£P X ) X P "v = 
 
 or equivalently 
 
 (2.2) P(C)t q (?) = a(?)T 1 (C) 
 
 where T (£) = det ( (I-?P~ 1 )P in )eT(I-?P~ 1 )~ P~ m v, which by Cramer's rule is 
 simply the determinant of (I-£P )P with the i-th column replaced by v. 
 Thus T . (O is a nonzero polynomial of degree q or less. Furthermore 
 
 ^(o = -a-^+V^u-rVV^v 
 
 oo 
 
 = (1-C) £ 5 e ,P v 
 
 j-o 
 
 = (-l)^C e Pv + lower order terms 
 
 and so t.. (?) is of degree at most k-1. Hence from (2.2) at least one of 
 
 the k roots of p(?) must be a root of G(£). This is a contradiction and 
 
 therefore the rows of T are linearly independent. □ 
 
 COROLLARY. The values L^p(t ), j = 0(l)k-m-l*, p(t .), 
 
 n n-m n-j 
 
 j = 0*(l)m-l, and p'(t _.), j = 0(l)m-l, uniquely determine p(t) if and 
 only if p(£) and a(?) have no common factor. 
 
 Remark 1. It is not enough to have formal order of at least 
 
 2 
 q. Consider p(£) = (£-1) and 0(O = 0; the values -p(t ) and 2p(t ) - 
 
 n n 
 
 p(t ., ) uniquely determine a polynomial of degree one or less. 
 n-1 
 
 Remark 2. Theorem 2.1 appears to be related to a result of 
 Dahlquist [1975] concerning the equivalence between a linear k-step 
 formula and the corresponding one-leg k-step formula under the assumption 
 that p(?) and a(£) have no common factor. 
 
 Apart from the need for a predictor there is another complication 
 that affects the amount of storage needed. In practice it has been found 
 
12 
 
 desirable to vary the stepsize, and thus a variable step form of the 
 multistep formula must be used. For any given fixed step formula there 
 are numerous variable step formulas, some of which require less storage 
 than others. This topic will be studied in a future paper. 
 
13 
 
 3. Construction of Linear Nordsieck Formulas. In this section linear 
 multistep formulas are described in terms of local polynomial approximations, 
 and an equivalent Nordsieck formula is constructed by representing the 
 polynomials in terms of their coordinates with respect to local monomial 
 bases. 
 
 Consider a linear k-step formula of order at least q where 
 q >_ k such that p(£) and G (£) have no common factor. (The case q = k-1 
 is treated at the end of this section.) Then as a consequence of the 
 corollary to Theorem 2.1 a unique polynomial of degree at most q is 
 determined by the q+1 conditions 
 
 L k-l (t n-l-J " S n-1- B - i - «»*-l'. 
 P„-l (t „-l-j )= " y n-l-j' 3-0*(Dm-l, 
 K-lK-l-S = K-i-y J-0(D-1. 
 
 The polynomial p , (t) can be regarded as an approximation to the solution 
 n-1 
 
 y(t) near t = t , and it can be used to obtain a predicted value 
 n— 1 
 
 y _ = p ,(t ) as an initial approximation to y . In terms of the values 
 n , n- In n 
 
 y« i y„ ir» and y' ,..., y' , the relation y = p (t ) becomes 
 
 n-1 n-k n-1 n-k n,0 n-1 n 
 
 an explicit linear k-step formula. 
 
 Advancing the numerical solution one step yields values which 
 
 determine the polynomial p (t). However, there is a more direct way of 
 
 n 
 
 expressing p (t) in terms of p ,(t). First, we examine how y and y 1 
 n n-1 n n 
 
 are determined. From the multistep formula we have 
 
 ay = h3 n y' + s ~ 
 On On n-1 
 
 and because L. annihilates p , (t) we have 
 h n-1 
 
14 
 
 a ny« n = he ny' n + L u~ P i < c i> 
 n,0 n,0 h n-1 n-1 
 
 where y' _ = p' 1 (t ) . It follows from the defining conditions that 
 n,U n-1 n 
 
 k-1 k-1 
 
 K p„ i (t « i^ = s „ T Hence a n(y n ~y n r) = he n ( yJ, _ yJ, J- Also 
 
 n n-i n-l n-1 U n n,U U n n,U 
 
 y' = f(t , y ). Together these last two equations define y and y'. 
 n n n n n 
 
 A more convenient way of expressing this is to introduce an increment 
 
 6 which satisfies 
 n 
 
 hy' n + a_6 = hf (t , y + 6 ) 
 n,0 On n n,0 n 
 
 and set 
 
 y n = y n,0 + 3 6 n ' 
 
 hy' = hy' + a_6 . 
 n n,0 n 
 
 Second, we construct d (t) := p (t) - p ,(t). For j = 0(l)k-m-l* 
 
 n n n-1 
 
 L h d n (t n J = L hPn (t n J " { ~\ ^r, 1 (t n J + h3 V iK 1 (t « J 
 
 h n n-m h n n-m k-j n-1 n-m k-j n-1 n-m 
 
 + Lrj" 1 ? .(t )} 
 
 h n-1 n-l-m 
 
 = s j - {-a. .y + hR .y' + s j "^ } 
 n-m k-j n-m k-j n-m n-l-m 
 
 = ; 
 
 for j = 0*(l)m-l 
 
 jVn if 3 = , 
 d n (t n _.) =< 
 
 3 1 otherwise; 
 
 and for j = 0(l)m-l 
 
 Hence 
 
 jVn if J = ° ' 
 
 hd n (t n- J = \ 
 
 I otherwise. 
 
 t-t 
 
 d (t) = A(-r^)6 
 n h n 
 
15 
 
 where A(x) is the unique polynomial of degree at most q which satisfies 
 
 L^A(-m) - , j - 0(l)k-m-l* , 
 
 3 , j = 0*(1)0 , 
 j = l(l)m-l , 
 
 I a n • n 
 
 A'(-j) -< ° « J = ° ' 
 
 \J) , j = l(l)m-l . 
 
 Briefly, we have shown that the numerical solution is advanced one 
 
 step by setting 
 
 t-t 
 
 (3.1) p (t) = p (t) + A (-t-^)6 
 
 n n-l n n 
 
 where 6 is chosen so that p (t) satisfies the differential equation and 
 n n 
 
 A(x) is characteristic of the multistep formula. 
 
 The polynomial formulation of the k-th and (k+l)-st order 
 Adams-Moulton formulas was discovered by Descloux [1963]. Schemes based 
 on general choices of A(x) are discussed by Skeel [1973] and by Wallace 
 and Gupta [1973], who give an interesting interpretation of polynomial 
 schemes in terms of polynomial predictive filters. They derive new 
 formulas for stiff problems by choosing A(x) to be a monic polynomial 
 which best approximates zero for x < 0. Different types of 
 approximations yield different formulas. Still more formulas are 
 presented in Gupta and Wallace [1975] and Gupta [1976]. The 1975 paper 
 uses the term modifier polynomial for A(x). (In the 1973 paper this 
 term is applied to A(x/h).) 
 
 A very simple identification of these polynomial schemes with 
 
 Nordsieck formulas becomes apparent if the polynomials are represented 
 
 by vectors of scaled derivatives. Let a = [p (t ), hp' (t ),..., 
 
 ~n n n n n 
 
 h q p (q) (t )/q!] T and I = [A(0) , A'(0),..., A (q) (0) /q ! ] T . Then from 
 n n ~ 
 
 (3.1) we have 
 
16 
 
 p (j) (t ) = P (j 1 ) (t n ) + h" j A (j) (0)6 , 
 n n n-± n n 
 
 from which it easily follows that 
 
 a = Pa , + 16 
 ~n ~n-l ~ n 
 
 -IT T 
 
 with 6 chosen so that h e,a = f(t , e.a ). 
 n ~l~n n ~0~n 
 
 For the k-th order backward differentiation formula 
 
 k 
 
 Z a.y , = hy' 
 j-0 3 n " J 
 
 the modifier polynomial A(x) must satisfy 
 
 A(-j) = , j = l(l)k-l , 
 A(0) = 1, A'(0) = a Q , 
 
 whence 
 
 k-1 
 A(-k) = (A'(0) - Z a.A(-j)}/a, 
 
 j=0 J 
 
 = . 
 
 Therefore 
 
 A(x) = ( X ^ k ) 
 
 3=0 k! ^ +1 
 
 where the (square) brackets denote Stirling numbers of the first kind 
 
 (see, for example Knuth [1968, p. 66]). The fact that A(x) vanishes 
 
 for x = -1 , -2 , . . . , -k means that the values y , , y „,...,y , are 
 
 n-1 n-z n-k. 
 
 "remembered" after advancing from t , to t . 
 
 n-1 n 
 
 For the (k+l)-st order Adams-Moulton formula 
 
17 
 
 y -y = h I B.y' . 
 n n_1 j=o J n ~J 
 
 the modifier polynomial of the (k+2)-value form must satisfy 
 
 A'(-j) = , j = 2(l)k-l , 
 
 A(0) = 3 Q , A»(0) = 1 , 
 
 whence 
 
 A(-l) = A'(-l) = , 
 
 k-1 
 
 A'(-k) = {A(0) - A(-l) - Z 3. A , (-j)}/3 1 
 
 J-0 J 
 
 = 
 
 Therefore 
 
 and so 
 
 A»(x) = ( X + k ) 
 
 A(x) - / ( 7 t k )dy 
 -1 k 
 
 = -/( y >y + / ( y > y . 
 
 
 
 In terms of Stirling numbers 
 
 / C k )dy = Z TTTT [ , ]x J . 
 
 o k j=i Jk ' J 
 
 These are the formulas used by the nonstiff options of the codes DIFSUB 
 of Gear [1971], GEAR Rev. 3 of Hindmarsh [1974], and EPISODE of Byrne 
 and Hindmarsh [1975] when the stepsize does not vary. 
 
 The (k+l)-st order Adams-Moulton formula could also be written 
 as a (k+l)-value Nordsieck formula. The modifier polynomial for such a 
 formula can be determined by applying the theorem that follows. 
 
18 
 
 THEOREM 3.1. Let A(x) be the modifier polynomial of the (q+1)- 
 value form of a linear k-step formula of order at least q where 
 q > k+1 and p(£) and a(£) have no common factors. Then A(x) := A(x) - A(x-l) 
 
 is the modifier polynomial for the q-value form of the multistep 
 formula. 
 
 Proof. For m <^ k 
 
 L m A(-m) = L A(-m) - Z (-a.A(-i) + 0.A'(-i)} 
 
 1=0 1 X 
 
 = , 
 and for m = k+1 
 
 3 1 k 
 A(0) = -f- A' (0) + -f- Z {-a.A(-i) + a.A'(-i)} 
 
 a o a o i=l 1 X 
 = V 
 
 Hence 
 
 L:JA(-m) = , j = 0(l)k-m , 
 
 A(0) = 3 Q , A'(0) = a Q , 
 
 A (-j) = A'(-j) =0, j = l(l)m-l , 
 from which the theorem follows. D 
 
 Therefore the (k+l)-value form of the (k+l)-st order Adams- 
 Moulton formula has 
 
 A(x) = / ( y + k )dy - f (^ k )dy 
 -1 R -1 k 
 
 = / ( y f)dy + / {(^ k ) - ( y+ ^ 1 )}dy 
 -1 K 
 
 ■ - "' C^ + { ?Z>y ■ 
 
19 
 
 These are the original formulas of Nordsieck. 
 
 Let us consider now the case m=0. The values <5 , y , and y' 
 
 n n n 
 
 are determined as before. However, for j = 0(l)k-l 
 
 I/}d (t ) = -a, .(y -y _) + hB, . (y'-y' n ) 
 h n n k-j n n,0 k-j n n,0 
 
 = (-a, .B n + 8, .cx.)6 . 
 k-j k-j n 
 
 Therefore the solution is advanced one step as given by (3.1) with A(x) 
 
 determined by 
 
 L^A(O) =-a k _. 3 Q + 3 k _.a Q , j = 0(l)k-l. 
 
 However, p (t) does not necessarily interpolate y and y 1 nor does 
 n n n 
 
 it generally satisfy the differential equation at t = t . Therefore, 
 
 in the scaled derivative form this scheme is a generalized Nordsieck 
 
 formula in the sense that 6 is determined by a condition other than the 
 
 n 
 
 satisfaction of the differential equation at t = t . Such formulas are 
 
 n 
 
 potentially useful because of their minimum storage property. Nordsieck 
 
 [1962, p. 27] considers the possibility of having Z = $ but JL f a 
 
 so that p (t) interpolates y but not y' ; however, the results of his 
 r n n n 
 
 experiments were not promising. Also, Wallace and Gupta [1973] mention 
 that "Other choices of 6 are rationally possible, and we hope to explore 
 some other choices later." Finally, it is worth noting that Theorem 3.1 
 extends to the case q = k. 
 
20 
 
 4. The Correspondence between Multistep and Nordsieck Formulas. In the 
 preceding section it was shown how to construct the modifier polynomial 
 A(x) of a Nordsieck formula from the polynomials p(£) and o(E,) of a linear 
 multistep formula provided that p(£) and a(£) have no common divisors. 
 In this section we show how to obtain p(£) and o(E,) from A(x), and we 
 establish a one-to-one correspondence between (i) the class of all (q+1)- 
 value linear Nordsieck formulas and (ii) the class of all linear multistep 
 formulas of formal order at least q and of stepnumber at most q (including 
 those for which p(£) and o(E,) have common factors). 
 
 For each linear Nordsieck formula we define a corresponding 
 linear multistep formula (p, a) by 
 
 (4.1) 
 
 p(0 := det(Cl-P)e] , (?I-P) h, 
 
 5(0 := det(£I-P)eJ(£I-P) 1 %. 
 
 T -1 
 Applying Cramer's rule to e.(^I-P) £ for j = 1 and j = yields 
 
 P(£) = det 
 
 5-1 
 
 o 
 
 
 
 -1 
 
 -2 
 
 -1 
 
 and 
 
21 
 
 5(0 - det 
 
 % 
 
 -l 
 
 -1 
 
 -l 
 
 *1 
 
 ?-i 
 
 -2 
 
 -2 
 
 h 
 
 
 
 5-1 
 
 -( q ) 
 
 • 
 • 
 • 
 
 • 
 • 
 • 
 
 
 
 £ 
 . q 
 
 
 
 
 
 5-1 
 
 Hence 5(C) is a polynomial of degree q or less, and p(C) is of degree q 
 since £ ^ 0. Let 
 
 p(5) =: Z a.? q 1 and 5(C) =: Z 3.C q 1 
 
 i=0 
 
 i=0 
 
 It might happen that a =3 =0, and hence express 
 
 q q 
 
 (4.2) 
 
 p(5) =: C m 1 P(?) and 5(C) =: C™ Vo 
 
 2 2 
 where a, + 3, > and m-1 = q-k. 
 k k 
 
 The following theorem shows that the p(C) and a(C) corresponding 
 to A(x) can be used to reconstruct A(x) by the process of §3 if it is 
 applicable. 
 
 THEOREM 4.1. With P(C) and a (5) derived from A(x) by (4.1) and 
 (4.2) we have 
 
 L^A(-m) =0,j= 0(l)k-l , 
 
 £ Q = A(0) = 3 Q , % 1 = A'(0) = a Q 
 
 A(-j) = A'(-j) = , j = l(l)m-l , 
 
 A(-m) = (-D q 3 k , A'(-m) = (-1)\ . 
 Proof. The relations £ = 3 n and £ = a follow from the 
 
 expressions for p(C) and a(C) as determinants. From (4.1) we have 
 
22 
 
 P(U = -(?-l) q+1 e^(I-?P X ) h X £ 
 
 oo 
 
 = -(?-D q+1 z Jr^p-j-^ 
 
 j=o ~ x 
 
 CXD 
 
 = -(?-l) q+1 Z A'(-j-l)? j . 
 J-0 
 
 Clearly, A'(-j) = 0, j = l(l)m-l, a k = (-l) q A'(-m), and 
 
 00 
 
 p(0 = -(?-l) q+1 Z A'(-j-m)C J . 
 j=0 
 
 Similarly A(-j) = 0, j = l(l)m-l, 3 R = (-l) q A(-m), and 
 
 oo 
 
 a(C) =-(C-D q+1 2 A(-j-m)? j . 
 3=0 
 
 From these expressions for p(£) and a(£) we get 
 
 OO 00 
 
 -p(0 Z A(-j-m)? j + o(0 Z A'(-j-m)? j = . 
 3=0 j-0 
 
 Equating coefficients of the powers of £ completes the proof .D 
 
 The next theorem shows that the multistep formula (p, 5) 
 
 corresponding to A(x) is of order at least q, and therefore it can be used 
 
 to reconstruct A(x) by the process of §3 if p(£) and a(£) have no common 
 
 factors. 
 
 THEOREM 4.2. The linear multistep formula corresponding to 
 
 a (q+1) -value linear Nordsieck formula has formal order at least q. 
 
 T 
 Moreover j it is of order at least q+1 if and only if b I = where the 
 
 components of b are the Bernoulli numbers defined by 
 
 00 b.z J 
 e Z -l j-0 J ' 
 
23 
 
 Proof. According to the definition of formal order and 
 equation (4.1) it must be shown that 
 
 z q+1 e[((l+z)I-P) _1 £ = log(l+z)z q+1 eQ((l+z)I-P)" 1 £ + 0(z q+1 ) , 
 
 and so it is enough to show that 
 
 e^(I-z _1 F) -1 = log(l-fz)eJ(I-z" 1 F)~ 1 + 0(z) 
 
 where F := P-I. Because F " = 0, it is not difficult to verify that 
 
 log(l+z)(I-z _1 F) -1 = log(I+F)(I-z~ 1 F)~ 1 + 0(z) , 
 
 T T 
 
 and so it suffices to show that e_ log(I+F) = e.. . This holds because 
 
 I+F = P = exp(Q) where 
 
 Q = 
 
 o 
 
 o 
 
 o 
 
 To prove the second assertion, retrace the first two steps of the proof 
 to get the following condition for formal order of at least q+1: 
 
 eJ{l-|F + ...+ ^ r (-F)lH = . 
 
 The matrix in braces is 
 
 log (i+F) = Q = 2 j_ n J 
 
 F ' exp(Q)-l ± 1 V- ' 
 
 from which the second part of the theorem follows. □ 
 
 Remark 1. The matrices P, F, and Q obviously represent the 
 linear operators shift, forward difference, and derivative, respectively, 
 for polynomials of degree q or less. Defining B = I-P to be the matrix 
 
24 
 
 representing the backward difference operator V it is not difficult to 
 
 T 
 show that the condition b Z = is equivalent to 
 
 q 1 1 
 
 Z ~r V J A(-1) = 
 
 j=0 JT± 
 
 as well as 
 
 q-1 , 
 A(-l) - Z y^. 1 V J A , (-1) = 
 j=0 3+i 
 
 where the y*., are the coefficients for the backward difference form 
 
 of the Adams-Moulton formulas. These conditions on the polynomial A(x) 
 
 arise in another situation in Henrici [1962, p. 342]. 
 
 T 
 Remark 2. The condition b I is obtained by Wallace and 
 
 Gupta [1973] although the coefficients b. are not identified as the 
 
 Bernoulli numbers. Instead a recursive definition is given for the 
 
 b . , which should read 
 
 j 1+1 
 b = 1, Z ( 3 7)b. = . 
 
 i=0 
 
 Remark 3. It is shown by Gear [1966] that a stable linear 
 
 Nordsieck method is convergent of order q in all components of a , and 
 
 ~n 
 
 it is shown by Skeel and Jackson [1977] that it is convergent of 
 
 T 
 order q+1 in all components if and only if b Z = 0. 
 
 It has thus been shown that to each linear Nordsieck formula, 
 
 which is uniquely specified in terms of £, there corresponds a linear 
 
 multistep formula of formal order at least q and of stepnumber at most q, 
 
 We now establish a correspondence in the opposite direction. We have 
 
 from (4.1) that 
 
25 
 
 5(l+z) - z q+1 eJ(zI-F) 1 £ 
 
 = z q e T I Z " j F j £ 
 J-0 
 
 where F := P-I is the matrix of the forward difference operator. Equating 
 
 coefficients of powers of z gives 
 
 2^ = A^A(O) . 
 
 and hence A(x) is determined from its corresponding 5(0 by Newton's 
 forward difference formula 
 
 q a (j) Q) x 
 (4.3) A(x) = I Q ./ ( X .) . 
 
 3=0 J ! q ^ 
 Since the linear multistep formula is uniquely specified by 5(?) , the 
 one-to-one correspondence is established. The inverse transformation 
 is conveniently expressed as 
 
 q 
 
 a(U = Z A J A(0)(C-D q J 
 j=0 
 
 and 
 
 q-1 . 
 p(?) = E A J A'(0)(^-l) q J . 
 j-0 
 
 Note from (4.3) that the popular normalization a(l) = 1 corresponds to 
 £ = 1/q!. 
 
 q 
 
 Remark. Theorem 3.1 can be generalized to any modifier 
 polynomial A(x) for which m >_ 2. Express d(£) =: ^d(^). Then 
 
 a (j) d) = a (j) d) + ja (j - 1} (i) , 
 
 from which it can be shown that 
 
26 
 
 q_1 S (j) (l) x+1 
 A(x) = I ., U ( X L ) 
 
 J-0 j! q ' J 
 
 and 
 
 q-i a (j) d) x 
 
 A(x) := A(x) - A(x-l) = I . ) L) ( * ,) . 
 
 j-o j! q_1 - J 
 
 Clearly A(x) is the modifier polynomial for the q-value formula, 
 
27 
 
 5. Equivalence of Linear Nordsieck Formulas to Linear Multistep Formulas, 
 Recall that a linear Nordsieck formula determines successive values by 
 the system of difference equations 
 
 (5.1) a = Pa . + £6 
 
 ~n ~n-l ~ n 
 
 where 6 is chosen so that y' = f(t , y ). We do not consider the more 
 n n n n 
 
 general formula a = Aa , + Z6 because formulas with A 4- P are of 
 ~n ~n-l ~ n 
 
 dubious value, and in any case, most of the results for A = P generalize 
 
 if minor restrictions are imposed on A. 
 
 For theoretical purposes a rewriting of (5.1) is often useful. 
 
 Premultiplying (5.1) by e? yields 6 = £~ 1 (hf -e^Pa ,). Let I := C^l 
 
 J ~1 n 1 n~l ~n-l ~ 1 ~ 
 
 ~ x 
 and S := (I-£e )P, and we get 
 
 (5.2) a = Sa , + h£f 
 
 ~n ~n-l ~ n 
 
 T 
 where f is chosen so that f = f(t , e.Sa - + h£_f ). Expressions for 
 n n n ~1 ~n-l On 
 
 p(£) and 5(5) in terms of S are given by the following theorem, whose 
 proof uses an idea from Osborne [1966, equation (4.3)]. 
 
 THEOREM 5.1. The polynomials p(£) and 5(£) defined by (4.1) 
 
 satisfy 
 
 and 
 
 p(£) = det(?I-S)e^(?I-S) _1 £ 
 
 5(0 = det(€l-S)eJ(£I-S) ^ . 
 
 Proof. We have 
 
 £l-S = £I-P + Ze^P 
 
 = (?I-P)(I + (Cl-P) ^P) , 
 
28 
 
 and so 
 
 det(£l-S)'(5l-S) -1 £ = det(Cl-P) det (I+C^I-P)" 1 ^?) 
 
 •(I+(?I-P)" 1 £e^P)" 1 (Cl-P)" 1 I . 
 
 This can be simplified by means of the identity 
 
 T T -1 
 det(I+xy )'(I+xy ) x = x , 
 
 which follows from Cramer's rule. Thus we get 
 
 det(£l-S)'(£l-S) -1 £ = det(5l-P)*(?I-P) _1 ^ , 
 from which the lemma immediately follows. □ 
 
 Note that e^S = T ; and so p(£) - det (^I-S)^" 1 ^ , which 
 gives the characteristic polynomial of S as 
 
 det(£I-S) = £~Vp(0 • 
 
 Thus the strict root condition is satisfied by the linear Nordsieck formula 
 (cf. Skeel and Jackson [1977]) if and only if it is satisfied by the 
 corresponding linear multistep formula. 
 
 The theorem that follows shows that in the case of all Nordsieck 
 formulas the zero-th and first components of the vectors a satisfy the 
 difference equation of the corresponding linear multistep formula (p, o) . 
 Hence all the limitations on the accuracy of multistep formulas (Dahlquist 
 [1956], [1963]) apply also to Nordsieck formulas. 
 
 THEOREM 5.2. The values y , y' , n - 0(1)N, computed from (5.1), 
 
 satisfy the linear multistep formula 
 
 p(E)y = ha(E)y' 
 
 n-k n-k 
 
 defined by (4.1) and (4.2) for n = k(l)N. 
 
 Proof. From (5.2) we have 
 
29 
 
 , • k-1 . 
 
 S k_:i a , = a . - Z S 1 J ilhf . 
 ~n-k ~n-i . . ~ n-i 
 i=J 
 
 and 
 
 k k-1 k-1 ._.- 
 
 p(S)a ,= Z a. a .- Z a. Z S J £hf 
 
 ~n-k . _ j~n-j s n J • j ~ n-1 
 j=0 J J j=0 i=j 
 
 T 
 Premultiplying by e_ and rearranging gives 
 
 k k-1 i _ 
 
 Z a.y . = Z el{ Z a.S J Hhf . + e„p(S)a . . 
 3=0 * °-J 1=0 ~° j-0 3 ~ n_1 ~° ~ n_k 
 
 It needs to be shown first that 
 
 q i • • 
 (5.3) a(0 = Z e^{ Z a.S 1 J H£ q X . 
 
 i=0 ~ U j-0 J 
 
 We have that 
 
 = pCOeJCl-r^)" 1 ! 
 
 = Z a.E, q 3 Z K e n S X i , 
 j=0 3 i=0 ~ ~ 
 
 which leads to (5.3). It remains to be shown that 
 
 Because S satisfies its characteristic polynomial, e p(S)S = . Also, 
 because S is of rank q, it has only one linearly independent left eigenvector 
 
 associated with the eigenvalue 0; and so e n p(S)S ' is a multiple of e.. . 
 
 T m-1 ~ T 
 This implies e p(S)S (I-let) = 0, and so by (5.3) 
 
 T ,_,,._m-l n T 
 
 ?0 P(S)S = e k+m-l!l • 
 
30 
 
 If m > 1, then e p(S)S = , and this argument can be repeated until 
 
 the assertion is established. □ 
 
 Remark 1. We have y 1 = f except possibly for n = 0. 
 
 n n 
 
 Remark 2. This theorem improves the result of Descloux [1963] 
 and Osborne [1966] in two respects: first, p(?) and cr(£) are given in 
 closed form, and second, the result is shown for n >. k rather than n >^ q. 
 
 The next theorem proves that Nordsieck and multistep methods 
 
 are equivalent in the sense of Gear, for it is shown that there is a 
 
 nonsingular matrix that relates the scaled derivative approximations a 
 
 ~n 
 
 to certain linear combinations of computed y and y' values. As a 
 consequence, in practically all cases linear Nordsieck methods cannot 
 be regarded as generalizations of linear multistep formulas (cf. Gear 
 [1971a, pp. 102, 136]). 
 
 THEOREM 5.3. Assume the polynomials p(£) and o(0 corresponding 
 to a given (q+l)-value Nordsieck formula have no common factors. Then 
 there exists a unique nonsingular (q+l)x(q+l) matrix T depending only on 
 H such that 
 
 a = T _1 y 
 n y n 
 
 for n = k-l*(l)N where 
 
 r k-m-1* , , , f -.T 
 
 y := [s ,..., s , y *,..., y , hy , . . . , hy J 
 in n-m n-m n-0 n-m+1 n n-m+1 
 
 Proof. By Theorem 4.2 the formula (p, a) is of order at least 
 
 q, and hence by the corollary to Theorem 2.1 there exists a unique 
 
 nonsingular matrix T such that 
 
31 
 
 P(t n ) 
 
 hp'(t n ) 
 
 hS (q) (t )/q! 
 
 = T 
 
 -1 
 
 LTp(t ) 
 h n-m 
 
 L"p(t ) 
 h n-m 
 
 "^Wi' 
 
 for any polynomial of degree q or less. Define values y. and y! for 
 j = l*-k(l)-l such that 
 
 The process of §3 may be applied to construct a (q+l)-value linear 
 
 Nordsieck formula which would compute vectors T y , n = 1(1)N. The 
 
 results of §4 imply that this formula would be identical to the given 
 
 Nordsieck formula, and hence a = T y . The uniqueness of T follows 
 
 ~n in 
 
 from the fact that the method is exact for all polynomials of degree 
 
 at most q if the starting values are exact. □ 
 
 Remark. For implicit formulas Wallace and Gupta [1973] give 
 
 an informal argument suggesting that the quantities a , and 6 . , 
 
 ~n-q n-j 
 
 j = l(l)q-l, can be expressed as linear combinations of y ., f ., 
 
 n-j n-j' 
 
 ■m— 1 
 
 j = l(l)q. Their conjecture is correct as long as E, p(K) and 
 
 E, O(C) have no common factors. In that case it is an immediate 
 
 consequence that the components of a are linear combinations of 
 
 y ., f ., j = l(l)q. 
 n-j n-j 
 
 Theorem 5.3 does not extend to Nordsieck formulas for which the 
 corresponding p and G polynomials have common factors. 
 
 THEOREM 4.4. Assume that p(£) and a(£) have a common factor 
 and that the formula (p, a) is of order at least q. Then a cannot be 
 
 expressed only in terms of y ., hy' ., j = 0(l)n. 
 
32 
 
 Proof. We begin by showing that there exists an eigenvector 
 
 T 
 v associated with a nonzero eigenvalue £, of S such that e~v = 0. From 
 
 (4.3) it follows that p(£.) and a(£) have the common factor 1 if and only 
 
 if £ =0. First, suppose that £ = 0. Then £ = 1 is a common root and 
 
 q q 
 
 T 
 so the formula must be of formal order at least k+1. Hence b £ = 
 
 where 
 
 T T log(P) 
 H fo P-I ' 
 
 and so 
 
 rp rp rp rp rp rp rp 
 
 (b +e|)S = b P = b + ej logP - b + e£ . 
 
 T T 
 Also £ = implies e S = e . Therefore the null space of S-I is at least 
 q ~q ~q 
 
 of dimension two, and so we can choose v ^ such that Sv = v and 
 
 T 
 e.v = 0. Second, suppose that £ ^ 0. Let 
 ~0~ q 
 
 v(K) = det(£l-S)-(a-S) _1 £ , 
 
 which is a vector of polynomials in £. Let £, be a common root of p(0 
 
 and a(£). Then v(£_) ^ 0, since e v(£) = (£-l) q £ . Also 
 
 (J ~ ~q~ q 
 
 (C I-S)y(£ Q ) = det(£ Q I-S).£ = , 
 
 T 
 and e n v(E. ) = a(£. ) = 0. Therefore in either case there exists E, £ 
 
 T T 
 
 and v 4- such that S = £_v and e~v = 0. Moreover, e,v = 0, because 
 ~ v 0~ ~0~ ~1~ 
 
 T 
 e.. is a left eigenvector for the eigenvalue 0. This means that if a_ 
 
 were changed to a A + v, the values of y and hy' would remain unchanged 
 ~U ~ n n 
 
 for all n, and yet a would become a + K„v. This proves the nonexistence 
 
 ~n ~n 0~ 
 
 of an expression for a in terms of the values y and y 1 .D 
 
 n n n 
 
 Theorems 5.3 and 5.4 do not cover the case where the order is 
 less than q due to a common factor of £-1, but extending the results to 
 such formulas does not seem worthwhile. 
 
33 
 
 6. Equivalence of Predictor-corrector Formulas. The use of a linear 
 
 Nordsieck formula requires the solution for <5 of the equation 
 
 n 
 
 hy' n + JL6 = hf(t , y n + 1.8 ) 
 ■'njO In n ^n,0 n 
 
 T T 
 
 where y _ = e^Pa , and hy 1 . = e.Pa , . In practice this must be 
 'n,0 ~0 ~n-l n,0 ~1 ~n-l 
 
 approximately solved by some iterative process: 
 
 6 n :- , 
 
 n,0 
 
 6 .. := 6 -[£l-£hJ ]~ 1 [hy* _+£_ 6 -hf ] , 
 n,y+l n,y 1 n,y n,0 1 n,y n,y 
 
 M = 0, 1,..., M(n)-1 , 
 
 where f ,1, and J , are defined as in §2 with y = y _ + &„6 
 
 n,y n,y-l n,y n,0 n,y 
 
 After determining 8 := 6 .,, N and adding the increment £6 a final 
 
 n n,M(n) ~ n 
 
 function evaluation may or may not be performed. For a P(EC)* Nordsieck 
 
 formula 
 
 a := Pa , + 18 
 ~n ~n-l ~ n 
 
 and for a P(EC)*E Nordsieck formula 
 
 a := Pa . + £6 + e.h(f -y ' ) , 
 ~n ~n-l ~ n ~1 n n 
 
 where hy' := hy' . + 1.8 „, , n . 
 n y n,0 1 n,M(n)-l 
 
 Remark. The predictor-corrector Nordsieck formulas of Gear 
 
 [1971a] are introduced independently of linear ("corrector only") 
 
 Nordsieck formulas. For this reason these predictor-corrector formulas 
 
 use l^l-l n hJ := 1. 
 
 1 n,y 
 
 T 
 For the P(EC)* formulas we have that hy' = e„Pa n + 16 , 
 
 n ~1 ~n-l 1 n 
 
 and so the recurrence can also be expressed as 
 
 a = Sa . + £hy' . 
 ~n ~n-l ~ n 
 
 Therefore the equivalence results apply to P(EC)* formulas as much as 
 
 they do to linear formulas. The underlying multistep predictor formula 
 
34 
 
 T 
 can be obtained from y . = e^Pa , by expressing a n as linear 
 
 n,0 ~0 ~n-l ~n-l 
 
 combinations of the components of y _. ; thus it has stepnumber of at most 
 
 k+0*. 
 
 The situation is quite different for P(EC)*E formulas. 
 
 To determine the equivalent multistep formula, one begins with the 
 
 recurrence 
 
 a = Sa , + (i-ejhy 1 + e,hf . 
 ~n ~n - i ~ ~1 n ~1 n 
 
 Proceeding as in the proof of Theorem 5.2, one obtains a difference 
 
 equation involving y . , j = 0(l)k, y 1 . , j = 0(l)k-l, and y' ., 
 
 j = l(l)k, which is not a true P(EC)*E multistep formula. 
 
 For example, consider the three-value Nordsieck formula. 
 
 By expressing y and y . as functions of y 1 , y' , , f ., , and a _ 
 n n-1 n n-1 n-1 ~n-2 
 
 and then eliminating h y" /2, one obtains the difference equation 
 
 n-z 
 
 \ + (2, 2- 2) Vl + (1 - U lK-2 
 
 - 4 h " y n + [(1 ^0 )hy ;-l + ^2- £ )h Cl ] + ( % +£ 2" 1)hy ;-2 > 
 which is not a P(EC)*E formula unless £ = £ . For the third order 
 Adams-Moulton formula this is 
 
 y„ - Vi ■ Ti K + l T2 K-i + n h K-i ] + s hy ;-2 • 
 
35 
 
 7. Applications. An important practical consequence of the equivalence 
 theory is that all multistep formulas are minimum storage formulas. 
 This idea is certainly implicit in the investigations by Gupta and 
 Wallace of new multistep formulas. It is also the rationale for the 
 computer research of Kong [1977] for k-step formulas of order k having the 
 smallest error coefficient for a given angle a of A(a)-stability. Previous 
 searches by Dill and Gear [1971] (the error coefficient for their formula 
 should be 14.0 rather than 0.0108) and by Jain and Srivastava [1970] 
 concentrated on formulas for which most of the trailing coefficients of 
 a(£) were set to zero. The computer results of Kong suggested formulas 
 which lead to a proof of the following result: for any a < tt/2 there 
 exists an A(a)-stable k-step formula of order k. This is also proved 
 by Grigorieff and Schroll [1977]. 
 
 The paper of Nordsieck mentions that "the potential advantage 
 of a more elaborate procedure in which the matrix hf is numerically 
 
 y 
 
 computed at every step and £ is made a chosen function of hf , implying 
 a nonlinear process tailored to the subject differential equation system, 
 is an interesting topic for future investigation," where we have 
 substituted our notation for his. The idea is developed further in a 
 report of Gear [1966], which proposes a formula for scalar implicit 
 differential equations of any order. For the equation y' - f(t, y) = 
 the formula becomes 
 
 I :- iVhf 
 
 y 
 
 and 
 
 £,l-Jl n hJ := l+h 2 f 2 
 1 n,u y 
 
36 
 
 where AM and BD refer to the (q+l)-st Adams-Moulton and q-th order 
 backward differentiation formulas, respectively. (The subscripts y 
 in the equation on page 21 of this report should read a . ) This idea 
 
 q 
 
 is extended to systems of differential equations by Skeel and Kong [1977], 
 who suggest 
 
 y 
 
 and 
 
 ItI-IAiJ : = (1-chf ) 2 
 
 1 n,y y 
 
 where y an< 3 c are free parameters. In §4 it was shown that p(£) and o(E,) 
 
 are linear transformations of Z, and so a linear combination of I 
 
 vectors corresponds to that same linear combination of the corresponding 
 
 linear multistep formulas. Thus, under the assumption of constant hf 
 
 y 
 
 the "blended" Nordsieck formula is equivalent to the same blend of linear 
 multistep formulas. 
 
37 
 
 BIBLIOGRAPHY 
 
 G. D. BYRNE & A. C. HINDMARSH (1975), "A polyalgorithm for the 
 numerical solution of ordinary differential equations," ACM Trans. Math. 
 Software, v. 1, pp. 71-96. 
 
 G. G. DAHLQUIST (1956), "Numerical integration of ordinary 
 differential equations," Math. Soandinavioa, v. 4, pp. 33-50. 
 
 G. G. DAHLQUIST (1963), "A special stability problem for linear 
 multistep methods," BIT, v. 3, pp. 27-43. 
 
 G. G. DAHLQUIST (1975), On Stability and Error Analysis for 
 Stiff Non-linear Problems } Part I, Report TRITA-NA-7508, Dept. of Computer 
 Sci. , Roy. Inst, of Technology, Stockholm. 
 
 J. DESCLOUX (1963), A Note on a Paper by A. Nordsieok, Report #131, 
 Dept. of Computer Sci., Univ. of Illinois at Urbana-Champaign. 
 
 C. DILL & C. W. GEAR (1971), "A graphical search for stiffly 
 stable methods for ordinary differential equations," J. ACM, v. 18, pp. 75-79. 
 
 C. W. GEAR (1966), The Numerical Integration of Ordinary 
 Differential Equations of Various Orders 3 Report //ANL-7126, Argonne National 
 Lab. , Argonne, Illinois. 
 
 C. W. GEAR (1971a), Numerical Initial Value Problems in Ordinary 
 Differential Equations } Prentice-Hall, Englewood Cliffs, N.J. 
 
 C. W. GEAR (1971b), "Algorithm 407: DIFSUB for solution of 
 ordinary differential equations," Comm. ACM, v. 14, pp. 185-190. 
 
 R. D. GRIGORIEFF & J. SCHROLL (1977), Uber A (a) -stabile Verfahren 
 hoher Konsistenzordnung , Nr.34, Fachbereich Mathematik (3), Technische 
 Universitat Berlin. 
 
 G. K. GUPTA (1976), "Some new high-order multistep formulae for 
 solving stiff equations," Math. Comp. , v. 30, pp. 417-432. 
 
 G. K. GUPTA & C. S. WALLACE (1975), "Some new multistep methods 
 for solving ordinary differential equations," Math. Comp., v. 29, pp. 489-500. 
 
 P. HENRICI (1962), Discrete Variable Methods in Ordinary 
 Differential Equations, Wiley, New York. 
 
 A. C. HINDMARSH (1974), GEAR: Ordinary Differential Equation 
 Solver, UCID-3001, Rev. 3, Lawrence Livermore Lab., Univ. of California, 
 Livermore, Calif. 
 
 M. K. JAIN & V. K. SRIVASTAVA (1970), High Order Stiffly Stable 
 Methods for Ordinary Differential Equations, Report #394, Dept. of 
 Computer Sci. , Univ. of Illinois at Urbana-Champaign. 
 
38 
 
 D. E. KNUTH (1968), The Art of Computer Programming, Vol. 1: 
 Fundamental Algorithms, Addison-Wesley, Reading, Mass. 
 
 A. NORDSIECK (1962), "On numerical integration of ordinary 
 differential equations," Math. Comp. , v. 16, pp. 22-49. 
 
 M. R. OSBORNE (1966), "On Nordsieck's method for the numerical 
 solution of ordinary differential equations," BIT, v. 6, pp. 51-57. 
 
 R. D. SKEEL (1973), Convergence of Multivalue Methods for 
 Solving Ordinary Differential Equations, Report TR73-16, Dept. of 
 Computing Sci. , Univ. of Alberta, Edmonton. 
 
 R. D. SKEEL & L. W. JACKSON (1977), "Consistency of Nordsieck 
 Methods," SIAMJ. Numer. Anal., v. 14, pp. 910-924. 
 
 R. D. SKEEL & A. K. KONG (1977), "Blended linear multistep 
 methods," ACM Trans. Math. Software, v. 3, pp. 326-345. 
 
 C. S. WALLACE & G. K. GUPTA (1973), "General linear multistep 
 methods to solve ordinary differential equations," Austral. Comput. J. , 
 v. 5, pp. 62-69. 
 
BIBLIOGRAPHIC DATA 
 SHEET 
 
 1. Report No. 
 
 UIUCDCS-R-78-940 
 
 3. Recipient's Accession No. 
 
 4. Title and Subt itle 
 
 EQUIVALENT FORMS OF MULTISTEP FORMULAS 
 
 5- Report Date 
 August 1978 
 
 7. Author(s) 
 
 Robert D. Skeel 
 
 8. Performing Organization Rept. 
 
 No - UIUCDCS-R-78-940 
 
 9. Performing Organization Name and Address 
 
 Department of Computer Science 
 222 Digital Computer Laboratory 
 University of Illinois at U-C 
 Urbana r Illinois 6 1801 USA 
 
 10. Project/Task/Work Unit No. 
 
 11. Contract /Grant No. 
 
 AFOSR-75-2854 
 
 12. Sponsoring Organization Name and Address 
 
 Air Force Office of Scientific Research/NM 
 Boiling AFB DC 20332 
 
 13. Type of Report & Period 
 Covered 
 
 Research 
 
 14. 
 
 15. Supplementary Notes 
 
 16. Abstracts 
 
 For uniform meshes it is shown that any linear k-step formula can be formula can be 
 formulated so that only k values need to be saved between steps. By saving an 
 additional m values it is possible to construct local polynomial approximations of 
 degree k+m-1, which can be used as predictor formulas. Different polynomial bases 
 lead to different equivalent forms of multistep formulas. In particular, local 
 monomial bases yield Nordsieck formulas. An explicit one-to-one correspondence is 
 established between Nordsieck formulas and k-step-f ormulas of order at least k, and 
 a strong equivalence result is proved for all but certain pathological cases. 
 Equivalence is also shown for P(EC)* formulas but not for P(EC)*E formulas. 
 
 17. Key Words and Document Analysis. 17a. Descriptors 
 
 linear multistep formula, multistep formula, linear multistep method, multistep 
 method, Nordsieck method, multivalue method, predictor-corrector method. 
 
 17b. Identifiers/Open-Ended Terms 
 
 17c. COSATI Field/Group 
 
 18. Availability Statement 
 
 Approved for public release; distribution unlimited 
 
 FORM NTIS-35 (10-70) 
 
 19. Security Class (This 
 Report) 
 
 UNCLASSIFIED 
 
 20. Security Class (This 
 
 Page 
 UNCLASSIFIED 
 
 21. No. of Pages 
 
 38 
 
 22. Price 
 
 USCOMM-DC 40329-P7I 
 
SECURITY CLASSIFICATION OF THIS PAGE (When Data Entered) 
 
 REPORT DOCUMENTATION PAGE 
 
 READ INSTRUCTIONS 
 BEFORE COMPLETING FORM 
 
 1. REPORT NUMBER 
 
 2. GOVT ACCESSION NO. 
 
 3. RECIPIENT'S CATALOG NUMBER 
 
 4. TITLE (and Subtitle) 
 
 EQUIVALENT FORMS OF MULTISTEP FORMULAS 
 
 5. TYPE OF REPORT & PERIOD COVERED 
 
 RESEARCH 
 
 6. PERFORMING O^G. REPORT NUMBER 
 
 UIUCDCS-R-78-940 
 
 7. AUTHOR(s) 
 
 ROBERT D. SKEEL 
 
 8. CONTRACT OR GRANT NUMBER(s) 
 
 AFOSR-75-2854 
 
 9. PERFORMING ORGANIZATION NAME AND ADDRESS 
 
 Department of Computer Science 
 222 Digital Computer Laboratory 
 Univ. of 111., Urbana, IL 61801 USA 
 
 10. PROGRAM ELEMENT, PROJECT, TASK 
 AREA 4 WORK UNIT NUMBERS 
 
 61102F 
 
 2304/A3 
 
 11. CONTROLLING OFFICE NAME AND ADDRESS 
 
 Air Force Office of Scientific Research/NM 
 Boiling AFB DC 20332 
 
 12. REPORT DATE 
 
 August 1978 
 
 13. NUMBER OF PAGES 
 38 
 
 14. MONITORING AGENCY NAME ft ADDRESSf// ditterent from Controlling Office) 
 
 15. SECURITY CLASS, (of this report) 
 
 UNCLASSIFIED 
 
 15a. DECLASSIFICATION/ DOWN GRADING 
 SCHEDULE 
 
 16. DISTRIBUTION ST ATEMEN T (of this Report) 
 
 Approved for public release; distribution unlimited, 
 
 17. DISTRIBUTION STATEMENT (of the abstract entered in Block 20, it different from Report) 
 
 18. SUPPLEMENTARY NOTES 
 
 19. KEY WORDS (Continue on reverse side if necessary and identify by block number) 
 
 linear multistep formula, multistep formula, linear multistep method, 
 multistep method, Nordsieck method, multivalue method, predictor- 
 corrector method 
 
 20. ABSTRACT (Continue on reverse side If necessary and identify by block number) 
 
 For uniform meshes it is shown that any linear k-step formula can be formulated 
 so that only k values need to be saved between steps. By saving an additional 
 m values it is possible to construct local polynomial approximations of 
 degree k+m-1, which can be used as predictor formulas. Different polynomial 
 bases lead to different equivalent forms of multistep formulas. In particular, 
 local monomial bases yield Nordsieck formulas. An explicit one-to-one 
 correspondence is established between Nordsieck formulas and k-step-f ormulas 
 
 DD , :° N RM 7 3 1473 
 
 EDITION OF 1 NOV 65 IS OBSOLETE 
 
 SECURITY CLASSIFICATION OF THIS PAGE (When Data Entered) 
 
ECURITY CLASSIFICATION OF THIS PAGEfWTifi Did Bnlurtd) 
 
 of order at least k, and a strong equivalence result is proved for all but 
 certain pathological cases. Equivalence is also shown for P(EC)* formulas 
 but not for P(EC)*E formulas. 
 
 , 1 4 1978 
 
 SECURITY CLASSIFICATION OF Tuic PAGE(TW)»n Dmta Entered) 
 
AUG itsn.