H 
 
 ^H 
 
 
LIBRARY OF THE 
 
 UNIVERSITY OF ILLINOIS 
 
 AT URBANA-CHAMPAIGN 
 
 510.84 
 
 IZ6t 
 
 no. 131 -140 
 cop. 3 
 
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 
 
 nov o i m 
 
 SFP 1 9 !<i 
 
 W 
 
 L161— O-1096 
 
Digitized by the Internet Archive 
 in 2013 
 
 http://archive.org/details/noteonpaperbyano131desc 
 
'5/064 
 Ll<br 
 
 no - ' 3 ' DIGITAL COMPUTER LABORATORY 
 
 " UNIVERSITY OF ILLINOIS 
 
 URBANA, ILLINOIS 
 
 REPORT NO, 131 
 NOTE ON A PAPER BY A. NORDSLECK 
 
 by 
 
 Jean Descloux 
 
 January 11^ 1963 
 
 This work was supported in part by the 
 National Science Foundation 
 under grant Gl6^-89 
 

 
 
 
 *)fO 
 
 
 J0.3 
 
 
 INTRODUCTION , 
 
 . . . . 
 
 1. THE ADAMS 
 
 METHODS 
 
 TABLE OF CONTENTS 
 
 Page 
 
 ...................... 1 
 
 ...................... 2 
 
 2. AN ASYMPTOTIC RELATION FOR THE ERROR OF INTEGRATION WITH THE 
 
 ADAMS METHODS . . . . ....... . . . . . . . . . . . . . k 
 
 3. SELF- STABILIZING PROPERTY OF THE ADAMS METHODS ......... l8 
 
 k. THE NORDSIECK METHOD ...................... 22 
 
 5. PREDICTOR FORMULAS ....................... 36 
 
 6. MODIFICATION OF THE LENGTH OF THE STEP ............. 37 
 
 7. STARTING PROCEDURES ...................... 42 
 
 8. CONCLUSION ........................... 56 
 
 9. TABLE OF RESULTS ........................ 57 
 
 APPENDIX 1 ............................ . 70 
 
 APPENDIX 2. SOME EQUIVALENT "METHODS" ............... 72 
 
 •1- 
 
INTRODUCTION 
 
 This report is an attempt to analyze Nordsieck's paper [l] with the 
 help of the asymptotic methods developed Toy Henrici in [2], 
 
 The definitions of the Adams -Bashforth formulas and the Adams -Moult on 
 formulas are contained in Section 1. 
 
 An asymptotic relation for the errors in the Adams methods is derived 
 in Section 2. Most of the arguments used are directly suggested by [2]. 
 
 Section 3 is concerned with the self-stabilizing property of the 
 Adams methods . 
 
 Section h shows how to get Nordsieck-like formulas from multistep 
 methods; applying this to the Adams -Moult on method,, one gets the Nordsieck 
 formulas; a modified Nordsieck formula is suggested „ 
 
 Predictor formulas for the Nordsieck method and the modified Nordsieck 
 formula are examined in Section 5 and interpreted in the language of multistep 
 methods . 
 
 The errors involved in the process of changing the step are analyzed 
 in Section 6. 
 
 Section 7 deals with starting procedures. 
 
 Section 8 gives some conclusions . 
 
 A table of results is presented in Section 9° 
 
 Appendix 1 is a complement of Section 2, 
 
 In Appendix 2, a formal generalization of the Nordsieck method is 
 discussed; it turns out that these formulas are "almost equivalent" to the 
 classical multistep methods. As a particular case, to every multistep method 
 corresponds an equivalent method using less "past values . " 
 
 No serious practical computations have been done for checking the 
 theoretical results of this report. However, some verifications of the 
 Nordsieck method and of the modified Nordsieck method, k = 2, with the dif- 
 ferential equation y 1 = y, have been satisfactory. 
 
 Note: In this report, all variables are real. 
 
 -1- 
 
SECTION 1. THE ADAMS METHODS 
 Let y(x) be a function and y'(x) be its derivative. 
 
 Adams -Bashforth Formulas (Predictor Type) 
 
 Let P(x) be a polynomial of degree k + 1 satisfying the relations 
 
 P(x Q ) = y(x Q ) 
 
 P'(x - ih) = y'(x - in), i = 0,1, ...,k; 
 
 th 
 P(x + h) is the result of the k Adams -Bashforth formula. It is easy to show 
 
 that there exist the constants (independent of h, y, x) (3 , (3 , „.., ft such 
 
 that 
 
 P(x Q + h) = y(x Q ) + h(p y'(x ) + P^'CXq - h) + ... e k y'(x Q - kh)); 
 
 if y(x) possesses k + 2 continuous derivatives, the difference 
 L(y(x ), h) - y(x + h) - P(x n + h) satisfies the relation 
 
 L(y(x Q ), h) = D k h k+2 y (k+2) (| ), x Q - kh < | < x Q + h, 
 
 where D, is a constant . 
 k 
 
 Adams -Moult on Formula (Corrector Type) 
 
 Let Q(x) be the polynomial of degree k + 2 satisfying the relations 
 
 Q(x Q ) = y(x Q ) 
 
 Q*(x Q - ih) = y'(x - ih), i = -1,0, ...,k; 
 
 Q(x + h) is the result of the k Adams -Moulton formula and satisfies a relation 
 of the form 
 
 -2- 
 
Q(x Q + h) = y(x Q ) + hCp^y'CxQ + h) + P Q y' (x Q - L) ... + ^y' (x Q - kh)) 
 
 where |3 , (3_. ..., (3 are constant. If y(x) possesses k + 3 continuous 
 derivatives, the difference L(y(x )_, h) = y(x + h) - Q(x + h) is given by 
 
 L(y(x Q ), h) = C k h k+3 (0, x Q - kh < | x Q + h, 
 
 where C, is a constant 
 k 
 
 -3- 
 
SECTION 2. 
 AN ASYMPTOTIC RELATION FOR THE ERROR OF INTEGRATION WITH THE ADAMS METHODS 
 
 Notations 
 
 Any letter with an arrow is a vector of dimension m; the components 
 
 m 
 of any vector y are denoted by y , y , ..., y . |y| = Z |y. | is the norm 
 
 -p "* i=1 
 
 of y. 
 
 We shall consider the integration of a system of m differential 
 equations on the interval b5 x S c > however, in order to allow the possibility 
 of initial values for x < b, we shall give the assumptions on the differential 
 equations for the interval a < x < c, where a < b. 
 
 h denotes the step of integration and we write x = b + nh, where n 
 
 n ' 
 
 is a positive or negative integer. For any function y(x) defined on the complete 
 
 interval [a,c] or only at the discrete points x , we use the notation y(x ) = y . 
 ' * n' n n 
 
 If t is a variable, 0(t) represents a quantity such that 
 
 |0(t)| < w|t|. where w is independent of t; t may be restricted to some 
 
 interval. 0(t) is a vector such that | |o(t)| - 0(t). We shall use in the 
 
 following notations of this type: 
 
 y(x) - z(x) + 0(h P ); 
 
 we always suppose that 0(h) is uniform in x for a < x < c, i.e., there exists 
 a constant w, independent of x and h such that |0(br )| < w|h | for a < x < c. 
 
 For any function y(x), y (x) denotes its i derivative. 
 
 Theorem 
 
 Let 
 
 |j 7<x) - ?(x,y) (1) 
 
 be a system of m differential equations and 
 
 4- 
 
g ±i (x,y) = — , i,J = 1,2,..., m 
 
 be the elements of the matrix G(x,y). 
 
 — » — > 
 We assume that for any y and y*, i,j = 1,2, . . „,m, a < x < c: 
 
 a) g. . (x,y) exists and is continuous in x and y; 
 
 b) |g i . (x,y)| < - , . n constant f2) 
 
 c) I g, .(*,?*) - g .(x,y)| <if\\y* - y||, if constant (3) 
 
 Let z(x) be a particular solution of (l); we assume that for 
 a < x < c. 
 
 d) z(x) possesses (p + 2 ) continuous derivatives (see (5)); 
 
 Let 
 
 e) g..(x,z(x)) has a continuous derivative in x, i,j = l ; ... ; m 
 
 L(y(x),h) = y(x + h) - y(x) - h(j3_ iy '(x + h) + ... + ^y' (x - kh)) (k) 
 
 be any formula described in Section 1 for k > 1; if p = 0, it is an Adams - 
 Bashforth formula; if 3 / 0, it is an Adams -Moulton formula; let C and p 
 (p > 2) be such that 
 
 L(y(x),h) = Ch P+1 y (P+l) (e ), x - kh < | < x + h (5) 
 
 Finally let u satisfy the difference equation 
 
 u ,. = u + h(p ,f ,, + p_f + ... a? . ) (6) 
 
 n+l n v -1 n+1 n k n-k v 
 
 -> ->. -»• . 
 where f = f (x , u ) for n > q: q is a given integer such that x , > a; the 
 n x n' n — U7 ^ ° q-k — 
 
 initial values u , f ,,..., f , satisfy the relations 
 
 q' q-1 7 ' q-k 
 
 f) u =z + h 6 + 0(h ), 5 constant, s integer > 1 (7) 
 
 f . = z 1 . + h T). + 0(h ). i = 1,2, ...,k. r\ . constant, r integer > 1, 
 q-iq-i l " ' ' ' ' i ' — 
 
 (8) 
 
 -5- 
 
Remark that u , . u _, .... u . are not necessarily defined. (6) is 
 q-1/ q-2' ' q-k J v y 
 
 an implicit equation for u ; in order to provide the existence of a unique 
 
 solution, we assume (see Appendix l) that h < h, where 
 
 g) s ^ssnr^i ; rsJi ♦ ... ie k n (9) 
 
 Then, under hypotheses a to g, the error d(x) = u(x) - z(x) is given by 
 
 d(x) = ^(e^x) + 0(h)) + h S (e n (x) + 0(h)) + h r+1 (e m (x) + 0(h)) (10) 
 
 where e , e , e are the solutions of the systems of differential equations 
 
 ^ e z (x) = G(x,z(x))e J (x) - Cz (p+l) (x), e^b) = 0; (ll) 
 
 ^ e n (x) = G(x,z(x))e i;[ (x), e (b) - 8; (12) 
 
 ~ e XII (x) = G(x,z(x))e Ilx (x), ^III (b) = ^ ^ 13 ^ 
 
 n - (^ + p 2 + . . . + e k )\ + O s + P 3 + . . . + P k H 2 + • • • + P k \- 
 
 Generalized Theorem 
 
 In the preceding theorem, we have assumed that z(x) is a fixed solution 
 of (l). In fact, the same conclusions hold if z(x) depend on h, provided we 
 add the two following assumptions: 
 
 h) there exists a constant U, independent of x and h such that for 
 a < x < c, < h < h 
 
 II? (P+2) (X)|| <U; 
 i) there exists a constant V such that 
 
 5 (x,z(x))| < V, i,j = 1, ...,m. 
 -'-J 
 
 •6- 
 
Remarks 
 
 1) The assumptions b and c are very strong since we ask them to be 
 
 -> — > 
 valid for all y and y*; however, it is possible a posteriori to 
 
 weaken them, since in the proof of the theorem we shall use only 
 
 a bounded region of y. 
 
 2 ) The proof of the generalized theorem is almost word for word the 
 same as the proof of the theorem itself and will consequently 
 not be given here. 
 
 Proof of the Theorem 
 
 All the symbols introduced so far in Section 2 keep their meaning 
 in the following , 
 
 Lemma 1 , 
 
 -* -> 
 Under assumptions b and c, for a < x < c and any y and y*, we have: 
 
 a) I |y* + y|| < ||y*ll + I l?l I J 
 
 b) | |f (x,y*) - f(x,y>| | < a\ |y* - y| |; 
 
 c) f(x,y*) - f(x,y) = G(x,y*)(y* - y) + 0(||y* - y|| 2 ) 
 
 d) I |G(x,y)y*| | < ft| \j*\ \ 
 
 Proof 
 
 Lemma la is a direct consequence of the properties of the absolute 
 values . For each component f.(x,y) of f(x,y), we have the mean-value theorem 
 
 m 
 f.(x,y*) - f. (x,y) = Z g,.(x,y**)(y* - y ) for some y**; (ik) 
 
 n m n 
 
 f.(x,y*) - f.(x,y)| <- Z |y* - y.| = - | |y* - y||; 
 
 j — -J- 
 
 ||f(x,y*) - f(x,y)|| < a||y* - y||, 
 
 which proves lemma lb„ 
 
 -7- 
 
In order to prove lemma lc, we start with equations (l^-); by (3): 
 
 g i(j (x,y**) = g i;j (x,y) + o(||7** - 711); 
 
 :ince y* - y. = 0( | |y* - y| | ) and | | y** -y||<j|y*-y||, it follows that 
 J J 
 
 f (x,7*) - f (x,y) = E g,.(x,y)(y* - y.) + 0(||y* - 711), 
 
 which is the desired result for the i component. 
 
 Lemma Id can be proved in the same way as lemma lb. 
 
 Lemma 2 
 
 Let y_ k , y_ k+1 .> •••, y Q , y 1} J 2 , ••• satisfy the inequalities 
 
 y . < y + h(b n y . + b_y + . . . + b. y , ) + 7 for n > 0: (is) 
 J n+1 - J n v -l^n+l J n k^n-k' - * v ' ' 
 
 < y. t < Y ; i = 1,2, ...,k; y Q > 0; b. > 0, i = -1,2, ...,k; 7 > 0; 
 
 let B = b , + b^ + ... + b. and suppose that 
 -10 k 
 
 h > 0, B > 0, h<|g; (16) 
 
 then 
 
 y n < (^(y + hBY + 7) + ^)e 3nhB for n > 0. (17) 
 
 Proof 
 
 From (l6), it follows that 
 
 1 + hB < 1.5; 
 
 -8- 
 
Thx ' 1 + r^hB 5 1 - 2hB < (1 + hB) 2 ; (18) 
 
 from (l6) and (l8), we get 
 
 o 
 
 y , < (l + hB) fy + h(b^y + b_,y _ + . . . + b. y , ) + 7} ; 
 J n+1 - v lJ n v Crn lrn-1 k J n-k y J ' 
 
 since y . < Y for i = 1,2, . . . ,k, we can write for n = 0,1,.. . ,k-l: 
 
 y _ < (1 + hB) {y + (b_y + b n y + ... + b y ) + hBY + 7}? 
 
 •'n+l — n v On l^n-1 n •" 
 
 from the two last inequalities, it follows by induction that 
 
 y < v for n > 0, where v., = y~, and 
 J n - n - ' J 0' 
 
 v . = (1 + hB) O + h(b_v + b_v _ + ... b v.) + hBY + 7} for < n < k-1, 
 
 n+1 n n 1 n-1 n — - 
 
 v . = (l + hB) fv + h(b_v + b n v . + . . . b v . ) + 7} for n > k 
 
 n+1 L n n 1 n-1 k n-k ' - 
 
 it follows that v < v < v ...; consequently 
 
 v , < (l + hB) 2 f(l + hB)v + hBY + 7} for < n < k - 1, 
 
 n+1 — n — — ' 
 
 < (1 + hB) 2 {(l + hB)v + 7} for n > kj 
 
 n+1 - L v n 
 
 2 ^ 
 
 by induction again, since (l + hB) < 3> (l + hB) < ^-, we have 
 
 v < w for n > 0, where w. = y^ and 
 n - n - y J 
 
 w = 4w + 3(hBY + 7) for < n < k - 1; (19) 
 
 w , = (1 + hB) 3 w + 37 for n > k; (20) 
 
 n+1 v / n - ? v 
 
 -9- 
 
(19) is a difference equation the solution of which is given (as it is easy 
 to prove by induction) by 
 
 4 n 1 
 w = 4 w„ + 3 i " n (hBY +7) for < n < k: 
 n 4-1 — — 7 
 
 in particular, for n = k 
 
 w = 4 k w + (4 k - l)(hBY + 7) < ^ k (w + hBY + 7); (2l) 
 
 in the same way, the solution of (20 ) is 
 
 ,. ^3(n-k) (1 + hB) 3(n " k) -l 
 = (1 + hB) v w + - ■ - 
 
 w - (1 + hBj^ v ' w, + ■* '- 37 for n > k 
 
 n 
 
 (1 + hB) 3 ■ 1 
 clearly, we can write for n > 0: 
 
 f 
 
 / f-, v,^^3n 0-v (1 + hB) 3n - 1 
 < (1 + hB) w. +37 - -z i 
 
 n - v ' k 
 
 (1 + hB) 3 - 1 
 
 from the inequalities 
 
 (1 + hB) 3n < e 3nhB , n > 0, 
 
 < 1 
 
 (1 + hB ) 3 - 1 " 3hB 
 and by (21 ), we get 
 
 w ■: (4 k (w n + hBY + 7) + r^)e 3nhB , n > 
 n — v v hB ' — 
 
 which is the desired result 
 
 -10- 
 
Lemma 3 
 
 If a function y(x) possesses a second continuous derivative on 
 
 a < x < c, then 
 
 L(y(x),h) = 0(h 2 ) 
 
 Proof 
 
 First consider the case y(x) = x„ By (5) it follows that L(x,h) = 0, 
 i.e., 
 
 P_ x + P + ... + 3 k = i; (22) 
 
 now, let y(x) be any function satisfying the hypothesis; since y'(x) is absolutely 
 continuous on a < x < c, y' (x + 9 ) = y' (x ) + 0(h) for any a < x - kh < < x + h < c 
 
 L(y(x),h) = y(x) + h(y'(x) + 0(h)) - y(x) - hy'txKp^ + P Q + ... + \ + 0(L)); 
 
 by (22), L(y(x),h) = 0(h 2 ), q.e.d. 
 
 Lemma h 
 
 If e(x) is the solution of the systems (ll) or (12), under assumptions 
 d and e, then e(x) possesses a continuous second derivative. 
 
 Proof 
 
 By lemma Id, the system of differential equations have a Lipschitz 
 condition which insures, under assumptions d and e a first continuous derivative 
 for e(x) (see for example [2]). The lemma h results by differentiating the 
 equations (ll) or (12). 
 
 Lemma 5 
 
 Let t(x,y) be a vector satisfying the Lipschitz condition 
 
 I I t(x,y*) - t(x,y) I < Il| |y* - y| I, a < x < c, II constant. 
 
 Suppose that the function y(x) satisfies the system 
 
 ■11- 
 
— y(x) - t(x,y(x)), a < x < c, 
 
 and that v satisfies the difference equation 
 n 
 
 v n = v + h(B ,t . + p.t + ... + at . ) + 7 . (23) 
 
 n+1 n v -1 n+1 n k n-k n' v J ' 
 
 — > — *. — * > , . 
 
 where t = t(x ,v j for n > j; j is a given fixed integer (positive or negative) 
 
 such that x , > a: let w = v - y and L(y(x),h) he the vector of components 
 j-k — ' n n n ww ' * 
 
 L(y,(x),h). If 
 
 1 
 
 1) I |w I j < W, n = j-1, j-2, . . .jj-kj W constant 
 
 2) I 7 I < T, n>j.r constant 
 1 ' n — — w * 
 
 3) I |L(y(x),h)| I A, a < x < c, A constant 
 
 then 
 
 ||w(x)|| < (+ k (||w. II + hBw + r +a) + L_±A) e 3( c - a ) B for x . < x < c 
 - j nB j - - 
 
 where B = Il( 1(3^1 + |p Q | + ... IpJ) and h < |^. 
 
 Proof 
 
 If we subtract the relation 
 
 y , = y + h(p n y' , + B^y' + ... + B, y' , ) + L(y ,h) 
 
 J n+1 J n v -l J n+l CTn K k J n-k y w n' y 
 
 from (23), we get 
 
 w _ = w + h(p _ (t , - y' ) + . . . + 6, (t . - y' , ) + 7 - L(y ,h); 
 n+1 n v -l v n+1 J n+1 k v n-k •'n-k n w n' f 
 
 Since t - y' < n| w |, we have by lemma la: 
 1 ' n J n - ' ' n' " 
 
 I |w . 1 1 < I |w 1 1 + h(n|B J I |w . 1 1 + ... + n|p. I I |w v | I ) + r + A; 
 
 ' n+1 — n 1 ' ' — 1 ' ' ' n+1 1 ' ■ k' ' n-k' ' ' 
 
 -12- 
 
by lemma 2, if we set y^ = | \^ n+ . | |, 7 = T + A, Y = V, B = Jld^J + ... \\\), 
 we get for n > 
 
 I -* ll^/),k/||-*|i 1T1T7 _ A \ T + Av 3nhB 
 w . < (4 ( w. + hBW + r + A ) + — -=— )e : 
 
 n+j ' ' - ' ' j hB ' 
 
 since a<x.<x .<c, it follows that nh < c - a, q.e.d. 
 J - n+j - - 
 
 Proof of the Theorem 
 
 We first consider a particular case. Suppose that for i = 1,2,,,, „.<,k 
 
 . -> 
 1} u . are defined; 
 q- 1 
 
 2) f . = f (x .oU . ); 
 q-i q-i q-i 
 
 3) u . = z . + 0(h ): (see 7). 
 ' q-i q-i " K 
 
 From the third condition, because of lemma lb, we have 
 
 f . = z' . + 0(h b ). 
 q-i q-i 
 
 Since by lemma lb, f(x,y; satisfies a Lipschitz condition with constant fl, we 
 
 — > — > — >. . — >. . — » -> — > —> 
 can apply lemma 5^ replacing t by f , II by Q, } y{x) by z(xj, v by u , w by d , 
 
 J by q, r by 0, W and | |w. | | by 0(h S ), A by 0(h P+ ) (see (5)); we get 
 
 J 
 
 ||d(x)|| •: (U k (0(h s ) + hB0(h s ) + 0(h p+1 )) + 0(h P ^l )e 3(c-a)B = Q(h s ) + o(h p ) = 
 
 — LID 
 
 (2k) 
 Now, let us subtract from (6) the corresponding equation for z(x): 
 
 d _ = d + h(p . (f . - z 1 . ) + p.(f - z' )+...+ ft (f . - z' )) - L(z ,h); 
 n+1 n -l v n+1 n+1 n n k n-k n-k v n' " 
 
 by lemma lc, for n > q - k: 
 
 f - z 1 = f(x,u ) - f(x,z ) = G(x ,z )d + 0( d ). 
 n n ' n ' n n'nn Ul n" 
 
 ■13- 
 
We introduce the notations 
 
 G(x) = G(x,z(x)); G n = G(x n ). 
 
 By (2^) we can write for n > q 
 
 d . - d + h(p ,G n d _+p.G d + ... +B G .d , ) - f(z ,h) + 0(h 2p+1 ) + 0{h 2s+1 ); 
 n+1 n -1 n+1 n+1 n n k n-k n-k v rr v ' v '> 
 
 on another side by assumption d, (5) and (22), we have 
 
 L(z ,h) - Ch P+1 (B n z ( P; l} + ... + p z (p + l} ) + 0(h P+2 ); 
 x n' -1 n+1 n n-k " 
 
 then 
 
 d , - t + h(p .(g n d - ch^P;^ + ... a (G v d v - ch p z (p + l} ) 
 
 n+1 n v -l v n+1 n+1 n+1 K k^ n-k n-k n-k 
 
 P+ 2 \ , 7t^2p+lx , ?t/v,2s+l- 
 
 + O(h^) + 0(h^ x ) + Oth"^). 
 
 We have to distinguish three cases : 
 
 -> -s— * ~~* 
 First case: p > s. Let d* = h d : d* satisfies the difference equation 
 n n 7 n 
 
 d* . = d* + h(p n G n d*+ ... +p,G .d* . ) + 0(h 2 ), 
 n+1 n -1 n+1 n k n-k n-k " 
 
 with the initial conditions 
 
 d* = 6 + 0(h), d*_ . = 0(1), i = 1,2, . . ,,k; 
 
 we compare d*(x) with the solution e(x) of the system of differential equations 
 
 f^ e(x) = G(x)e(x), e(b) = 5; 
 
 Let w = d* - e : we have 
 n n n' 
 
 -ik- 
 
w = 0(h), w . = 0(1), i = 1,2,... ,k; 
 
 by lemmas 3 and k, L(e(x),h) = 0(h ); we can apply lemma 5 with t(x,y) = G(x)y, 
 n = fi (lemma Id), v n = d*, y n = e n , W = 0(l), r = 0(h 2 ), A = 0(h 2 ), j = q, 
 w = 0(h): one gets 
 
 | |w(x)| | < (^(0(h) + hBO(l) + 0(h 2 ) + 0(h 2 ) + 0(h2) h ;° (h2) -)e 3(c - a)B = 0(h), 
 
 i.e., d*(x) = e(x) + 0(h) and d(x) = h S (e(x) + 0(h)), x < x < c. 
 
 ->■ -p-» -> 
 Second case: p = s. Let d* = h d ,° d* satisfies the difference equation 
 n n^ n 
 
 d* = d* + h(p n (G > - Cz (p ! l} ) + ... + MG J* . - Cz (p : l} )) + 0(h 2 ); 
 n+1 n v -l v n+1 n n+1 K k v n-k n-k n-k \ '* 
 
 using the same arguments as in the first case, we get 
 
 d(x) = h P (e(x) + 0(h)), x < x < c; 
 
 where e(x) satisfies the differential equation 
 
 ^ e(x) = G(x)e(x) - Cz (p+l) (x); e(b) - 5. 
 
 Third case : p < s. This case is identical to the second one, except that we 
 have to replace 5 "by 0. Clearly, the following relation is valid for the three 
 cases 
 
 d(x) = h P (e x (x) + 0(h)) + h S (e n (x) + 0(h)), 
 
 where e (x) and e (x) satisfy respectively (ll) and (12); this finishes the 
 proof of the theorem for the particular case. 
 
 Now, it is easy to prove the theorem for the general case. Subtracting 
 from (6) the corresponding equation for z(x) for n = q we get 
 
 ■15- 
 
Vl - d q + h ( p -l< f q+ l " Z Q + l' + P 0< f q " * q > + ••• + M f q-k " Z q-k» " L \>^> 
 
 (25) 
 
 by lemma lb, |f - z ! | | < fi| |d | | ; by lemma la and equations (5), (l), 
 
 (8), we have 
 
 ||d 1 1 < o(h s ) + h(|p Jn||d , || + 0(h S ) + 0(h 2 )) + 0(h P+1 ): 
 
 by (9), 1 - h|p . |ft > .5; consequently 
 
 l|d q+1 M = 0(h S ) + 0(h r+1 ) + 0(h P+1 ); 
 using this relation in (25) we get by (7) and (8): 
 
 Vi = h "* + *f*\ + *4**% + ••• + V r+1 \ + 3(hS+1) + 3(h " 2) + ° (hP+1) - 
 
 — * 
 
 In the same way, we can obtain for d : 
 
 v = v + e 2 hr+1 ^i + h^X + ••• + v r+1 \-i + 3(hS+1) + °< hr+2 > + °( hP+1 » 
 
 = h s f + h r+1 (p 1+ e 2 )n 1+ ...h r+1 (p k _ 1+ p k « k _ 1+ h I ' +1 p k r lk+ o(h 3+1 ) + o(h r+2 ) + o(h p+1 ); 
 
 repeating the same process for d } ..., d ; we finally get 
 
 d = h 6 + h + p + ... + p^ + h o + ... + p k )n 2 + ... + n p k i } 
 
 7*/,s+l N -^/. r+2 v -*/.p+ls 
 
 + 0(h ) + 0(h ) + 0(h^ ) 
 
 , s^ tT+1-* 7?/i_ s+ l\ 7*/, r+l\ t?/vP + 1n 
 
 = h 6 + h ^ + 0(h ) + 0(h ) + 0(h^ ), 
 
 -16- 
 
d q+k _. = 0(h S ) + 0(h r+1 ) + 0(h P+1 ), i = 1,2,... ,k, 
 
 We are now in position to apply the theorem proved above for the particular 
 case. It suffices to replace q by q+k and to remark that the system of 
 differential equations 
 
 ^e IX (x) = G(x)e n (x) 
 
 is linear, 
 
 ■17- 
 
SECTION 3- SELF- STABILIZING PROPERTY OF THE ADAMS METHODS 
 
 We use the same notations and conventions as in Section 2. Let 
 y = f (x,y) be a system of m differential equations satisfying the Lipschitz 
 condition 
 
 | |f(x,y*) - f(x,y)| | < a\ |y* - y| | f or a < x < c. (l) 
 
 ->. . -> 
 
 Let z(x) be a particular solution of the system and u the approximation 
 
 obtained by using the k Adams-Moulton formula (see Section l): 
 
 u . = u + h(p ,? , + . .. + a f , ). (2) 
 
 n+1 n -1 n+1 k n-k ' v 
 
 where x = b + nh, f » ?(x„,u^) for n > 0, < h < nCj( \ a \ — — | Q ' r\ 1 H. 
 
 :.u ;iorn^u, u^-n^ ._, / i • 1 1 _ ' i > 
 
 We suppose that the initial conditions for u in b are such that 
 
 u(x) = z(x) + 0(h), b < x < c; (2a) 
 
 let v(x) be the exact solution of the differential equations such that 
 
 v(c) = u(c); v(x) depends on h; we suppose that v(x) possesses k+3 continuous 
 
 derivatives and that there exists a constant W such that 
 
 | |v^ k+3 ^(x)| | < W for b < x < c, h < h. 
 
 Theorem 
 
 For any fixed integer r independent of h 
 
 u(c - ih) - v(c - ih) = 0(h k+3 ), i - 0,1, . ..,r (3) 
 
 Remark 
 
 The self-stabilizing property allows us to apply the so-called 
 "Milne's method for evaluating the local discretization error" even if the 
 initial conditions are "bad." 
 
 ■18- 
 
Let x = c (t depends on L) and u* be the value of the predictor 
 
 \th 
 obtained by the (k+l) Adams -Bashforth formula 
 
 H " Vi + h(p 5 ? t-i + p I ? t- 2 + • • • + ^.m'* (*) 
 
 in Section 1, we have mentioned that 
 
 v t = vi + «tfu ♦ - ♦ pfca-i. M ) ♦ w k+37(k+3) <*> (5) 
 
 Subtracting (5) from (k), we get by (l) and (3): 
 
 in the same way, comparing (2) with the analogous equation for v(x): 
 
 vvvrvr c / 3 ' (kt3) (^ 5 ^' <?> 
 
 where x . , , < £ , f] < x , . 
 
 t-k-1 — ~ ' — t 
 
 Subtracting (7) from (6): 
 
 u* - v t - h (k+3) (C k v( k+ 3)(,) - D k+1 v( k+ 3)( 5) ) + 0(h k+U ) ; (8) 
 
 if we suppose that v (x) is continuous in x and bounded on b < x < c, 
 < h < h, then 
 
 c k 7 (k+3) (^) - \ + / k+3) (i) - (c k - D k+1 )7( k+ 3) ( ,) + 3 {h); 
 
 replacing in (8), one gets an approximative value for v ( T l) which provides 
 an asymptotical evaluation of the truncation error C h v (n) in (7)« 
 
 ■19- 
 
Proof of the Theorem 
 Lemma 
 
 Under hypotheses (l) and (2a), 
 
 u(x) - v(x) = 0(h), b < x < c. 
 
 Proof 
 
 By (2a), it suffices to show that z(x) - v(x) = 0(h). By (l), we can 
 write 
 
 x 
 
 z(x) - v(x) = z(c) - v(c) + / (f(z(u),u) - f(v(u),u))du; 
 
 J c 
 
 x 
 
 I |z(x) - v(x)| I < ||z(c) - v(c)|| + / I |f(z(u),u) - f(v(u),u)||du 
 
 J c 
 
 c 
 
 < I |z(c) - v(c) I I + n / |z(u) - v(u) I |du; 
 
 x 
 
 from this last inequality, it is known that 
 
 ||z(x) - v(x)|| < ||z(c) - v(c)||e fi(c " x) . 
 
 Since z(c) - v(c) = z(c) - u(c) = 0(h), the conclusion follows. 
 
 Now we go to the proof of the theorem itself. Subtracting the two 
 equations 
 
 VVi th(p .iV ••• * Vt-i-i'' (9) 
 
 \ " Vi + h ^-A + • • • + \%,*i) + c k h k+3 7 (k+3) (i ), (10) 
 
 where x = c, x < f < x , we get by the lemma and (l) 
 
 -20- 
 
= u t _ x - ? + 0(h 2 ) + 0(h k+3 ) 
 
 i.e. 
 
 * i^ 1 i Af V> —' 
 
 V X - Vi = ° (h > + ° ( * } " (11) 
 
 Changing t in t-1 in (9) and (10), one gets because of (ll) 
 
 u t _ 2 - v ' = 0(h 2 ) + 0(h k+3 ); 
 
 "by induction, it follows that u - v = 0(h ) + 0(h ) for i = 0,1, ...,s, 
 
 Li *~ X — X 
 
 where s is any fixed number independent of h; with this result in mind, we 
 subtract again (10) from (9) and get this time: 
 
 3> , T*fJs*3 
 
 - u - v + O(h^) + O(h^); 
 
 by induction it follows that 
 
 k+3 
 
 u, 
 
 . - v = 0(h J ) + 0(h °) for i - 0, 1,2, . . ., s-k-1 
 
 X Xi ~* X 
 
 Again by induction for j = 2 } J,,k } . . . 
 
 = 0(h J ) + 0(h k+3 ) for i - 0,1,2,..., s-('j-2)(k+l); 
 
 u, . - v, . 
 t-i t-i 
 
 for J = k+3, we get (3) 
 
 ■21- 
 
SECTION k. THE NORDSIECK METHOD 
 
 Note 
 
 In order to simplify the notations and avoid confusion with different 
 kinds of vectors, we shall consider in this section and in the following ones 
 only scalar differential equations; generalizations for systems are always easy 
 and in most of the cases obvious. 
 
 Let y 1 = f(x,y) be a differential equation to be integrated for x > b 
 with the help of the general multistep method 
 
 W ■ Vn + °lVl + ••• + a /n-j + h ( p -l f n + l + Vn + •■• + P k f n-k>> n 2 ° 
 
 (1) 
 
 where 0c. } $. } k, j, (j > 0) are the constants which characterize the formula; 
 
 some of the a.'s or B.'s may vanish; furthermore 
 1 i ' 
 
 x = b + nh; f = f(x ,y ) for n > 1; 
 
 n ' n n' n — ' 
 
 (2) 
 
 Y , y ± , >-, y_y f , t. ml f • -., f_ k are the initial values; f Q , f_ 1 , 
 
 may 
 
 be different from f(x Q ,y ), f(x_ 1 ,y_ 1 ) 
 
 We operate some formal transformation on (l). Let u be the column 
 
 -*■ 
 vector of elements y,y ,, .,., y ., f.f n . .... f , and p the constant 
 
 n' •'n-l n-j' n' n-1 ' n-k 
 
 vector of elements a_. ol , .... a., hB_, hB, , .... hB, . (l) can be written in 
 m 0' 1 J 1 k 
 
 r-*T 
 
 the form (p being the transposed of p): 
 
 y . = p u + hB _f . 
 J n+1 ^ n -1 n+1 
 
 u . = Au + f , a 
 n+1 n n+1 
 
 f n+l = f(x n+r y n+l ) 
 
 y for n > 
 
 (3) 
 
 M 
 
 (5) 
 
 where 
 
 -22- 
 
'% a i •'• a j| hp o hp i • 
 
 •■ h V 
 
 10 1 
 
 
 
 
 
 10 
 
 
 
 
 
 
 1 1 
 
 
 | * 
 
 
 1 
 
 1 o . 
 
 rvi 
 
 a = 
 
 
 1 
 
 
 the matrix has been subdivided into four submatrices; the upper-left submatrix 
 is square of order j+1; the lower-right submatrix is square of order k+l„ 
 
 For given u , (3)> (^); (5) define an iterative process; for n - 0,, 
 (3) and (5) allow to compute f and (k) u , etc.... 
 
 Let T be a regular matrix of order k+j+2, the elements of which depend 
 only on h; we consider the linear transformation 
 
 u = Tv : 
 n n' 
 
 replacing in (3), (M, (5)j we get for v the iterative process 
 
 -»T- 
 
 y . = q v + hp _f . 
 n+1 n -1 n+1 
 
 v , = Bv + f _,b > for n > 
 n+1 n n+1 ' 
 
 f n+l = f (V y n + l ] 
 
 (6) 
 (7) 
 (8) 
 
 ■* _T-* -1. -+ -l-> 
 where q = T p, B = T AT, b = T a. 
 
 -23- 
 
The process (6), (7), (8) is equivalent to the process (3); (M; (5) 
 — > — » — > — > 
 
 in the following sense: if u and v are chosen such that u = Tv , then both 
 
 processes provide the same values for y and furthermore u = Tv for all n > 0, 
 
 r n n n - 
 
 In order to get a particular matrix T, let us consider the polynomial 
 P (x) of degree k+j+1 satisfying the relations 
 
 n 
 
 P (x . ) - y . for i = 0,1, . ..,j; 
 n n-i' n-i ' 7 
 
 P'(x . ) = f . for i = 0,1, . . .,k; 
 n n-i n-i ' ' 7 7 . 
 
 (9) 
 
 this polynomial exists and is uniquely defined; in order to prove that, it 
 suffices to show that the only polynomial Q (x) of degree k+j+1 satisfying the 
 conditions 
 
 S/Vi^ = °' i = °* 1 '-"'J 
 
 Q*(x . ) = 0, i = 0,1, ...,k 
 tl n-i ' 7 ' ' 
 
 is the null polynomial; by Rolle ' s theorem, there exist | , | , ..., |. such 
 that 
 
 x > |_ > x . > 6_ ... > 6. > x ., 
 
 n 1 n-1 2 j n-j' 
 
 Q^(ft ± ) = 0, i = i,...,j; 
 
 QlCx) is of degree k+j and possesses at least k+j+1 zeros; consequently 
 Q^(x) s and, since Q^xJ = 0, Q^x) = 0, q.e.d. 
 
 Following Nordsieck's notations, let us write P (x ) in the form 
 
 ' n 
 
 p ( * U) h ( (2)^ X " ^ (3)/ X " X nV (k+j+2)/ X ' X n^ + 
 
 P (x) = v x + h(v v I — + v — ) + ... + v 
 
 nn K n\hjn\hy n V h J 
 
 (10) 
 
 "• (l) fk+i+2) ~» 
 
 let v denote the vector of components v ..... v and recall that u 
 n n ' ' n n 
 
 has the components y, ..., y . , f , .... f ,; then the system (9) can be written 
 
 n' ' n-j' n ; ; n-k' 
 
 in the form 
 
 -24- 
 
Tv = u : 
 n rr 
 
 since the expressions x . - x depend only on h, so does T^ T is regular, since 
 the system Tv = has only the trivial solution. 
 
 The particular choice of T we have got with the help of the polynomial 
 P (x ) gives to the iterative process (6), {l) } (8) two advantageous features: 
 
 l) Predictor formula . Because of (9)> P (x ) can be considered as 
 a polynomial of approximation of the true solution in the 
 neighborhood of x ; consequently P (x ) provides an approximate 
 value for y ; i.e., a predictor value: 
 
 xW \ (l) u t ( 2 ) (3) (k+j+lK 
 
 P(x ,) = v v +h(v v/ +v w +...+v v ° ); 
 n n+1 n n n n 7 
 
 2) Modification of the length of the step . Suppose we use the 
 
 step h for x < x and the step h* for x > x ; all the quantities 
 
 related to h* will be denoted with a star: x* = x , 
 
 n n' 
 
 — * — > 
 
 x* = x + h*, . . .; q, B, b also depend on h and will be 
 
 be replaced by q*, B*, b*; it is natural to choose v* such that 
 
 P*(x) = P (x), where 
 n v n v " 
 
 n n x n V h* / n V h* J " 
 
 by comparing with (10), one gets the simple rule: 
 
 (1)* (1) 
 
 v = v 
 
 (2)* (2) 
 
 v = v 
 
 (3)* h* (3) 
 
 n h n ' 
 
 (k+j+2)* .. /h*\k+j (k+j+2) 
 
 n Ui n ' 
 
 -25- 
 
by (9) it is easy to interpret the meaning of this rule in 
 terms of the original method (l); denoting by y*, y* , . .., 
 
 y* ., f*. .... f* the initial values necessary to the 
 
 n-j n' ' n-k J 
 
 integration with step h* from x , we get 
 
 y* . - P (x - ih*) i - 0,1, ..,,,}, 
 n-i nn ' ' ' ' ' 
 
 f* . - P'(x - ih*) i - 0,1, ...,k; 
 
 n-i n n ' ' 7 7 
 
 one can observe that generally f* . ^ f (x - ih*. y* . ) 
 
 n-i - n n-i 
 
 The Nordsieck Formulas 
 
 If (l) is the Adams-Moulton method described in Section 1 to which 
 one applies the transformation described above, then (6), (7), (8) represent the 
 Nordsieck formula. The formula corresponding to the k Adams-Moulton method 
 will be called "k Nordsieck formula;" the corresponding matrix B in (7) is 
 of order k+2 and the corresponding polynomial P (x) which is of degree k+1 is 
 defined by the relations 
 
 P (x ) = y 
 n n n 
 
 P'(x . ) m f ., i = 0,1, . ..,k 
 n n-1 n-i' ' 7 ' 
 
 The Modified Nordsieck Formulas 
 
 We can write the k Adams-Moulton formula in the following way: 
 
 y . = y + hp n f + hp . f . , . . + hft f , + Of . . + p . f . : 
 n+1 ''n n -1 n-1 k n-k n-k-1 -1 n+1' 
 
 we simply add a zero term; we can apply the same transformation on this formula 
 as we did for the Nordsieck formula; the result will be called the "k 
 modified Nordsieck formula;" the corresponding matrix B in (7) is of order k+3 
 and P (x) which is of order k+2 is determined by the relations 
 
 -26- 
 
P (x ) = y 
 
 n v n' J n 
 
 P»(x . ) = f ., i = 0,1, ...,k+l, 
 
 n n-i n-i 
 
 Remark 
 
 The k Nordsieck method and the k modified Nordsieck method are 
 
 both equivalent to the k Adams-Moulton method in the following sense: denoting 
 
 by P (x) the polynomial corresponding to the Nordsieck method and P*(x) the 
 
 polynomial corresponding to the modified Nordsieck method, then all three 
 
 methods will give the same values for y , n > if 
 
 n' — 
 
 P^x.i) - Pi*(*_ ± ) - t_ ± , 
 
 i = 0,1, . ,.,k; 
 
 the value of PA*(x , , ) has no influence on y , n > 0, 
 x -k-1' J n' — 
 
 Example 
 
 For k = 1, the Adams-Moulton method is given by 
 
 8h _h 5h 
 A n+1 y n + 12 n 12 n-1 + 12 n+l J 
 
 in order to get the Nordsieck formula, we have the sequence of intermediate 
 stages : 
 
 12' ~ 12 -1 " 12' 
 
 -T 
 
 u = (y , f , f .); 
 n KJ n' n' n-l y ' 
 
 A - 
 
 A 
 
 8h 
 12 
 
 h 
 12 
 
 
 
 
 
 
 
 \o 
 
 1 
 
 -27 
 
 12 ' 
 
 a = 
 
t> f \ d) J (2)/ " n \ (3)/ 
 
 P (x) = v + h iv — + v. ' 
 
 X - x \ /^/x - x\2 
 n \ ( 
 
 \ — Z — + v 
 n n n \ h ) n 
 
 -T ( (1) (2) (3). 
 
 n n ' n ' n 7 
 
 the relations P(x)=y.P'(x)=f.P'(x .) = f . lead to the relations 
 n n n 7 n n n 7 n n-1 n-1 
 
 where 
 
 T = 
 
 u n = Tv n (11) 
 
 
 
 
 °\ 
 
 
 /l 
 
 
 
 
 
 
 
 1 
 
 ol 
 
 T- 1 - 
 
 (° 
 
 1 
 
 
 
 p 
 
 1 
 
 J 
 
 
 Vo 
 
 1 
 
 2 
 
 1 
 " 2 
 
 -T -T_ ,. 7h hs 
 q = p T = (1, — , p 
 
 Th Ir 
 
 12 Z 
 
 ' X AT = ( ) , b = T _1 a 
 
 -I ° 
 
 (6) and (7) become explicitly 
 
 (l) 7h (2) h (3) 5h - 
 V = v + - — V + T V + - — f 
 ^n+l n 12 n 6 n 12 n+1 7 
 
 p (l) = y (l) ^ 7h „(2) ^ h „(3) ^ 5h 
 
 n+1 n 
 
 12 n 6 n 12 n+1' 
 
 n+1 n+1 7 
 
 r (3)..i T (a) + i, . 
 
 n+1 2 n 2 n+1' 
 
 -28- 
 
(1) 
 
 .(2) 
 
 (3) 
 
 as it appears also from (ll), v^ - y and v^ ' = f ; replacing uy* J by a , 
 one gets for the first Nordsieck formula 
 
 7h _ h 5h „ 
 y n + l = y n + 12 f n + 6 & n + 12 f n+l' 
 
 1 f 1 f 
 
 a n+l 2 n 2 n+1" 
 
 (12) 
 
 J 
 
 In the same way, we derive the modified Nordsieck formula for k = 1; 
 by writing the Adams -Moult on formula in the form 
 
 8h _ h _ __ 5h „ 
 y n + l = y n + 12 f n - 12 f n-l + ° f n-2 + 12 W 
 
 we get successively 
 
 -*T f . 8h h n v _ 5h 
 
 -KP 
 
 = (y , f , f , f ), 
 
 n n' n n-1' n-2 ' 
 
 A = 
 
 8h 
 12 
 
 h 
 "12 
 
 
 
 
 
 1 
 
 
 
 
 
 1 
 
 a = 
 
 the polynomial P (x) = v' 1 ' + h(v' 2 M — - — -J + v' 3 M — =- — -) + v 
 
 n n v nVh/ n \ h / nVh 
 
 (U^L^V) 
 
 satisfies the conditions P (x ) = y . P'(x .) = f .. i = 0.1,2; it follows 
 
 n n n' n n-i n-i' ' ' y 
 
 that 
 
 u = Tv 
 n n 
 
 (13) 
 
 where 
 
 -29- 
 
T = 
 
 1 
 
 
 
 
 
 °\ 
 
 
 
 1 
 
 
 
 \ 
 
 
 
 
 
 1 
 
 -2 
 
 3 / 
 
 
 
 1 
 
 -k 
 
 */ 
 
 T 
 
 
 
 
 
 
 
 1 
 
 
 
 
 
 3 
 
 
 1 
 
 k 
 
 -1 
 
 h 
 
 1 
 
 1 
 
 l 
 
 Z 
 
 " 3 
 
 6 
 
 -»T -*T m /., 7h h h- 
 
 B = T AT= 
 
 q 
 
 - p i 
 
 - K±j 
 
 12' 
 
 IS' " 
 
 " k' 
 
 i 
 
 
 7h 
 
 h 
 
 £\ 
 
 
 
 
 
 M 
 
 12 
 
 6 
 
 + \ 
 
 \ 
 
 
 
 
 /l2\ 
 
 
 
 
 
 
 
 1 
 
 b = 
 
 -- t" 
 
 l-» 
 
 a = 
 
 " 
 
 3 
 
 l 
 
 3 
 
 
 
 
 
 1 } J 
 
 h 
 
 2 
 
 5 
 
 / 
 
 
 
 
 \ * / 
 
 1 
 
 1 
 
 i / 
 
 r 
 
 
 
 
 \ i/ 
 
 6 
 
 " 3 
 
 2/ 
 
 
 
 
 
 W 
 
 (6) and (7) become explicitly 
 
 (l) 7h (2 h (3) 
 
 J n+1 n 12 n 6 n 
 
 5 V n + 12 n+1' 
 
 (1) (1) 7h (2) h (3) 
 
 n+1 n 12 n on 
 
 5 V n + 12 n+1' 
 
 v^ = f , 
 
 n+1 n+l ; 
 
 r (3)_ 3 v (2) 1 (3) + 3 v 0O + 3 f 
 
 T _ — - T- V - — V +T-V +T-I _,. 
 
 n+1 4 n 2 n 4 n k- n+1' 
 
 (M 1 (2) 1 (3) 1(h) 1 . 
 
 n+1 6 n 3 n 2 n 6 n+1' 
 
 setting v^ = a and v^ = b . we get for the first modified Nordsieck method: 
 n n n n' 
 
 7h „ h h , 5h _ 
 
 y . = y + — f +7^a -rt +f--f , . 
 J n+1 J n 12 n 6 n + n 12 n+1' 
 
 -30- 
 
3 13 3 
 
 a = - f f - ± a +r b + r f u 
 n+1 4 n 2 n 4 n 4 n+1' 
 
 b , = - ? f - -r a + 77 b + z f . 
 
 n+1 on 3 n 2 n 6 n+1 
 
 Another Way for Deriving the Nordsieck Formulas 
 
 A comparison of (12) with Reference 1 shows that Nordsieck gets its 
 formula in a slightly different form. In fact, by working directly on the 
 polynomial P (x), one obtains the "true" Nordsieck formula much faster. 
 
 a) Nordsieck formula 
 
 Let p < q be integers and S(x:p,q) be the Lagrangian 
 polynomial of degree (q - p): 
 
 (x - x )(x - x , ) . . . (x - x , ) 
 
 p /v P+l _ 3" 1 
 
 q 
 
 th 
 
 s(x:p ' q) = (x - x ')(x - x''") ••• 5 - * 75 ' 
 
 q P q p+l q q-1 
 
 the polynomial P (x) of degree k+1 for the k Nordsieck formula 
 is defined by the relations 
 
 p (x ) = y , (13) 
 
 P'(x ) = f ., i = 0,1, ...,k; (Ik) 
 
 n v n n-i' ' ' ' 7 
 
 P ,(x) will be determined by 
 
 P (x , ) = y , . 
 
 n v n+l y J n+1' 
 
 P' (x . . ) = f _ ., i = 0,1, ...,k; (15) 
 
 n+1 n+l-i n+l-i' ; ' ' 7 ' 
 
 P'(x) and P' (x), polynomial of degree k, are entirely 
 determined by (ik) and (15); one checks easily that 
 
 -31- 
 
P' _(x) = P'(x) + (f _ - P'(x . ))s(x:n + 1 - k. n + l): (l6) 
 
 n+1 n n+1 n n+1' v ' J ' v ' 
 
 let f P = P*(x . ); by (10) 
 n n+1 ' 
 
 f P . y (2) + 2v (3) + (k) ( k + l)v (^). 
 
 n n n v n ' 
 
 equalling terms of some degree in (l6) allows to compute 
 
 (2) (k+2) . . .. .. _ (2) (k+2) 
 
 v ,, .... v' as linear combination of v ..... v ' 
 
 n+1' ' n+1 n ' ' n 
 
 and (f _ - f P ). 
 K n+1 
 
 It remains to compute v^ _. = P . (x . ) = y , ; by 
 
 n+1 n+l v n+l' "'n+l' J 
 
 Section 1, y . = Q(x ), where Q(x) is the polynomial of 
 
 degree k+2 satisfying the relations 
 
 Q(x ) = y . 
 
 n n 
 
 Q'(x .) = f ., i = -1,0, . ..,k; 
 
 n-i n-i ' ' ' * 
 
 one checks easily that 
 
 Q'(x) = P'(x) + (f . - P'(x . ))s(x:n - k, n + l), 
 
 Q(x) = P n (x) + (f n+1 - P;(x n+1 ))J S(y:n - k, n + l)dy; 
 
 x 
 
 n 
 
 consequently 
 
 - „(1) _ 
 
 x , 
 n+1 
 
 W ■ v £i P n< x n+ 1> + (f n + l " fP) j S(y:n " k ' n + l)dy (17) 
 
 X 
 
 n 
 
 ■32- 
 
b) Modified Nordsieck formula 
 
 The polynomials P (x ) and P , (x ) of degree k+2 for the 
 * J n v n+l v ° 
 
 "th 
 k modified method are determined by the relations 
 
 P (x ) = y , P'(x .) = f ., i = 0,1, . ..,k+l; 
 
 n x n n n n-i n-i' ' ' ' ' 
 
 P . (x ,) = y ,, P' (x _ .) = f . ., i = 0,1, ...,k+l: 
 
 n+1 n+1 J n+1' n+l v n+l-i' n+l-i' ' ' ' ' 
 
 then 
 
 P; +1 (x) = P;(x) + (f n+] _ - f P )s(x:n - k, n + l), (l8) 
 
 with f P = P'(x .) = v (2) + 2v (3) + 3^ k) + ... + (k + 2)v (k+3) i 
 n n+1 n n n n * 
 
 equalling the coefficients of some degree in (l8) gives 
 
 (2) (k+3) . . .. .. _ (2) (k+3) 
 
 v ',..., v -, as linear combination of v , .... v and 
 
 n+1 7 ' n+1 n ' ' n 
 
 (f . - f P ). 
 v n+1 
 
 By Section 1, y = Q(x ), where Q(x) is the polynomial 
 
 of degree k+2 which satisfies the relations: 
 
 Q(x ) = y , Q'(x .) = f ., i = -1,0,. ..,k; 
 n J n n-i n-i' ' ' ' ' 
 
 then 
 
 Q'(x) = P'(x) + (f . - P'(x . ))s(x:n - k, n + l); 
 
 a comparison with (l8) shows that Q(x ) = P -.(x), since 
 
 Q(x , ) = y , = P , (x . ) : 
 v n+1 7 J n+1 n+l v n+1 
 
 P_-,(x) ■■■■ Q(x) - Pjx) + (f_ n P:(x_J)J S(x:n - k, n + l)dx (19) 
 
 x 
 
 X n + 1 
 y n+ l = F J\^ + < f „ + l - f P )/s(x:n - k, n + l)dx. (20) 
 
 X 
 
 n 
 -33- 
 
Remark 
 
 A conparison of (l6) with (l8) and (17 ) with (19) shows that the k 
 
 modified Nordsieck formula is identical to the (k+l) Nordsieck formula, except 
 
 for the coefficient of (f , - f ) in the expression of y .. which is the same 
 , n+1 'n+l 
 
 in both k formulas . 
 
 Example . We consider again the case k = 1. 
 Nordsieck formula 
 
 p (x) . r* 1 ' ♦ h^'hr 5 ♦ ^'{ __£ ); 
 
 n n n \ h / n V A / 
 
 (1) ., (2)/ X " X n \ (3) 
 
 = v + h(v v : + v A 
 
 n n V h J n 
 
 by (13) and (14), v^ = y , v' 2 '' = f : we set a = v^ 3 '; by (16) 
 
 n' n n n n 
 
 x-x x-x x-x 
 
 P' (x) = f . + 2a . r2i± = f + 2a -^—2 + (f - f p ) — — » 
 
 n+1 n+1 n+1 h n n h n+1 y h 
 
 it follows that a . = a + -(f _ - f P ): by (17): 
 n+1 n 2 V n+1 " 
 
 x 
 
 p (y - x )(y - x ) 
 
 y . = P (x _ ) + (f . - f P ) / ^~ — dy, 
 
 J n+1 n^ n+1 7 v n+1 'J n1 2 J ' 
 
 2h 
 x 
 n 
 
 f 
 
 y , = y + h(f + a ) + |£ (f , - f P ); 
 J n+1 J n v n n' 12 v n+1 " 
 
 we get for the first Nordsieck formula the expressions 
 
 y _ = y + h(f + a + t| (f , - f P )), 
 J n+1 J n v n n 12 v n+1 '" 
 
 f P = f + 2a , 
 n n 
 
 a = a + \ (f . - f P ) 
 
 n+1 n 2 v n+1 
 
 •3+- 
 
Modified Nordsieck formula 
 
 As before, we substitute: v = y » v^ ' = f , v = a . 
 
 n rr n rr n n 
 
 v^ ' = b . By (19): 
 
 n n 
 
 /x-x \ /x-x \2 /x-x \3 r (y-x ., )(y-x ) 
 
 ■ ^ n ff % -r *n(-T s ) ^(^ 1 -* p ) J ""Jg P *■ 
 
 ^ ' x 
 
 n 
 
 by taking the successive derivative at x for the two members^ 
 one gets 
 
 where 
 
 y = y +h(f +a +b + ^ (f - f P )), 
 n+1 n n n n 12 v n+1 
 
 = 2a^ + 3b + 2 (f _ - f P ), 
 
 n+1 n n 4 n+1 
 
 n+1 n n+1 " 
 
 f P = f + 2a + 3b . 
 n n n 
 
 -35- 
 
SECTION 5. PREDICTOR FORMULAS 
 
 a) Nordsieck formula 
 
 For the k Nordsieck formula, the predictor value for x _, is given 
 
 ' n+1 & 
 
 by P (x ), where P (x) is the polynomial of degree k+1 determined by the 
 relations 
 
 P (x ) = y , P (x .) = f .. i = 0,1, ...,k; 
 n n n' n n-i n-i * ' ' 
 
 a conparison with Section 1 shows that P (x ) is nothing else then the value 
 
 th n k+2 (k+2) 
 
 of the k Adams -Bashforth formula with a truncation error D h y (!)• 
 
 b) Modified Nordsieck formula 
 
 th 
 For the k modified Nordsieck formula, the predictor value for 
 
 x is given by P (x ), where P (x) is the polynomial of degree k+2 determined 
 
 by the relations 
 
 P (x ) = y , P*(x .) = f ., i = 0,1, . ..,k + 1; 
 
 n n n' n n-i n-i 7 ' ' '■ 7 
 
 it follows that P (x ) is the value of the (k+l) Adams -Bashforth formula 
 
 with a truncation error D n h y (l)i 
 
 k+l 
 
 Remark 
 
 The modified formula has the advantage that the predictor and the 
 corrector formulas have a truncation error of the same order. 
 
 • 36- 
 
SECTION 6. MODIFICATION OF THE LENGTH OF THE STEP 
 
 Suppose we integrate the differential equation y' = f(x,y) from a to 
 c with the step h for x < b and h* for x > b. 
 
 Let u(x) be the discrete approximative solution for x < b and u*(x) 
 be the approximative solution for x > b. Let z (x ) be the exact solution of the 
 differential equation such that z(b) = u(b); as h varies, z(x) varies also. 
 
 We use the notations 
 
 x = b + nh, x* = b + nh*: 
 
 n ' n 
 
 u = u(x ), u* = u*(x*) 
 n n ' n v n • 
 
 f = f (x ,u ) for n < 0: 
 n n' n — 
 
 f* = f(x* u*) for n > 0: 
 
 n ^ n' n — 7 
 
 z = z(x ), z* - z(x*) 
 n n n n 
 
 We suppose 
 
 1) p = — ■ is a constant as h •+ 0: 
 
 h* 
 
 2) all the hypotheses necessary to apply the self -stabilizing 
 property and the generalized theorem in Section 2 are 
 filled. 
 
 Because of the self-stabilizing property, it is suitable to compare 
 u*(x) with z(x). Let d*(x) = u*(x) - z(x) for x > b. Replacing d(x) by d*(x), 
 h by h*, q by 0, the generalized theorem gives the asymptotical relation 
 
 d*(x) = h^( ei (x) + 0(h)) + h* S (e IZ (x) + 0(h)) + h r+1 (e 1T1 (x) + 0(h)), x > b; 
 
 the term e (x ) corresponds to the accumulation of the local truncation errors 
 and is comsequently independent of the starting procedure for the new h*; the 
 
 -37- 
 
term e (x ) corresponds to the error at x = b; by hypothesis d*(b) =0, i.e., 
 s = _oo and the second term of the asymptotical relation vanishes; only the 
 term e (x ) depends on the starting procedure. 
 
 Nordsieck Method 
 
 Consider the k Nordsieck formula; according to Section k, P^(x) is 
 the polynomial of degree k determined by the relations 
 
 P^(x_.) = t , i = 0,1, ...,k; 
 
 the starting values f*. for the new h* are given by 
 
 f ?i = P c,( x ?i)> i = 0,1, ...,k; 
 
 because of well-known properties of the interpolation there exists the constant 
 coefficients (independent of h*) a. . such that 
 
 f*. = E a. .f ,, i = 0,1,. .,,k, (l) 
 
 -i , =0 iJ -J 
 
 z'* = Z a. .z' . + R., i = 0,1, . . .,k, (2) 
 
 -1 ij -J x 1 
 
 where the interpolation error R. is given by 
 
 (*;, - *p )(*?!- mj ••• ^!j - *_ k ) (k+2 ) 
 \-- ( k ; i): z bi.*' i = o,i,...,k 
 
 with min(x*.,x ) < | < x . Since h = ph* and z^ k+2 ^(| . ) = z^ k+2 '(b) + 0(h), 
 we can write 
 
 R . . h *k + l (-i)(p - i)(gp ■ O ... (kp- 1) z (k + 2) (b) + 0(h# k + 2 ); . . 0jl; ... jk . 
 
 subtracting (2) from (l), one gets 
 
 -38- 
 
f*. - z 1 * = -R. + 0(h* ), 
 -1-1 1 v ' 
 
 since by the self-stabilizing property f . - z' . = 0(h* ) 
 
 J J 
 
 Using the notation 
 
 . (-i)(p - l)(2p - i) ... (kp - i) 
 \ ~ (k + i): 
 
 we get by the generalized asymptotical theorem 
 
 d*(x) = h^ k+2 ( ei (x) + 0(h)) + h* k+2 ( e;[II (x) + 0(h)), b < x < c (3) 
 
 where 
 
 ^ e x (x) = G(x, Z (x))e I (x) - C k z (k+3) (x), e^b) = 0, 
 
 dx e HI ( ' X ^ = G ( x > z ( x )) e iii( x )> 
 
 e Ii;r (b) = - {(p x + P 2 + ...+ P k )E x + (0 2 + P 3 +...+ P k )E 2 +...+ P k E k )z (k+2) (b), 
 
 th 
 we recall that |3 , f3 , P , . . . , p are the coefficients of the k Adams- 
 
 ~ _L U _l_ K. 
 
 Moult on method. 
 
 Remarks 
 
 1) The error due to the modification of the step is non-negligible in comparison 
 with the error due to the accumulation of the truncation errors. 
 
 2) In (3), G(x, z(x)) is a scalar," for systems of differential equations, it 
 suffices to replace d*(x), e (x), e (x), z(x) by the corresponding vectors 
 and G(x,z(x)) by a matrix. 
 
 ■39- 
 
Modified Nordsieck Formula 
 
 Consider the k modified formula; according to Section k, PJL(x) 
 is a polynomial of degree k+1 determined by the relations 
 
 Pq(x_.) = f_ ± , ± = 0,1,.,,, k+1; 
 
 the initial values f*. for the new step h* satisfy the relations 
 
 f !i = F ^ x -^' i = 0,1,. .„k; 
 
 there exist the constant coefficients b . . (independent of h*) such that 
 
 k+1 
 
 f* = Z b,,f ,, i = o,i,,.,,k, (k) 
 
 j=o 
 
 k+1 
 
 z'* = Z b. .z' . + r., i = 0,1, ...,k, (5) 
 
 -i Q ij -j i' ' 
 
 where R. is the interpolation error 
 
 „ (X -1 - X )(X -1 - X l } ••• (X -1 - X -k-l } (k + 3)„ , , „ , v 
 
 with min(x*., x . , ) < 6 . < x_. We can write 
 -i' -k-1 — " i — 
 
 R. = h* k+2 H.z (k+3) (b) + 0(h* k+3 ), 
 
 II. 
 
 - (-i)(P - i) ... ((k + l)p - i) 
 
 1 (k + 2).' 
 
 because of the self-stabilizing property, if we subtract (5) from (h) t one gets 
 
 f* - z'* = -R. + 0(h* h+3 %- 
 -l -l l v ' 
 
 ■+o- 
 
by the generalized asymptotical theorem, one obtains the relation 
 
 where 
 
 d*(x) = h* k+2 ( ei (x) + 0(h)) + h* k+3 (e ni (x) + 0(h)) (6) 
 
 gj ei (x) = G(x,z(x)) - C k z (k+3) (x), e x (b) = 0; 
 
 dx e HI^ X ^ = G ( x ^ z ( x )) e III ( x )^ 
 
 e IT1 (h) = - {(P 1 + 3 2 + . .. + P k )H 1 + (P 2 + ... + P k )H 2 + ... + p^Jz^^^b); 
 
 th 
 we recall that 8 _ , 8_. 8. . .... 6 n are the coefficients of the k Adams-Moulton 
 
 -1' 0' 1' k 
 
 formula . 
 
 Remarks 
 
 1) Except for x closed to b, the error due to the modification of the step is 
 negligible in (6); 
 
 2) the same relation (6) holds for systems of differential equations if one 
 replaces the scalars by the convenient vectors or matrix. 
 
 -kl- 
 
SECTION 7- STARTING PROCEDURES 
 
 Although the starting procedure is not the most characteristic feature 
 of the Nordsieck method, it is most painful to analyze. In order to make the 
 job easier, we shall introduce a so-called "simplified starting procedure." 
 
 Given a differential equation 
 
 y' = f(x,y) (1) 
 
 and initial value 
 
 the problem is to find a "good" initial v~ (or equivalently an initial poly- 
 nomial P (x) ) . 
 
 We introduce the notations : 
 
 x = b + nh, f = f(x .y ), n > 0; 
 
 the norm | |x| of any vector x is the sum of the absolute values of its components 
 We suppose : 
 
 1) f(x,y) satisfies the Lipschitz condition |f(x,y*) - f(x,y)| < ft|y* - yrj., 
 b < x < c; 
 
 2) the solution z(x) of (l) with z(b) = y_, possesses (k+k) continuous 
 derivatives. 
 
 Starting Procedure for the k Nordsieck Method 
 
 Starting with any polynomial P n ( x ) such that P n ( x n ) = y n and 
 P' n (x n ) = f n , k steps of integration are executed forwards; then h is replaced 
 by -h (by the standard procedure of modification of the step) and k steps are 
 executed backwards; -h is replaced by h; let Q n ( x ) he the polynomial obtained 
 at that stage; the new polynomial P ,(x) is determined by the conditions 
 
 \J 9 -L 
 
 -42- 
 
p o,i<*o> = V p i,i (x o> " f o> p o,i( x) ■ 9 0,0< X '' 
 
 the same procedure is repeated with P -,(x); one gets a sequence of polynomials 
 
 p o,o (x) < V (x) ' p o,i (x) ' %,i M > ■■'■■ 
 
 th 
 Simplified Starting Procedure for the k Nordsieck Method 
 
 Starting with any polynomial P n (x) such that P n n (x n ) = y n and 
 
 P' _(x_) = f_, k steps of integration are executed forwards; let P. ~(x) "be 
 0,0 Cr ^ ° ' k, X 
 
 the polynomial obtained at this stage; P ,(x) is then defined by 
 
 V< x o> ■ V V (x) ■ p k,o( x) ' 
 
 the same procedure is repeated for the new initial polynomial P n , (x); a sequence 
 
 p o,o (x) ' P k,o (x) > p o,i (x) > P k,i (x) > •■• is defined ' 
 
 For the Nordsieck formulas , we have the following theorem. 
 
 Theorem 1 
 
 For h sufficiently small, all the polynomials, P_ . (x) and Q .(x) 
 for the starting procedure, P . (x ) and P . (x ) for the simplified starting 
 
 U, 1 K, 1 
 
 procedure, coverge as i goes to infinity to the same polynomial U(x) of degree 
 k+1 satisfying the relations 
 
 U(xJ = y n , (2) 
 
 U'(x.) = f(x.,U(x.)), i = 0,1,. ..,k. (3) 
 
 th 
 Starting Procedure for the k Modified Nordsieck Formula 
 
 P_ . (x ) and CL . (x ) are defined in the same way as for the Nordsieck 
 0,i 0,i J 
 
 formula, except that (k+l) steps of integrations are executed forwards and 
 backwards at each iteration. 
 
 43- 
 
Simplified Starting Procedure for the k Modified Nordsieck Formula 
 
 It is identical to the simplified starting procedure for the k 
 
 Nordsieck method except that (k+l) steps of integration are executed at each 
 
 iteration and define the sequence P. ~(x), P. . „(x), P^ , (x), P n . n (x). ... 
 
 0,Cr '' k+l,O v " 0,l v " k+1,1 
 
 For the modified Nordsieck formula, we have the following theorem. 
 
 Theorem 2 
 
 For h small enough, all the polynomials, P . (x ) and Q . (x) for the 
 starting procedure, P . (x) and P . (x) for the simplified starting procedure, 
 converge as i goes to infinity to the same and unique polynomial U(x) of degree 
 k+2 satisfying the relations 
 
 U(0 - Yn> 
 
 U'(x.) = f(x.,U(x.)), i = 0,l,...,k+l, 
 
 Proof of the Convergence of the Starting Procedures 
 
 We shall consider only the case of the simplified starting procedure 
 for the k Nordsieck method; all the other cases can be treated analogously. 
 
 We first recall that v is the vector, the components of which are 
 the coefficients of the polynomial 
 
 p , y (1) h , (2){ X ~ X n\ (k+2)^_l^n\ k+ \ 
 
 P (x) = V v + h(V : + ... + V : ) 
 
 n n v n\h/ n \h/ 
 
 and that u is the vector of the components y,f,f ,.,.., f , . Since 
 
 n * J n } n } n-1' ' n-k 
 
 P (x ) = y , P 1 (x . ) = f . , i = 0,1, . . . ,k, there exists a regular matrix T 
 n n "'n' n v n-i n-i' ' ' 
 
 such that u = Tv . Let u . and v , be the components of P . (x ) 
 n n n, l n, k n, l 
 
 (n = 0, 1, ...,k); proving the convergence of P .(x), (or P . (x)) for i -» °° 
 
 U . 1 K « X 
 
 — > -> ' 
 
 is equivalent to the proof of the convergence of v n . or u n . . 
 
 —» —* .—* . 
 
 Let for any u_, k(u ; be the new initial vector obtained after one 
 
 iteration of the simplified starting procedure. From Appendix 1, it follows 
 
 -kk- 
 
that the convergence of u_ . towards a limit vector independent of the initial 
 vector u n _ will be insured if we can show that there exists a constant K < 1 
 such that for any two u and u* (u!l u q Yq; u q = u o f o^ we have 
 
 k(u*) - k(u q )\ I < K||u* - u Q | 
 
 We apply to u and u* the simplified starting procedure; the quantities con- 
 
 — > 
 cerning u* will be denoted by a star. Let 
 
 y , = y + h(B _f , + ...+ ft f . ) 
 J n+1 J n v -1 n+1 K k n-k y 
 
 be the Adams -Moult on formula; for y and y*, we get 
 
 y i = y o + h ^-i f i + p o f o + "• + V-k*' 
 yJ = vg + h(p_ 1 f* + :P f S+ ... + P k f ? k ); 
 
 but |f* - f | < ft|y* - y |; if 1 - |p |hfl > 0, then 
 
 lv* - y I < t— — i — ( I v* - v I + h ( 1 8 I I f * - f I + , . . + 1 8 If* -f \))i 
 
 u l y l l - 1 - |B |hn uy y l a vlP l l x o 1 l |P k' ' -k -k Uh 
 
 from now on. we suppose that h < h = -rr-r — -r— . let B = 6_ + \b.\ + ... + \&, \ ; 
 
 — 2 1 6 1 * ' ^0 ' "i. ' k ' ' 
 
 since y* = y _, we have 
 
 |y* - yj < 2hB| |u* - u Q | | and | f * - f ± \ < 2hftB| |u* - u Q | |; (k) 
 
 in the same way, we get for the second step: 
 
 |y* - v I < 2fl v* - v + h( Ib I f * - f + I B I If * - f I + + I B I If* -f 1)1 
 u 2 • y 2 l - u,y l y l' "^IPo 111 ! I i" + -l p l» 1*0 I I IP k M 1-k l-k" 
 
 by (h): 
 
 |y* - y 2 l < 2{2hB-||u* - u || + h(2hBa||uJ - u Q | | + b||u* - ujl)}; 
 
 -4 5 - 
 
\y* - y 2 | < Bh(3 + 2hn)| |u* - u Q \\; 
 
 since B and ft are constant, we can write 
 
 y* - y | : 0(h) ||u* 
 
 u.. 
 
 2 j^ _ v-/i !- 
 
 t 
 
 repeating the same process k times, we finally get 
 
 |y* - y.| = 0(h)| |u* - u Q | |, i = 0,1,,..,*, (5) 
 
 |f J - f.| = 0(h)| |ug - u Q | |, i = 0,1,. ..,k. (6) 
 
 Let R(x) and R*(x) be respectively the polynomials of degree k such that 
 
 R(x.) = f ., i = 0,1, ...,k, 
 
 R*(x.) = f *, i - 0,1, ...,k; 
 
 P l' '"' P k+2' p l' '••' P k+2' com P° nents respectively of k(u q ) and k(u*) are 
 given by 
 
 p x = pJ = y , 
 
 p 2 = R(x Q ) = p* = R*(x Q ) = f Q , 
 
 P 3 = R(x_ 1 ); p* = R*(x_ 3 ), 
 
 P k+2 ■ R(x )j P J +2 = R*(x_ k ); 
 
 then, there exist the constant coefficients (independent of h) such that 
 
 46- 
 
p. = Z a. .f ., 
 
 1 j=0 1J J 
 
 k 
 p* = Z a. .f*: 
 
 1 j=o 1J J 
 
 by (6) 
 
 p. - p*| = 0(h)||u* - u |.|, i = 0,1,. ..,k 
 
 |R(u*) - k(u q )| I = 0(h)| |ug - u Q | 
 
 this is the desired result since for h small enough 0(h) < 1. 
 
 Remark that we have proved so far only the convergence of the 
 polynomials P n . (x); "but the same arguments apply also for P, . (x)„ 
 
 Proof of Theorem 1 
 
 We have already proved the convergence of the simplified starting 
 procedure and of the starting procedure. 
 
 Let us call U(x) the limit polynomial for the simplified starting 
 procedure; by construction U(x) satisfies (2) and (3) for i = 0; it remains 
 to prove (3) for i = 1,2, ...,k. Let y , y , ..., y be the values obtained by 
 
 , , 12 K. 
 
 integrating with the k Nordsieck formula and starting polynomial U(x). 
 Because U(x) is the limit polynomial, it follows that 
 
 U'(x.) = f(x.,y.), i = 0,l,...,k; (7) 
 
 since the k Nordsieck formula is equivalent to the k Adams -Moulton formula, 
 it follows that y = R(x ), where R(x) is the polynomial of the degree k+2 
 satisfying the relations 
 
 R(* ) = y = u( x ^ 
 
 47- 
 
R'(x.) = U'(x.), 0,-1,. ..,-k, 
 
 R'(x 1 ) = f( Xl , yi ); 
 
 ^y (T); R'OO = U'(x ), and this implies that R(x) = U(x), y = U(x ), which 
 proves (3) for 1=1, For i = 2,3* . ..>k the proof follows the same pattern. 
 
 We now discuss briefly the part of theorem 1 which concerns the 
 starting procedure. Let P n (x) and Q n (x) be the limit of the polynomial P .(x) 
 and Q (x). We know that P n (x) and Q n (x) are independent of P n ( x ); let us 
 choose P n (x) = U(x), where U(x) is the polynomial which satisfies the relations 
 (2) and (3) (its existence has been established above.)/ by arguments similar to 
 those used previously, it follows easily that Q_ n (x) = U(x) = P n n (x) and 
 consequently Pq( x ) = ^n^ x ^ = u ( x )- 
 
 In order to avoid fastidious repetitions, we renounce to give the 
 details of the proof of theorem 2. 
 
 Starting Errors in the Nordsieck Method 
 
 From theorem 1, it follows that the limit initial polynomial U(x) is 
 the same for the starting and the modified starting procedure. Let P (x), 
 P (x), ... be the polynomials corresponding to the steps 0,1,2,.../ we know that 
 U(x) = P (x) = P k (x) (and, in fact, also U(x) = P (x ), P (x), . . . ,P ,(x)). Let 
 z(x) be the exact solution of the differential equation with z(x ) = y ; we want 
 to apply the asympototical theorem of Section 2 with q = k/ we have to determine 
 
 U(x k ) - z(x k ), (8) 
 
 U'(x.) - z'(x.), i = 0,'l,...,k. (9) 
 
 Let Q(x) be the polynomial of degree k+1 (as U(x)) satisfying the relations 
 
 Q(x ) = y , (10) 
 
 Q'(x.) = z'(x.), i = 0,1,. ..,k; (11) 
 
 ■kQ- 
 
let 
 
 R. = z(x.) - Q(x.), i = 0,1, ...,k; (12) 
 
 we have the relations 
 
 U(x Q ) = y 
 
 U'(x.) = f(x.,U(x.)), i = 0,1, ...,k, (13) 
 
 Q(x Q ) = y , 
 
 »' (x.) = f(x.,Q(x.) + R.), i = 0,1, ...,k; (lV 
 
 let 
 
 H(x) = U(x) = Q(x), (15) 
 
 t. - H'(x.), i = 0,1, ...,k; (16) 
 
 since H(x n ) = 0, there exist the constants a., independent of h such that 
 
 k 
 H(x) = h Za t, 1 = 0,1,...,^ (17) 
 
 j=0 1J J 
 
 we define the quantities g., i = 0,1, ...,k by the relations 
 
 f(x.,U(x.)) - f(x.,Q(x.) + R.) = g.(u(x.) - Q(x.) - R.), i = 0,1,. ..,k; 
 
 (18) 
 
 if U(x. ) - Q(x. ) - R. = 0, we set g. = 0, . 
 
 i ^ i 7 i * & i 
 
 We know that |g.| < Q; subtracting (lh) from (13), we get by (15), 
 (16), (17), (18), 
 
 -k 9 - 
 
k 
 
 t. = g.h 2 a tj - g.R., i = 0,1, ...,k; (19) 
 
 J — \J 
 
 this is a linear system for the t.'s if we consider R. as known; we can regard 
 t. as function of h; suppose that there exists a constant m such that 
 
 R. = 0(h ), i = 0,1, ...,k; it follows the relati 
 
 on 
 
 1 
 
 t. = -g.R. + 0(h m+1 ). (20) 
 
 111 v ' 
 
 We only sketch the proof of (20 ). Since g. is bounded, a. . independent of h, 
 
 (20) follows from (19), if we show that t. = 0(h m ), i = 0,1, ...,m. Let t be the 
 
 -> 
 vector of components t., w the vector of components -g.R. and B the matrix of 
 
 elements g.a..; (19) becomes 
 1 ij 7 
 
 (I - hB)t = w; (21) 
 
 we introduce the following notations: for any matrix C, the norm | |c| | of C 
 is the sum of the absolute value of its elements; norms of vectors and matrices 
 have the properties: 
 
 a) I |Cx| I < I |c| I I |x| |; 
 
 b) I |c + D| I < I |c| I + I |d| I; 
 
 c) ||cd|| < ||c|| • ||d|| ; 
 
 by the Gerschgorin theorem, if h < 1 11 1 , all the eigenvalues of hB lie in the 
 
 — 2 I I B \ I 
 
 circle of radius 0.5 around the origin; the matrix I - hB is regular and its 
 inverse possesses the convergent development 
 
 (I - hB)" 1 = EhV; 
 n=0 
 
 by b ) and c ) 
 
 00 
 
 n=0 
 -50 
 
by (21) and a) 
 
 — x - n ±5 — 
 
 Now we can determine the errors (8) and (9); "by (12) and (15), we have 
 
 U(x.) - z(x.) = H(x.) - R., i = 0,1,... ,k; 
 
 by (17) and (20 ) 
 
 m+2 
 
 U(x ) - z(x ) = -h Ea g R - R + 0(h ), i = 0,1, ...,k; (22) 
 
 by (ll), (16) and (20 ) 
 
 U'(x.) - z'Oc.) = t = -g.R. + 0(h m+1 ), i = 0,1, ...,k. (23) 
 
 In order to get more useful expressions, let us suppose (as in the asymptotical 
 
 df 
 theorem) that G(x,y) = — exists and satisfies the following requirements: 
 
 l) |G(x,y*) - G(x,y)| < ty|y* - y|, \|/ constant, b < x < c; 
 
 2) G(x,z(x)) possesses a continuous derivative j 
 
 it follows that g. = G(x.,y*), where y*€[U(x.), z(x.)]; by (22), U(x, 
 0(h ),• because of the above two requirements, we have 
 
 z(x.) 
 
 g. - G(x.,z(x.)) + 0(h m ) - G(x Q ,z(x )) + 0(h) + 0(h m ) 
 
 Replacing in (22) and (23), we get the final formulas 
 
 *\ 
 
 U(x.) - z(x.) = 
 
 U'(x.) - z'(x.) - 
 
 1 1 
 
 ■hG(x n ,z(x n )) Ea R - R + 0(h m+2 ), 
 j=0 J J 
 
 0' v 0' 
 
 > i. - 0,1, . . .,k 
 
 ■G(x ,z(x ))R. + 0(h m+1 ); 
 
 J 
 
 (2k) 
 
 -51- 
 
 UNIVERSITY OF 
 11 1 imdiq I IRRARY 
 
We recall that 
 
 G(x,y) - 21^1 
 
 Q(x) is the polynomial of degree k+1 which satisfies the relations 
 
 Q(x Q ) = z Q , Q'(x.) = z'(x.), i = 0,1,.. .,k; 
 
 R. = z(x.) - Q(x.); 
 
 a. . are the coefficients of the relations 
 
 K(x ) - K(x ) + h Za K'(x ), i = 0,l, ...,k, (25) 
 
 j=0 J 
 
 where K(x) is any polynomial of degree k+1, 
 
 Remark 
 
 We have considered the errors of the starting procedure for the limit 
 L U(x); 
 some j; we write 
 
 polynomial U(x); but practically, we must deal with the polynomial P . for 
 
 P. (x. ) - z(x ) = (P (x ) - U(x )) + (U(x ) - z(x )), i = 0,1,. ..,k; 
 
 by choosing P n (x) = y n + f n (x - x„), one gets easily from the proof of the 
 convergence that 
 
 P (x ) - U(x ) = 0(h 2+J ), i = 0,1, ...,k. (26) 
 
 This result is valid for both the starting procedure and the modified starting 
 procedure; for the simplified starting procedure, there is no reason for thinking 
 that the order of the error in (26) may be greater than h ,° however, it is 
 very likely that for the starting procedure we have 
 
 -52- 
 
P (x ) - U(x ) - 0(h 2+2j ), i = 0,1, ...,k; 
 
 no proof of this statement has "been established; if that is true, it means 
 that the speed of convergence in comparison with the number of steps of 
 integration to execute is the same for the starting and the modified starting 
 procedure; the simplified starting procedure has the advantage of using only 
 forward integration which remains closer to the true solution than the backward 
 integration. 
 
 Example 
 
 We consider the starting procedure for the first Nordsieck method and 
 
 suppose 
 
 1) enough iterations have been executed for using (2k) } when U(x) 
 
 is replaced by P .(x), 
 V> J 
 
 2) all the assumptions necessary to the application of the 
 asymptotic theorem are filled. A little arithmetic leads to 
 the relations 
 
 A - 1 
 
 a 00 " a 01 " ° ; a i0 " a il " 2 ; 
 
 R = 0; R x = - X2 z '"(* ) + 0(h^.); 
 
 h 3 „,,h 
 
 y 
 
 P 0> .( Xl )- Z (x i )=^Z-(K ) + 0(h') 
 
 p i,j (x o' ■ z ( x o' ■ °- 
 
 if u(x) denotes the approximative solution obtained with this 
 starting procedure and the first Nordsieck formula and 
 d(x) = u(x) - z(x), we have by the asymptotical theorem 
 (with q = 1 ) 
 
 d(x) = h 3 ( ei (x) + 0(h)) + h 3 (e I]; (x) + 0(h)), 
 
 -53- 
 
where 
 
 §- e_(x) - G(x,z(x))e T (x) + 5J z^\x), e_(b) = 0; 
 
 dx I ' I 24 
 
 ^ e n (x) = 0(x,«(x))e (x), a (b) - j| z (3) (b); 
 
 the term e JT (x) represent the initial error due to the starting procedure, 
 
 Errors for the Starting Procedures for the Modified Nordsieck Formula 
 
 All we have said above can be repeated word for word for the k modified 
 Nordsieck formula except that we must replace everywhere k by k+1. 
 
 h 
 
 Example 
 
 For the first modified Nordsieck formula, we have: 
 
 a 00 = a 01 = a 02 = ° ; 
 
 _5 _8 zl 
 a l0 "' 12' a il " 12' a i3 " 12' 
 
 14 1 
 
 a 20 3' a 21 "" 3' a 22 " y 
 
 4 5 
 
 R - 0, R x = 1^ z ( ^(b) + 0(h 5 ), R 2 = - ^ z (5) (b) + 0(h 6 ); 
 
 P ^.(x 2 ) - z(x 2 ) = - -^ h 5 G(b,z(b))z (1+) (b) + ^ z (5) (b) + 0(h 6 ); 
 
 p 6,j (x o ) - z,(x o } = °' p o,j (x i } - z ' (x i } = - G ( b > z ( b )) h z(i+)(b) + 0(h5) ' 
 
 if u(x) denotes the approximative solution and d(x) = u(x) - z(x): 
 
 ■54- 
 
d(x) = h 3 ( ei (x) + 0(h)) + h 5 (e i;[ (x) + e TI1 (x) + 0(h)), 
 
 where 
 
 |j e x (x) - G(x,z(x)) e] .(x) + g| z (J+) (x), ej (h) - 0; 
 ^ e IX (x) = G(x,z(x)) eiI (x), e n (b) = - j| G(b,z(b))z (i+) (b) + ^ z (5) (b); 
 ^ e Ilx (x) = G(x,z(x))e ni (x), e m ( b ) = 2B8 G(b, z(b) )z (U) (b ); 
 
 the starting error, represented by the term e + e is negligible in comparison 
 with the term e T « 
 
 •55- 
 
SECTION 8. CONCLUSION 
 
 The Nordsieck method, the modified Nordsieck method and all Nordsieck- 
 like formulas considered in Section k look very elegant. 
 
 Let us consider once more the matrix formulation of the multistep 
 method (3), (k), (5) and its Nordsieck-like equivalent (6), (7), (8). The 
 classical multistep method is such that A and a in (k) are the simplest possible, 
 i.e., each normal step involves the least possible arithmetics. If v is the 
 vector of the iterative process (7) and v* the vector obtained after the modifi- 
 cation of h, the Nordsieck-like method is such that the matrix D which transforms 
 
 — »■ — * 
 
 v in v* is diagonal; the corresponding matrix for the classical multistep 
 
 method has not this property. 
 
 In view of the fact that keeping the step constant is the rule and 
 modifying it is the exception, one may prefer the old and classical formulation 
 of the multistep method. 
 
 -56- 
 
SECTION 9- TABLE OF RESULTS 
 
 For the Nordsieck method and the modified Nordsieck method (k = 1,2, 3**0 
 we give successively: 
 
 1) The Nordsieck formulation of the method. 
 
 2) The corresponding multistep corrector formula: 
 
 y(x + h) = y(x) + h(p_ 1 f(x + h) + ... + 3 k f(x - kh)); 
 the truncation error is denoted by 
 
 T = Ch P+1 2 (p+l) , 
 
 meaning that if z(x) is a function which possesses (p+l) continuous 
 derivatives, then 
 
 z(x + h) = z(x) + h(p z'(x + h) + ... + az'(x - kh)) + T, 
 
 T = Ch P+1 z (p+l) (| ), x - kh < I < x + h. 
 
 3) The corresponding multistep predictor formula with the truncation 
 error. 
 
 k) The error involved in the modification of the length of the step. 
 Suppose the differential equation y' = f(x,y) is integrated by the 
 method we are considering with the step h* for x < b and the step 
 h for x > b. Let u(x) be the approximative solution, z(x) be the 
 exact solution such that u(b) = z(b) and d(x) = u(x) - z(x) for 
 x > b; the error satisfies the asymptotic relation 
 
 d(x) - h P (s(x) + 0(h)) + h q (t(x) + 0(h)); 
 
 the term s(x) represents the accumulation of the truncation errors; 
 the term t(x) is due to the modification of the step; in the defini- 
 tions of s(x) and t(x) we shall use the notations 
 
 -57- 
 
(x) . M^l/y = z(x ). (1) 
 
 h* 
 
 p = T' 
 
 5) The error involved in the starting procedure. The differential 
 
 equation y 1 = f(x,y) is integrated for x > b with prescribed initial 
 value at x = b; if u(x) is the approximative solution, z(x) the 
 exact solution such that u(b) = z(b), d(x) =f u(x) - z(x) ; then for 
 x > b 
 
 d(x) - h P (s(x) + 0(h)) + h r (w(x) + 0(h)); 
 
 the first term represents the accumulation of the truncation errors 
 and the second the error due to the starting procedure. The 
 notation (l) will be used. 
 
 Nordsieck Formula k = 1 
 
 y(x + h) = y(x) + h(f (x) + a(x) + -| (f(x + h) - f P )}, 
 
 f P = f(x) + 2a(x), 
 
 .(x + h) = a(x) + I (f(x + h) - f P . 
 
 Corrector Formula 
 
 y(x + h) = 3I (5f(x + h) + 8f(x) - f(x - h)); 
 
 1 ~2% z 
 
 -58- 
 
Predictor Formula 
 
 y(x +■*)• y(x) + - (3f(x) - f(x - h)); 
 
 1 i 2 z • 
 
 Modification of h 
 
 d(x) = h 3 (s(x) + 0(h)) + h 3 (t(x) + 0(h)); 
 s'(x) - g(x)s(x) + -^ J k) (x), s(b) = 0; 
 
 f(x) =» g(x)t(x), t(b) -i^£ z (3) (b). 
 
 Starting Procedure 
 
 d(x) - h 3 (s(x) + 0(h)) + h 3 (w(x) + 0(h)); 
 
 s'(x) = g(x)s(x) ± ^ z {k) (x), s(h) = 0; 
 
 w'(x) = g(x)w(x), v(b) =~ z (3) (b). 
 
 Modified Nordsieck Formula k = 1 
 
 y(x + h) = y(x) + h{f(x) + a(x) + b(x) + -| (f(x + h) - f P )}, 
 
 f P = f(x) +■ 2a(x) + 3b(x), 
 
 a(x + h) = a(x) + 3b (x) + ^ (f(x + h) - f P ), 
 
 b(x + h) = b(x) + | (f(x + h) - f P ). 
 -59- 
 
Corrector Formula 
 
 y(x + h) = j| (5f(x + h) + 8f(x) - f(x - h)); 
 
 — £- w - 
 
 Predictor Formula 
 
 y(x + h) = -| (23f(x) - l6f(x - h) + 5f(x - 2h)); 
 
 ,.#.»>. 
 
 Modification of h 
 
 d(x) = h 3 (s(x) + 0(h)) + h (t(x) + 0(h)); 
 
 (x) = g(x)s(x) + ^z ( ^(x), s(b) = 0; 
 
 t'd) - g(x)t(x), t(b) « -< 2p2 - T f + H .<«(»), 
 
 Starting Procedure 
 
 d(x) = h 3 (s(x) + 0(h)) + h 5 (w(x) + 0(h)); 
 
 (x) = g(x)s(x) +-i z (U) (x), s(b) = 0; 
 
 2K 
 
 w'(x) = g(x)w(x), w(b) = g g(b)z (4) (b) + ^ z (5) (b), 
 
 -60- 
 
Nordsieck Formula k = 2 
 
 y(x + h) = y(x) + h{f(x) + a(x) + b(x) + | (f(x + h) - f P )}, 
 
 f P = f(x) + 2a(x) + 3b(x), 
 a(x + h) = a(x) + 3b (x) + | (f(x + h) - f p ), 
 b(x + h) = b(x) + i (f(x + h) - f P )„ 
 
 Corrector Formula 
 
 h , 
 
 y(x + h) = y(x) + ^ (9f(x + h) + 19f(x) - 5f(x - h) + f (x - 2b)), 
 
 720 z 
 
 Predictor Formula 
 
 y(x + h) = y(x) + -| (23f(x) - l6f(x - h) + 5f( x - 2h)) ; 
 
 T = ^ z (^) 
 
 Modification of h 
 
 d(x) = h\s(x) + 0(h)) + h\t(x) + 0(h)),- 
 
 '(x) = g(x)s(x) + -|2 z (5) (x ), s(b) = . 
 
 t'(x) = g(x)t(x), t(b) =^^ z (4) (b). 
 
 -6i- 
 
Starting Procedure 
 
 d(x) = h*(s(x) + 0(h)) + h 5 (w(x) + 0(h)); 
 
 s'(x) = g(x)s(x) + ^|§ z (5) (x), s(b) = 0; 
 v'(x) - g(x)w(x), w(b) = - ^ g(h)z (4) (b) + ^ z (5) (b)„ 
 
 Modified Nordsieck Formula k = 2 
 
 y(x + h) = y(x) + h{f(x) + a(x) + b(x) + c(x) + | (f(x + h) - f p )}, 
 
 f P = f(x) + 2a(x) + 3b (x) + kc(x), 
 
 a(x + h) = a(x) + 3b (x) + 6c(x) + j| (f(x + h) - f p ), 
 
 b(x + h) = b(x) + i+c(x) + | (f(x + h) - f p ), 
 
 c(x + h) = c(x) + ^ (f( x + h) - f P ) 
 
 Corrector Formula 
 
 y(x + h) = y(x) + ^ (9f(x + h) + 19f(x) - 5f(x - h) + f(x - 2h)). 
 
 I9hf (5) 
 
 720 Z 
 
 Predictor Formula 
 
 y(x + h) = y( x ) + ^ (55f(x) - 59? (x - h) + 37f(x - 2h) - 9f(x - 3h)); 
 
 _ 2$lh^ (5) 
 1 '" 720 Z ' 
 
 -62- 
 
Modification of h 
 
 d(x) = h (s(x) + 0(h)) + h 5 (t(x) + 0(h)); 
 
 (x) = g(x)s(x) + -& z (5) (x), s(b) = 0; 
 
 f(x) = g(x)t(x), t(b) = - £- - f 2 + 1 z (5) (h), 
 
 Starting Procedure 
 
 d(x) = h (s(x) + 0(h)) + h 5 (w(x) + 0(h)); 
 
 (x) = g(x)s(x) + -g z (5) (x), s(b) = 0; 
 
 w'(x) = g(x)w(x), w(b) = ^ z (5) (b). 
 
 Nordsieck Formula k = 3 
 
 y(x + h) = y(x) + h(f(x) + a(x) + b(x) + c(x) + ||i (f(x + h) - f P )} 
 
 f P = f(x) + 2a(x) + 3b(x) + kc(x), 
 
 a(x + h) = a(x) + 3b(x) + 6c(x) + i| (f(x + h) - f P ), 
 
 b(x + h) = b(x) + kc(x) + i (f(x + h) - f P ), 
 
 (x + h) = c(x) + -i (f(x + h) - f P ) 
 
 2¥ 
 
 •63- 
 
Corrector Formula 
 
 y(x + h) = y(x) + =|^(251f(x + h) + 6k6f(x) - 26kf(x - h) + 106f (x - 2h) - 19f (x - 3h)); 
 
 1 ISo ' 
 
 Predictor Formula 
 
 y(x + h) = y(x) + g| (55f(x) - 59f(x - h) + 37f(x - 2h) - 9f(x - 3h)); 
 
 _ 2$lh^ (5) 
 1 720 z * 
 
 Modification of h 
 
 d(x) = h 5 (s(x) + 0(h)) + h 5 (t(x) + 0(h)); 
 
 s'(x) = g(x)s(x) + jgg z (6) (x), s(b) = 0; 
 
 f(x) = g(x)t(x), t(b) = - 10P j^g - 9 z (5) (b) 
 
 Starting Procedure 
 
 d(x) = h 5 (s(x) + 0(h)) + h 5 (w(x) + 0(h)); 
 
 s'(x) = g(x)s(x) 4- ^ z (6) (x), s(b) = 0; 
 V(x) = g(x)w(x), w(b) = g§ z (5) (b). 
 
 -6k- 
 
Modified Nordsieck Formula k = 3 
 
 251 
 
 y(x + h) = y(x) + h{f(x) + a(x) + b(x) + c(x) + d(x) + |g (f(x + h) - f P )}, 
 f P = f(x) + 2a(x) + 3b (x) + kc(x) + 5d(x) 
 
 ,(x + h) = a(x) + 3b (x) + 6c(x) + 10d(x) + || (f(x + h) - f P ), 
 
 b(x + h) = b(x) + kc(x) + 10d(x) + || (f(x + h) - f p ), 
 
 c(x + h) = c(x) + 5d(x) + 5! (f (x + h) - f P ), 
 
 d(x + h) = d(x) + ^ (f(x + h) - f p ) 
 
 Corrector Formula 
 
 y(x + h) = y(x) + =|rr (251f(x+h) + 6k6f(x) - 26kf(x- h) + 106f (x - 2h) - 19f(x-3h)). 
 
 720 
 
 T - - 4- z (6) 
 1 I£o z 
 
 Predictor Formula 
 
 y(x+h) = y(x) + T=!pr(l901f(x) - 277^f(x-h) + 26l6f (x - 2h) - 127^f (x - 3h) ) 
 
 720 
 
 + 251f (x - kh) } 
 
 T -22L (6) 
 1 "" 288 z 
 
 Modification of h 
 
 d(x) = h 5 (s(x) + 0(h))+ h 6 (t(x) + 0(h)); 
 
 -65- 
 
s'(x) = g(x)s(x) + rjk z (6) (x), s(b) = 0; 
 
 I5o 
 
 4.1/ > f ur \ 4-^^ JJ8P 11 - 7P 2 - 180P + 67 (6),,x 
 t'(x) = g(x)t(x), t(b) = ^ z v '(b), 
 
 Starting Procedure 
 
 d(x) = h 5 (s(x) + 0(h)) + h 7 (v(x) + 0(h)); 
 
 (x) = g(x)s(x) + J- z (6) (x), s(b) = 0; 
 
 160 
 
 w'(x) = g(xMx), w(b) = - j^_ g(b)z (6) (b) + J^ z (7) (b) 
 
 Nordsieck Formula k = h 
 
 y(x + h) = y(x) + h{f(x) + a(x) + b(x) + c(x) + d(x) + ^|| (f(x + h) - f P )}, 
 
 f P = f(x) + 2a(x) + 3b(x) + kc(x) + 5d(x), 
 
 a(x + h) = a(x) + 3b (x) + 6c (x) + 10d(x) + |? (f(x + h) - f P ), 
 
 b(x + h) = b(x) + 4c(x) + 10d(x) + || (f(x + h) - f P ), 
 
 c(x + h) = c(x) + 5d(x) + jj| (f(x + h) - f P ), 
 
 d(x + h) = d(x) + ^ (f(x + h) - f P ) 
 
 -66- 
 
Corrector Formula 
 
 y(x + h) = y(x) + j^ {^75f(x + h) + l427f(x) - 798f(x - h) + 482f(x - 2h ) 
 
 - 173f(x - 3h) + 27f(x - hh)}, 
 
 _ 863hl (7) 
 1 " ~ 60480 z • 
 
 Predictor Formula 
 
 y(x + h) = y(x) + ^|q (I901f(x) - 277^f(x - h) + 26l6f(x - 2h) 
 
 - 127^f(x - 3h) + 251f(x - kh)) 
 
 T _ 95h^ z (6) 
 1 288 z 
 
 Modification of h 
 
 d(x) = h 6 (s(x) + 0(h)) + h 6 (t(x) + 0(h)); 
 
 S ' (x) ^ (x)s(x)+ » z(T) ^ s(b) = 0; 
 f(x) = g(x)t(x), t(b) = - j^ (2i+0p^ - 35P 2 - 1792P + I587)z (6) (h), 
 
 Starting Procedure 
 
 d(x) = h 6 (s(x) + 0(h)) + h 7 (w(x) + 0(h)); 
 
 s'(x) = g(x)s(x) + ^g_ z ^)( x) , s(b ) = Oo 
 V(x) = g(x)w(x), v(b) = - J^ g(b)z (6) (b) + ^ z (7) (b)„ 
 
 -67- 
 
Modified Nordsieck Formula k - h 
 y(x + h) - y(x) + h{f (x) + a(x) + b(x) + c(x) + d(x) + e(x) + ^|| (f(x + h) - f p )}, 
 f P = f(x) + 2a(x) + 3b(x) + 4c(x) + 5d(x) + 6e(x), 
 a(x + h) = a(x) + 3b (x) + 6c(x) + 10d(x) + 15e(x) + ^ (f(x + h) - f p ), 
 b(x + h) = b(x) + 4c(x) + 10d(x) + 20e(x) + | (f(x + h) - f P ), 
 c(x + h) = c(x) + 5d(x) + 15e(x) + U (f(x + h) - f P ), 
 d(x + h) = d(x) + 6e(x) + ^ (f (x + h) - f p ), 
 e(x + h) = e(x) + -|- (f(x + h) - f P )„ 
 
 Corrector Formula 
 
 y(x + h) = y(x) + ^^ (U75f(x + h) + lk-27t(x) - 798f(x - h) + 482f(x - 2h) 
 
 173f(x - 3h) + 27f(x - kh)) 
 
 _ 863^1 Z (T) 
 1 " " 60480 z ' 
 
 Predictor Formula 
 
 y(x + h) = y(x) + ^q (U227f(x) - 7673f(x - h) + 9^82f(x - 2h) 
 - 6798f(x - 3h) + 2627f(x - kh) - 425f(x - 5h)). 
 
 _ 19087h 7 (7) 
 1 " 60480 z 
 
 -68- 
 
Modification of h 
 
 d(x) = h 6 (s(x) + 0(h)) + h 7 (t(x) + 0(h)); 
 s'(x) . g(x)s(x) t ^g z (7) (x), s(b) = 0; 
 
 t'(x) = g(x)t(x), t(b) = - —jL (i6p 5 - 3P 3 - Mp + 28)z (7) (b) 
 
 Starting Procedure 
 
 d(x) - h (s(x) + 0(h)) + h 7 (w(x) + 0(h)); 
 
 •7 
 
 s'(x) = g(x)s(x) + fg^ z (7) (x), s(h) = 0; 
 
 w'(x) = g(x)w(x), w(b) =I |g_ z ^) (b)< 
 
 -6 9 - 
 
APPENDIX 1 
 
 Let x be a vector of m components x, , x^, .... x . f(x) be a vector 
 
 e 1' 2' ' m 
 
 function of x with components f (x), ..., f (x). For any vector z f 
 
 m 
 
 Z I z . I is the norm of z . 
 1-1 X 
 
 Theorem 
 
 If for any two x and x*, there exists a constant c < 1 such that 
 
 ||f(x*) - f(x)| | < c||x* - x||, (1) 
 
 then 
 
 a) the equation 
 
 x = f(x) (2) 
 
 — * 
 possess one and only one solution y; 
 
 b) for any x n , the sequence x_, x , ..., x , ... defined by 
 
 x , = f(x ) (3) 
 
 n+1 n' v ' 
 
 converges towards y. 
 
 Proof 
 
 This theorem is only a particular formulation of classical results 
 in Banach spaces theory. The proof results from the two statements: 
 
 a) the equation (2) has at most one solution; 
 
 b) the sequence defined by (3) converges towards a solution of (2). 
 
 a) Let y and y* be two solutions of (2). By (l): 
 
 My* - y|| ■ l|f(y*) - ?(y)| | < c||y* - y| I, 
 
 -70- 
 
and consequently | | y* - y| | - 0. 
 
 . -» 
 
 b) Let x -, , x ...... x be the components of x . In order to 
 
 ' n, 1' n,2' n,m n 
 
 prove the convergence, it suffices to show that for every i, 
 1 s 1. . .., m, the sequence x . is a Cauchy sequence. Let 
 p > q > N: 
 
 i i I i -* ~ * i i l I ~* ~~ * II I I - * ~ * I I I I ~* — * l l 
 
 |X p,i - X q,i' " l|X p- X q M - l|X q- X p-l M + ' 'Vl " X p-2 ' ' + °°° + "Vl"^ 1 ' 
 
 (k) 
 
 by (1): 
 
 ||Xg - x 1 || = ||f(x 1 ) - f(x )|| < cl|x 1 - x ||, 
 
 ||x 3 - x 2 || = ||f(x 2 ) - Z&JW < c||x 2 - xjl < c 2 !!^ - x ||, 
 
 |x . - x I I < c j Jx n - x„J I; 
 
 1 n+1 n n - ' ' 1 ' " 
 
 replacing in (h). we get 
 
 p-1 oo n 
 
 n c M-* -> 
 
 x . - x . < x - x Z c < x_ - xJ Z c = - — 
 p, i q, l — 1 O 11 — ' ' 1 O' 1 „ 1-c'l 
 *> *' n=q n=N 
 
 since c < 1 and | |x - x | | is constant, for given e, it is possible to 
 
 N such that x . - x 
 p,i c 
 
 Let y = lira x : by (3) 
 
 find N such that llx . -x .11 < €, p,q > N. 
 
 p,i q,i' - 
 
 n 
 
 y = lim x . = lim f (x ), 
 J n+1 v n'- 
 
 n-x» n-*» 
 
 by (l), f(x) is continuous and consequently lim f(x ) = y, q„eod, 
 
 n-*» 
 
 -71- 
 
APPENDIX 2. SOME EQUIVALENT "METHODS' 
 
 Let y' = f(x,y) be a differential equation to be integrated at 
 x„. x, ...,x, ...,x = x_ + nh. We consider h as a fixed number and denote 
 
 '0> ' V 
 
 n 
 
 by y the approximative solution at x . We define two processes for computing 
 
 Method 1 
 
 For n > 0: 
 
 -vp-> 
 
 y , = p v + yf . , 
 J n+1 * n n+1' 
 
 v ., = Av +f _, a , 
 n+1 n n+1 ' 
 
 f . = f (x _,y , ); 
 
 n+1 n+1' n+1 ' 
 
 (1) 
 (2) 
 (3) 
 
 A is a matrix of order m with constant elements; p and a are constant vectors; 
 
 — > — * 
 
 y is a scalar constant; v„ is the initial vector,, For given v_, this method 
 
 allows to compute successively y , y , 
 
 ., provided that the non-linear system 
 of equations (l), (3)> has always a unique solution (we shall always suppose 
 that it is the case). 
 
 Method 2 
 
 For n > 0: 
 
 y . = s , + t , . 
 
 J n+1 n+1 n+1' 
 
 GO 
 
 S n+1 = Vn + Vn-1 + •" + a jVj + Vn + Vn-1 + "" + ¥n-k + 6f n+l 
 
 (5) 
 t = for all n, except for n = w, w+1, ...., w+p, where w = max(0,k~j); 
 
 -72- 
 
j, k, p, Q!_, ...,0L., 3 n , . ..,3 , 5 are the constants which characterize the 
 
 method 2. Suppose that j > -1, k > -1, p > -1; if j = -1, there are no a-terms 
 
 in (5); if k = -1, there are no (3- terms in (5); if p - -1, all t-terms are zero 
 
 and the method (2) becomes the classical multistep method; furthermore a./ 
 
 if j > 0, B k / if k > 0; t v+v t w+2 , ...,, t w+p+] _, fl Q| s_ i; ..., s_., f Q , f_ ± , 
 
 f , .... f , are the initial values: let u_ denote the vector of dimension 
 -2 } -k '0 
 
 p+j+k+3, the components of which are the initial values. 
 
 Equivalence 
 
 A particular method 1 and a particular method 2 are said equivalent if 
 there exist two rectangular matrices P and Q with constant coefficients (i„e. ; 
 independent of f(x, yj, u_, v ) such that if v = Pu for any given u or 
 u = Qv n for any given v., then methods 1 and 2 give the same values for y } 
 n = 1,2,.-.. 
 
 Theorem 1 
 
 a) For any method 1, there exists an equivalent method 2. 
 
 To ) For any method 2, there exists an equivalent method 1. 
 
 c) Two methods of type 2 are equivalent if and only if all their 
 characteristic parameters are identical. 
 
 Theorem 1 has the following practical meaning: all methods 1 are 
 "essentially equivalent" to a classical multistep method. These methods 1, 
 which are not rigorously equivalent to a classical multistep method (i.e., which 
 are equivalent to a method 2 with p > 0), do not seem to have much interest. 
 Therefore all "interesting methods l" can be discussed in terms of multistep 
 formulas . 
 
 The proof of theorem 1 which is rather painful and long, will not 
 be given here. 
 
 We now consider two other methods, method III (multistep method) and 
 method IV. 
 
 ■73- 
 
Method III 
 
 For n > 0: 
 
 y , ■ a.y + a. y ,+...+ a . y .+ B n f + 0, f . + ..'.+ . 'f . + 7 f 
 J n+1 CTn l J n-l J n-j K n K l n-1 j+p n-j-p 
 
 n+l J 
 
 f n+1 " f ( x n+ l^n + l ) ' 
 
 i. D, a , .... a , B~» .... B. 3 7 are the constants which characterize the 
 
 o} r> q> > j^ t-Q, , J+P 
 
 method; j > 0, a / 0; p > 0; if p > 0, then B. +p / 0; y Q , y_ ± , .,., y n _. 
 •J-P 
 
 f^, f ,, .... f .' are the initial values 
 0' -1' ' -J-T 
 
 Method IV 
 
 For n > 0: 
 
 Z n + 1 = Vo,n + Vl,n + "' + Vj,n " 6 j + l g n ' & j+2 S n-l "" "' " 6 j+P g n-p+l + 7 W 
 
 g n+1 = f ^ n+ l^n + l ) 
 
 z^ . = <x.z_ + a n z n + ... + a.z. + (8. + y)g , 
 0,n+l 0,n 1 l,n J J,n v j /to n+l^ 
 
 z. . = z. . + 5. .g ., i = 1,2, . „ „,j; 
 
 i,n+l i-I, n j-i°n+l' ' ' *°' 
 
 j> P» Qvn, • ••> Ct., 5 , ..., 6., 5. , . 6. , 7 are the constraints which 
 o> *>> q> > j> o' J J+l J+P 
 
 characterize the method; j > 0, a. £ 0, p > 0; if p = 0, there are no 6- terms 
 
 J 
 
 in the expression of z^ ± ; if p > 0, 5 J+p / 0; z Q Q , . ..., z^. Q , g Q , g_^ . 
 
 1 ^ 
 
 g are the initial values. 
 
 "*- IT 
 
 Theorem 2 
 
 If the coefficients B_, ..., B. of method III and the coefficients 
 
 0' J+P 
 
 S~f &.,..., 5. of method IV satisfy the relations 
 0* 1' ' J+p J 
 
 -Ik- 
 
J 
 
 P. = Zq 6. . - 8. ... i = 0,l,...,p+j; (6) 
 
 K i _ r j+i-r j+l+i 
 
 where one defines 6. . = for i > 0, and if the relations 
 
 j+p+i 
 
 S n-i = f n-i' 1 = OAj.-.iP-I^ (7). 
 
 p+i 
 
 z. = y . + E 6. . f , i= 0,1, ...,j, (8) 
 
 i,n J n-i _ J-i+r n-r^ ' ' '° 3 v 
 
 3 r=0 
 
 are satisfied for n = 0, then they also hold for n = 1. 2. . . . and z = y „ 
 
 3 J 3 n n 
 
 The proof of theorem 2 by induction with respect to n is straightfor- 
 ward and will not be given here. 
 
 Remarks 
 
 1) For given p., (6) is a linear system of equations for 8„, . . . , 8 . ; as it is 
 
 easy to show the matrix of the system is triangular; all diagonal elements 
 
 are equal to a.; by hypothesis a. f and the matrix is regular. This 
 
 J J 
 
 establishes the equivalence of methods III and IV. 
 
 2) One iteration in methods III or IV involves practically the same number of 
 
 operations. However method III needs 2j + p + 2 "past values" y 3 y ,..., 
 
 y ., f ,...,f . (supposing p > 0), but method IV needs only j + p + 1 
 
 "past values," z^ , z., ..... g ,...,g , . That curious property would 
 ' , r 0, n' ' IfXx' n' ' to n+l-p * * J 
 
 have perhaps been useful in a time where the computers had only small 
 memories. One can prove that for two different sets of initial values in 
 method IV, the corresponding sequences of values z are also different; 
 
 this is not the case for method III, since the equations (7)> (8) considered 
 
 as a linear system for y,...,y ., f,...,f . have more unknowns than 
 
 n' 3d n.-y n 3 ' n-j-p 
 
 equations; however, the corresponding matrix has the maximal rank, 
 
 -75- 
 
REFERENCES 
 
 1. Arnold Nordsieck. On numerical integration of ordinary differential 
 
 equations. Mathematics of Computation. January, 1962. 
 
 2. Peter Henrici. Discrete variable methods in ordinary differental equations. 
 
 John Wiley and Sons, Inc., New York, London. 1962. 
 
 3. W. E. Milne. Numerical solution of differential equations. John Wiley and 
 
 Sons, Inc., New York. 1953- 
 
 -76-