TO Xk if ISM ■ 1 SHI hHSnII H Bsfli 111 ■ HI M nmmfflh 8SIW ■i Ml LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 510.84 IJKor "0.698-702. £ 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. UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN RE L161 — O-1096 J_#(tAS UIUCDCS-R-75-702 COO-2383-0017 ANALYSIS OF FIXED-STEPSIZE METHODS • by Robert Skeel February 1975 THI LlbRARY OF MAY 7 19/5 DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN URBANA, ILLINOIS Digitized by the Internet Archive in 2013 http://archive.org/details/analysisoffixeds702skee UIUCDCS-R-75-702 COO-2383-0017 ANALYSIS OF FIXED-STEPSIZE METHODS* by Robert Skeel February 1975 DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN URBANA, ILLINOIS 6l801 * This work was supported in part by the Atomic Energy Commission under Grant US AEC AT (ll-l) 2383. PREPRINT: Submitted to SIAM Journal on Numerical Analysis. ABSTRACT The unified theory developed by Henrici [3] for one-step methods is extended to the more general case where the order of the system of difference equations can exceed the order of the system of ordinary differential equa- tions. The analysis is applicable to every fixed-stepsize fixed-formula method known to the author. For many of these methods the concept of consis- tency is inadequate. A more appropriate concept, termed quasi-consistency, is introduced. 1. Introduction . Most numerical methods for solving initial value problems in ordinary differential equations discretize the system of differential equations and then solve the resulting system of difference equations. Henrici [3] has developed the theory for one-step methods. This paper studies the more general case where the order of the system of difference equations may exceed the order of the system of differential equations. Just as ordinary differential equations are more easily studied when written as a first order system, it is likewise profitable to study difference equations as a first order system: u = S u + h ±(t , , u ; h). -n -n-1 ■*■ n-1 -n-1 This one-step formulation is quite powerful, for it renders the analysis more transparent . Two basic properties of methods are of practical concern: convergence, which means that the accumulated discretization error goes to zero as the step- size h ■*■ 0, and stability, which ensures that the accumulated roundoff error does not grow "too fast" as h •+ 0. Some authors (for example, Chartres and Stepleman [l]) combine the concepts of convergence and stability. This approach is cumbersome, and it does not reflect current practice in the sense that generally the roundoff error does not tend to zero like some positive power of h. Convergent one-value (one-step) methods generate difference equations which are discrete analogs of the differential equation, and as a result the stability properties of the differential equation are preserved by the difference equation for small enough h. For multivalue methods the difference equations are not analogs of the differential equations and so the stability properties of the differential equation might not be faithfully represented. In order to ensure that this is not the case, we ought to require that a multivalue method be "as stable as" a one-value method. The notion of being "as stable as" is defined in §3 in a quite natural way. And then in §U it is demonstrated that being as stable as a one-value method is equivalent to requiring that the matrix S have all of its eigenvalues inside the unit circle except for an eigenvalue at 1 whose index is one and whose multiplicity is equal to the number of initial values required for the differential equation. There are common methods, like Milne's method, which are not as stable as a one-step method but which can be expressed in a form for which they satisfy the strict root condition: all the roots of the minimal polynomial of S are inside the unit circle except possibly for a root at 1 . It is determined that for methods satisfying the strict root condition convergence is equivalent to a property termed quasi-consistency. If d , n = 0(l)N, denotes the local dis- cretization error, then quasi-consistency means that max i , i _ , _ d -* as h ■*■ £ n £ N ' — n ' and max j E ( +cL+"'+d)|-»-0 as h + $ n £ N ' -0 -^ -n ' where E = S . On the other hand consistency usually means that l<^l + |d 1 ' + , *' + |cLl~ >0 as h + (see Chartres and Stepleman [l], Stetter [7» pp. 5 5 75]). It is shown that the order of convergence is always equal to the order of quasi-consistency. For an often-studied class of methods which Stetter [7, p. 310] calls straight multi step methods , quasi-consistency is equivalent to consistency. This is not true for more general methods. For example, Nordsieck'i [5] method is quasi-consistent, and hence convergent, of order 6 and yet only consistent of order 5- Other examples are given in §5. To have convergence of order p, it suffices to have that d = 0(h ) and that the essential local — n truncation error E d = 0(h ). In the case of linear multistep methods E d is just the local truncation error divided by 8„ + 8 + •" + B (cf. Gear [2, U 1 q p. 118]). In §6 the leading term in the asymptotic expansion of the global error of a method satisfying the strict root condition is determined in terms of the leading term(s) in the asymptotic expansion of the local discretization error. 2. Fixed-stepsize methods . Let s be a fixed positive integer. We consider a fixed class of systems of s first order differential equations (2.1a) y' = f(t, y) where f is continuous in t and uniformly Lipschitz continuous in y for s S t < 1, y€=R. For any initial condition (2.1b) y(0) = y Q a unique solution y(t) is guaranteed. As an example, the class of systems under consideration might consist of those equations having the form (y 1 )' = y 2 , (y 2 )' = f(t, y\ y 2 ) . Let q and k be fixed positive integers with k >, s. The methods under consideration have three components: (i) a mapping which to each f assigns a starting procedure o(y; h)^-R k continuous for y€^ R S and £ h £ h where < h <, 1/q. ( ii ) a mapping, called a formula, which to each f assigns a forward step procedure S u + h _^(t, u; h) G R k where S is a k x k matrix independent of f and the increment function _£ is uniformly Lipschitz continuous in u and continuous in t and h for 0$ t Jl, ugH , and S h S h . (iii) a mapping which to each solution y(t) of a problem assigns a correct value function v_(t; h) = A y(t) + h z.(t; h) £ R k where A is a k x s matrix of rank s independent of y(t) and _z(t; h) is bounded for $ (q-l)h S t S 1. Stetter [7, p. 7] introduces a similar mapping, which he denotes by A . n Note that y(t) can be recovered "by means of y(t) = (A T A) -1 A T v_(t; 0). Given a problem and an h ^ h , a method generates a uniform grid and a discrete problem on that grid. The grid is (q - l)h = t Q < t 1 < ... < t H S 1 where t = t + nh and N is the largest integer such that t £ 1. Since there is no danger of confusion we do not exhibit the dependence on h of such quantities as t and N. The discrete problem on the grid is (2.2a) Uq = a(y Q ; h), (2.2b) u = S u . + h j)(t . , u . ; h) s n = l(l)N, — n — n-1 ■*■ n-1 — n-1 which can be solved to yield a numerical solution U= (uq, u 1} ..., %) T . The discrete problem (2.2) is normalized so that the unknown data vector u — n appears on the left hand side. This normalization serves to define o_, S, and i£_ uniquely. The correct value function y_(t; h) assigns a meaning to the numerical solution U by means of 1= (zo> z v •••> i^) T where y = y_(t ; h) . Hence the global discretization error is defined to be U - Y . Remark. An unfortunate aspect of methods with k > s is that the data within u are normally inconsistent with the differential equation. For T example, if u = [u , u ] , then there is usually no solution y(t) of the differential equation such that y(t )= u and y(t ., ) = u n n n-1 n-1 Example I. A straight q-step method (Stetter's [7, p. 310] terminology) determines one new value for each step: u = S u + h e n i|i(t _ , u , ; h) -n -n-1 —1 n-1 —n-1 where S = u = — n [ V Vi' • • > u n -q+l J ' " a i -a 2 ... -a q- i q 1 l k and e ± = [1, 0, ..., 0] T . T Also v_(t; h) = [y(t), y(t-h), ..., y(t - (q-l) h) ] . For the special case of a linear q-step method i|>(t, u; h) = g Q 4> + i f(t-j h, u J ) where solves the equation 4> = f(t, h U + { -I n (- a . u J + h 6. f(t-j h, u j ))}). Example 2. Perhaps the most important class of methods which are not straight multistep methods are P(EC) multistep methods. For example a PEC method with a second order Adams-Bashforth predictor and a third order Adams- Moulton corrector can be written u u* n-1 10 10 1 u n-1 n-1 u n-2 + h 5/12 2/3 -1/12 3/2 -1/2 ■n-1 "n-2 where ,P n-2 f- „ = f(t_ oJ u? J, n-2 9 n-2' f 9 . = f(t . , u P , ), n-1 n-1 n-1 >P = f(t , u , + |h f P , - i-h fP _). n n' n-1 2 n-1 2 n-2 T Naturally y_(t; h) = [y(t), y(t), y(t-h)] . Alternatively it could be written u n 1 2/3 -1/12 n-1 h u' n = h u 1 . n-1 h u 1 . n-1 1 h u' n-2 5/12 • h f(t , u , + | h u' - \h u' ) n n-1 2 n-1 2 n-2 with y_(t; h) = [y(t), h y'(t), h y'(t-h)] T . Example 3. Butcher's three-stage formula of order four (cf. Stetter [7, p. 275]) is given by M — *" u n 1 U n-1 + h h U n-l/2 h U n-3/2 1/6 2/3 1/6 10 "n-1 n-1/2 f 5 n where r n-l " f(t n-l- "„-!>' n-1/2 = f(t u + fh f =. n t n-1/2' n-1 k n-1 " h n-3/2 ), = f(t , u + 2h f . ._ - 2h f _ + h u' _ /0 ] n n-1 n-1/2 n-1 n-3/2 and by y_(t; h) = [y(t), h y'(t-h/2)] . Example 4. A modified linear multistep method equivalent to a fourth order Adams -Moult on method (Gear [2, p. 150] with M = <=°) is given by u n 1/2 1 1/2 1/8 u . n-1 3/8 h u' n h u' , n-1 1 n-1 1/2 1/3 1/2 5/2U Va + h -1/2U h u' , n-1 1 hU n-2 ij^(t , , u n ; h) n-1 -n-1 where y_(t;h) = [y(t), h y*(t), y(t-h), h y'(t-h)] T and i|;(t, u.i h) = «p the solution of ip = f(t, \ u 1 + u 2 + | u 3 + | u U + | h ^) i i k 3. Convergence and stability . Let | * | be a norm for R , and define I I tt v I I max I I II H-I IL- « n * N I Hfc-Zn I- DEFINITION 3.1. A method is said to be convergent for the problem (2.1) if II U-Y_ll O o" ) ' 0ash " > ' - If the method is convergent for all infinitely smooth problems (2.1), then it is simply said to be convergent. In practice U cannot be computed exactly; we must take into account the local computing error which arises from doing operations in finite precision. Hence we consider the perturbed numerical solution LL = (Hq> Has •••> Ujj) defined by Uq. = ^0 ; h) + £o' u — n =Su . + h *(t , , u _ ; h) + r , n = l(l)N, -n-1 *- n-1 -n-1 — n m a for any perturbation R_ = (r^, r_ , . . ., r ) . The dependence of U on R is not made explicit in our notation. Various definitions of stability as an absolute concept have been proposed but most of them seem somewhat arbitrary. (An exception is the use of Spijker's [6] norm to define stability for straight multistep methods.) It is more natural to define stability as a relative concept. DEFINITION 3.2. Method 2 is said to be at least as stable as method 1 if for each problem there is a fixed number K such that II " (2) -u (2) il * k || i (1) -u (1) il for any h $ h and any perturbation R ^ IT + . If , in addition, method 2 is at least as stable as method 1, then the methods are said to be equally stable ; otherwise method 1 is said to be move stable than method 2. 10 Since good stability properties are important for practical algorithms, it is desirable that a method satisfy the most stringent stability requirement that is compatible with convergence. DEFINITION 3.3. A convergent method is optimally stable if there is no convergent method which is more stable. Note. In general it is not possible to compare the stability of two arbitrary methods, and hence two optimally stable methods are not necessarily equally stable. Remark. Spijker's [6] definition of optimal stability differs from ours in that it requires an optimally stable method to be at least as stable as every other method in the class of methods under consideration. It is shown in §U that the optimally stable methods are those convergent methods that satisfy the very strict root condition, which is defined below. DEFINITION 3.^. A method is said to satisfy the root condition (RC) if the roots of the minimal polynomial of S are either inside the unit circle or on the unit circle and simple. A method satisfies the strict root condition (SRC) if the roots of the minimal polynomial are inside the unit circle except possibly for a simple root at 1. A method satisfies the very strict root condition (VSRC) if it satisfies the SRC and if S has at most s eigenvalues equal to 1. It is shown by Theorem 3.9 that the matrix S must have at least s eigenvalues equal to 1 in order for the method to be convergent; however any additional eigenvalues on the unit circle impart a degradation to the stability of the method. There is no increase in the complexity of the theory if instead of 11 assuming the VSRC we merely require that the strict root condition be satisfied. Furthermore, it is then possible to treat marginally useful methods like Milne's method. Hence the theory is developed for methods satisfying the SRC. The following lemma gives the best possible qualitative bound on the difference U - U in terms of R for methods satisfying the RC. First we introduce some notation. For any k x k matrix T let [T] denote the (N+l) x (N+l)k matrix I -T I -T I Then the quantity max J - M oo :? n $ N ' j=0 n n-1 L, T r ■J (N+l)k is a norm for R LEMMA 3.5- Given a method satisfying the RC and a problem (2.1), there exist positive constants c and C such that for any h $ h„ and any R €R S M U - U I I * C O I I — I I 00 ,(N+l)k c II [S]" 1 R [S] _1 R Proof. Set 6 = u - u . Then -n -n -n 6 =SS +h6 +r — n —n-1 —n-1 — n (3.1) where 5 n = lk.(" t n> u -,;h)-\J;(t ,,u ,;h) —n-1 *- n-1 -n-1 *- n-1 -n-1 By assumption there exists a constant L such that 4-i I s L 4-i 12 Solving the difference equation (3.1) gives (3.2) n n-1 ~ n n-1 By the root condition there exists a constant B such that i «n (3.3) Thus (3.2) "becomes L <: B. n-1 6 $ h B .E n 6, -n ' j=o ' -j [S]" 1 R By induction on n it follows that 6 — n -1 $ (1 + h B) n || [S]" 1 R || .< e* || [S] 1 R | |, which proves the second inequality in the lemma. From (3.2) and (3.3) ,n-J max i j£o S— r. | S (1 + n h B) Q ^ „ | ^ |; « (1 + B) | | U - U ||„, and thus I | S _1 R j ! $ (l + B) | I U - U j . Q.E.D. II — I ico 'II — — I loo When a norm, such as- | I [S] • M , permits an upper hound on U - U, then the method is said to he stable with respect to that norm, and if the norm also permits a lower bound on U - U, then the norm is said to he a minimal stability functional for the method. These ideas are discussed in more detail by Spijker [6] and Stetter [7, pp. 81-8U] . Of considerable importance to the theory is the component of S (cf. Lancaster [U, p. 162]) corresponding to the eigenvalue 1, which we denote by E. Assuming S satisfies the RC, we have that E = p"(l) _1 p"(S) where p(£) = (g - l)p(^) is the minimal polynomial of S. If S does not have an eigenvalue at 1, then we define E = 0. It is easily verified that E is idempotent and that S E = E = E S whence S n = E + (S - E) n , n = 1, 2, ... . If S satisfies the SRC, the eigenvalues of S - E are inside the unit circle. 13 The following theorem is for methods satisfying the SRC, and it gives bounds on the difference U - U in terms of the more useful norm I I r-,-1-1 „ I I _ max I n ~ 1 I LE J R = „ . T E .Z. r. + r l|LJ -i'oo s n £ N ' j=0 -j -n ' THEOREM 3.6. Given a method satisfying the SRC and a problem (2.1), there exist positive constants c' and C such that for any h £ h and any EeR (H + l)k c 1 [E] _1 R | ! £ | U - U £ C [E] -1 R Proof. The theorem follows from Lemma 3.5 if it can be shown that [El R I can be bounded by some constant times ' [S] R and vice versa. First of all II [E]" 1 R II.- || [E] -1 [S] [S]" 1 R || = || [S - E] [S] _1 R II. = S T SH l Jo 4 - S z^ * fa lit,,.!, j^; h) - ^ , d . 1(1)1. DEFINITION 3.7. A method satisfying the RC is said to be quasi- consistent with the problem (2.1) if max 1 , 1 ~ , _ ..id M- as h -+ ■ $ n <: N ' -n ' and max I ^ +d+...+d)|+0 as h ■* 0. < n $ N ' —0 —1 -n .' A method is consistent with (2.1) if |d^|+|d| + ...+ |d_|-»-0 as h -► . Clearly consistency implies quasi-consistency. For a method to be quasi-consistent with a problem it is sufficient that (i) d = o(l) , n = 0(1)N, — n (ii) E d = o(h) , n = l(l)N where we employ the convention that whenever the little-oh or big-oh notation is used, the implied limit or bound is uniform w.r.t.n. The quantity Ed is the essential local discretization error. — n THEOREM 3.8. A method satisfying the SRC is convergent for (2.1) if and only if it is quasi- consistent with (2.1). Proof. When R_ = - D_, then U = _Y ; i.e., the true solution is a perturbed numerical solution. (This unified treatment of discretization error and roundoff error is due to Chartres and Stepleman [l].) Thus from Theorem 3.6 convergence is equivalent to „ ■■ E(cU +d+...+d , ) + d +0 as h + 0, $ n $ N ' '—0 —1 -n-1 -n ' ' which is equivalent to quasi-consistency. Q.E.D. 15 We extend a theorem of Henrici [3, p. 71] stating necessary and sufficient conditions for convergence. THEOREM 3.9. A method satisfying the SRC is quasi-consistent with the differential equation (2.1a) if and only if (i) o_(y; 0) = Ay, (ii) S A = A, (iii) E ±{t t Ay; 0) = E A f(t, y) . Proof. Using the fact that in the region $ t S 1, $ h <: h continuity is equivalent to uniform continuity, we have (3.10 d^ = a(7 ; 0) - A 7 Q + o(l), d =SAy -Ay + 0(h) -n n n (3.5) = (S A - A)y n + 0(h), n = l(l)N, and E ^ = E {v^ + h i(t n _ l5 v^; h) - yO = E { ^-l " ^n + h ^ (t n-l' A y n-l ; ° )} + ° (h) ' n = l(l)N ' Therefore } h « 4, = 1 % - ^ ♦ jlo h 1ft,. A yj; 0)} + 0(1) t = E {A(y - y ) + / n j,( T , A y( x ) ; 0) dx} + o(l) U n Z t (3.6) = E Q / n {-A f(x, y(x)) + i(x, Ay(x); 0)}dx + o(l). Assume (i), (ii), and (iii) hold. Then it clearly follows that ^0 = o(l) , d^ = 0(h) , n = 1(1)N n ,2. E d. = o(l) , n = l(l)N, 0=1 —J which are sufficient for quasi-consistency. Now assume the method is 16 quasi-consistent with (2.1a). Then d = o(l), n = 0(l)N, and so from (3.*0 and (3.5) we have a(y •, 0) - A y Q = and (S A - A)y = 0. Parts (i) and (ii) immediately follow because y is arbitrary and y(t) is continuous. By convention equation (3-6) means that there is some function n(h) tending to as h ->• such that t 1 /T* N I E X n {_A f(T ' y(T)) + i(T ' A y(Th ° )}dT I * n(h) * Hence for fixed £ t S 1 | E _/* {-A f(x, y(x)) + £(t, A y(t); 0)}dx | < n( TIT") » (J n - q + 1 n = q, q + 1, •••? and so the left-hand side must be zero. Since t is arbitrary and the integrand is continuous, E {-A f(t, y(t)) + ^(t, A y(t); 0)} = By varying y , y(t) can be made to assume any value and so (iii) follows. Q.E.D. For a convergent method satisfying the VSRC part (ii) of Theorem 3.9 implies the existence of a unique k x s matrix M such that T E = A M . The rows of M are linearly independent left eigenvalues of S corresponding to the eigenvalue 1. Part (iii) of Theorem 3.9 thus reduces to (3.7) M T Mt, A y; 0) = f(t, y). Remark. Conditions (i)-(iii) of Theorem 3.9 are equivalent to convergence for methods which merely satisfy the RC if we strengthen the smoothness conditions on z_(t; h) as follows: (i) | _z(t; h) - z_(t; 0) | <: n(h) where n(h) ■> as h ■* 0, (ii) z_(t; 0) is Riemann integrable (cf. Chartres and Stepleman [l, p. U8U]). We can use the fact that I -1 (S- E)J| IT is bounded independently of n to show that for each problem (2.1) there exist positive constants c" and C" such that for any h S h and R_^ R c" | | [E]" 1 R I 1^ < | | U - U | 1^ $ C" | | R | | where I I O I I maX I r T3 I I I v I I I I £ I I = S n < N I jio E ^ I + I ^0 I + nil ! ^n " ^n-1 ! ' Just as before conditions (i) - (iii) of Theorem 3.9 are necessary for convergence. On the other hand, the additional restrictions on z(t; h) make it possible to show that N I | d - d | = o(l), n=l ' — n — n-1 ' and hence that conditions (i) - (iii) are sufficient for convergence. Example I. Let us examine the quasi -consistency conditions for a straight multistep method satisfying the SRC. We have A = e = [1, 1, ..., 1] T and so part (ii) of Theorem 3.9 requires that (3.8) -a -a - ... -a =1. 12 q Since the minimal polynomial of S is £ + ot_ £ + • • • + a , the SRC is 1 q T T equivalent to the VSRC, and hence E = _e m where m is the left eigenvalue of T S normalized so that m e_ = 1. It is easily verified that m =p(l) [l, 1 + a-.,..., 1 + a. + . . . + a , ] — 11 q-1 where p"(0 = £ q + (1 + aj £ q_1 + ■•• + (1 + o. + ••• + a .). 1 1 q-1 Because i^(t, u; h) = e_ ^(t, u; h), we have from (3.7) that (3.9) i|>(t, A y; 0) = p"(l) f(t, y). Together with part (i) of Theorem 3.9» (3.8) and (3-9) constitute the quasi- consistency conditions for straight multistep methods. In this case 18 quasi-consistency is equivalent to consistency "because the local discretization error e., d and the essential local discretization error p(l) e d are both —In — n multiples of the same quantity d . Example 2. The SRC does not exclude those weakly stable linear multistep methods based on numerical quadrature. For example, two steps of Milne's method, u =u + §• ( \t + -t , 4 f o)> n n-2 2 3 n 3 n-1 3 n-2 with stepsize h/2 are equivalent to one step of the method 1 u n-1/ 2 1 n-1 u n-3/2 + h l f 4f +if 6 n 3 n-1/2 6 n-1 if 4-i^ . 1 - 6 n-1/2 3 n-1 6 n-3/2 which does satisfy the SRC though not the VSRC. The above pair of difference equations is the discretization of a system of two identical differential equations, which is not surprising because the derivation of Milne's method involves two numerical integrations over each subinterval. Example S. The VSRC does not exclude linear multistep methods for the special second order differential equation y" = f(t, y) as long as the difference equation is written in a form that does not have bad roundoff properties. For example the summed form of Stormer's method, + h F u = u n n-1 n-1' F = F . + h f(t , u ), n n-1 n n satisfies the VSRC. 19 h. Optimally stable methods . We wish to characterize those methods that satisfy the most stringent definition of stability which is compatible with convergence. Spijker [6] has shown that among convergent straight multi- step methods which satisfy the RC there are some which are more stable than the others and that these optimally stable methods are those which satisfy the strict root condition (which for straight multistep methods is equivalent to the very strict root condition). Theorem U.2 extends Spijker's result by considering a much more general class of methods. In the proof of the theorem every optimally stable method is shown to be equally stable with a form of Euler ' s method . LEMMA k.l. If a method satisfies the VSRC and if a second method is at least as stable as the first method, then the second method must satisfy the VSRC and the E matrices of the two methods must satisfy E E = E . Proof. Call the first method Method 1 and the second method Method 2. By Definition 3.2 and Theorem 3.6 there is a constant K independent of N such that || U (2) - U (2) || .< K || [E J" 1 R || . Let m be an arbitrary fixed nonnegative integer, and let R = (iV, r_ , ..., T r , , .... ) . Then -m — — I "(2) (2) I yr max i „ , , \ u' - u £ K .. „ E. r.+ ...+r . +r ' -m -m ' £ n S m ' 1 — — n-1 -n for < h < min(h , l/(m + q - l)). Since the left-hand side is continuous in h, we can let h ■* to get (k.l) | ? n S^" n r UK* a * iLlr-.+ M. + r J+ r |. 1 n=0 2 -n ' <: n <: m ' 1-0 -n-1 -n ' This inequality holds for every nonnegative integer m. Choosing r = 0, n = l(l)m, shows that S must be power-bounded, which implies the RC 20 is satisfied. Suppose, however, Method 2 does not satisfy the SRC. Then there must be an eigenvalue £ / 1 of S of modulus 1 and a corresponding eigenvector v. Let .n-m r = £ v , n = l(l)m. — n — Then from (U.l) we get (m + 1) | v <: K _n-m -m max | £ =_£_ E v + 5 n-m £ n < m ' £ - 1 1- * - 2| E v | * k( i s - i r i 2 i } - Since this inequality cannot hold for every m, we conclude that the SRC must "be satisfied. Therefore Theorem 3.6 may he applied to Method 2, and so there is a constant c > such that (k.2) c || [Ej- 1 R || < || U V ^ ~ U V 2) TT (2 < K | | [E r 1 R Choose r = (E. - I) v, n = 0(l)N, where v is arbitrary. Then (k .2) becomes — m. 1 — — c | N E (E - I) v + (E - I) v | S K | (E - I) v | . Because both c and K are independent of N, E_(E -I) v = 0_, and because v is arbitrary, E 2 El = E 2 . Therefore rank E = rank E E £ rank E £ s showing that the VSRC is satisfied. Q.E.D. THEOREM k.2. A convergent method is optimally stable if and only if it satisfies the VSRC. Proof. Let a convergent method be given. The global discretization errors are Lq = ^y ; 0) " A y o + o(1) ' % = s i(y ; o) - a y Q + o(i). 21 Convergence implies that £(y;0) = A y, S o(y; 0) = A y for all y; hence S A = A. This means that the columns of A are linearly independent right eigenvectors of S corresponding to the eigenvalue 1. There- fore there exists a k x s matrix A of rank s such that T T A o s " V Define o^y; h) = A Q y, {km3) s o ■ V A o A o )_1 A o> ^(t, u; h) = A Q f(t, S Q u), v^(t; h) = A Q y(t). The method (4.3) satisfies the VSRC and is convergent by Theorem 3.9» Let us show that (4.3) is at least as stable as the original method. Setting 6 = u - u , we have that — n — n -n ^0 = V r=6-S6-h6 , n= l(l)N, ^i —n —n-l — n-± where 6 n $ L 6 _ I . Thus 1 -n-l ' ' —n-l ' Jo sn o" J *-r s l±o + ih S T 3 % - s i)-i - h ij-i>> 2 and since S = S and S S = S , I jo S V r -i I ■ I 4 + (s o - s > 4-i " h( ji s o" J Ij-i' I « (1 ♦ I S Q - S I + t n I Sq I L) Q ~, „ | ^ Therefore by Theorem 3.6 22 II U {0) -U (0) || (0. || [SJ- 1 R || II— — I'oo — °° .< C Q (1 + I S Q - S I + I S Q I L) || U - U II., which shows that (U.3) is at least as stable as the original method. Having established this result, let us show the necessity of the VSRC for optimal stability. Assume the original method is optimally stable. Then it must be at least as stable as (U.3), and so by Lemma U.l the original method must satisfy the VSRC. Next we show that the VSRC is sufficient for optimal stability. Assume the original method satisfies the VSRC. Let any other convergent method (Method 2) be given which is at least as stable as the original method. We must show that the original method is at least as stable as Method 2. From Lemma k.l it follows that Method 2 satisfies the VSRC and that E E = E . Hence the null space of E is a subspace of the null space of E . Because both methods are convergent, rank E = s = rank E; so the null spaces of E and E are identical. From E E = E it follows that E 2 (E E 2 - I) = 0, and since E and E have identical null spaces, which simplifies to Therefore E(E E 2 - I) = 0, E E = E. n „ i ™„ v n n-1 max | " n-j i max , " n-j , ) ji-j < n $ N ' j = E h ' $ n S N ' j=0 E 2 -j U E 2 J j = E 2 ^ < (1 + I V w M maX I ? F n ~ j r I $ (1 + | E 2 - E |) % n g N | .l Q E 2 r_. [, and so 23 II H - H I L * C || [E]" 1 R | 1^ $ C (1 + | E 2 - E |) | | [E 2 ] _1 R | | ^ C (1 + | E - E | ) c" 1 | | U (2) - U (2) || , showing that the original method is at least as stable as Method 2. Q.E.D. 2U 5. Order of convergence . In this section it is shown that the order of convergence is exactly equal to the order of quasi-consistency if the SRC is satisfied. Then examples are given of methods for which the order of quasi-consistency is greater than the order of consistency. DEFINITION 5.1. A method is said to be convergent of order p for the problem (2.1) if II U - I I loo = °( hP ) • If a method is convergent of order p for all infinitely smooth problems (2.1), it is simply said to be convergent of order p. DEFINITION 5-2. A method satisfying the RC is said to be quasi- consistent of order p with the problem (2.1) if .< n .< N I ^ I " ° (h } and "*i l«flo*S^"-*V l = °< hP > If this is true for all infinitely smooth problems (2.1), then the method is simply said to be quasi-consistent of order p. For a method to be quasi-consistent of order p with a problem it is sufficient that (i) d^ = 0(h P ) , n = 0(1)N (ii) E d^ = 0(h P+1 ) , n = 1(1)11. THEOREM 5.3. A method satisfying the SRC is convergent of order p for the problem (2.1) if and only if it is quasi- consistent of order p with (2.1). The proof is omitted since it is similar to that of Theorem 3.8. Example I. If the Adams-Bashforth-Moulton PEC method is written in the first way shown in Example 2 of §2, we have 25 E = 10 10 10 d = h" — n -5y< 3) /12 + o(h*), and E d = 0(h ), shoving that the method is quasi-consistent of order three ~ n Written in the second form, the method is also consistent of order three. Example 2. Butcher's three-stage procedure of order four has . , d = h k *° y r 1 "* E = 5f W V ^ 3> + 0(h 5 ) and E d = 0(h ). There seems to be no form of this method for which the — n order of consistency is four. Example 3. The modified linear multistep method has E = 1/2 5/6 1/2 1/6" " »i fc W k 1/2 5/6 1/2 1/6 d = — n = h 4 V>« + 0(h 5 ) and E d = 0(h P ). -n Example 4. The optimally stable two-cyclic method which alternates Milne's formula with the third order Adams -Moult on formula (see Stetter [7» p. 217]) is convergent of order four yet only consistent of order three. 26 6. Asymptotic behavior of the global discretization error . The principal term in the asymptotic expansion of the global error is given in Theorem 6.1 for methods satisfying the SRC. There is some simplification if, in fact, the VSRC is satisfied. This is discussed after the proof of the theorem. THEOREM 6.1. Assvme ij> (t, u; h) is uniformly Lipschitz continuous f" k in u and continuous in t and h for s t < 1, u^ R , and £ h £ h^. Assume ^q = h P y + o(h p ), ^ ■ h P lp(t n )+ h P+1 ^ +1 (t n ) + o(h P+1 ), n = 1(1)N, where (t) and j> -.(t) are continuous and E jL(t) = 0. 27zen u = y + h P (e(t ) + (I - S + E)" 1 6 (t ) + (S - E) n c) + o(h P ) — n j ti — n -*p n — uniformly for < n $ N where c_ = (I - E) y - (I - S + E)" 1 4(0) and the function e_ (t) satisfies £(0) = E Y e» = E G(t) (e + (I - S + E) A (t)) + E A^Ct) £■ = fi ti^"c ; ^_e + u - b + £J~ with G(t) = j^(t, A y(t); 0) Proof. Let y^ + h-(e(t n ) + (I - S + E)' 1 ip(t n ; u : + h P (e.(t_) + (I - S + E) 4_(0 + (S - E) n c_). Then we must have £o = Uq - i(y s h), r =u -Su , - h , n . i(t _, u _; h) s n = 1(1 )N. — n — n — n-1 - 1 - n-1 — n-1 It is easily verified that _c has been chosen such that r_ = o(h ). By the 27 mean value theorem and the Lipschitz continuity of $> r =u -Su - - h £ ( t -, » y -, ; h ) -n — n -n-1 *■ n-1 lL n-l ■* vv *»-i s h) ( --i - ^-i' + 0 + o(h P+1 - n — n — n-l — n-l n-l — = -h P+1 E G(t .) (S -E) 11 " 1 c + o(h P+1 ). n-1 — Therefore since E n (S - E) < °° n=l ' ' we must have n-1 P' .Z n E r. + r = o(h^), j=0 -j -n ' 28 which by Theorem 3.6 is equivalent to U - U | = o(h ). Q.E.D. In the case of a convergent method satisfying the VSRC we have that T E = A M . Differentiating (3-7) with respect to y yields M T 4 (t, A y; 0) A = f (t, y), and hence E G(t) e_(t) = E G(t) E e_(t) = A f (t, y(t)) M T e(t). T If we define e(t) = M e_(t), then the magnified error function is A e(t) + (I - S + E)" 1 6 (t) where e(0) = M T y and e' = f (t, y(t)) e + M T [G(t) (I - S + E)" 1 A (t) + 6 +1 (t)]. Note that the stability of the differential equation for e(t) depends only on the problem and not on the method. COROLLARY 6.2. Let the hypotheses of Theorem 6.1 be satisfied, and let < £ 1 he given. Then u = y + h P (e(t ) + (I - S + E)" 1 - By the remarks following Theorem 6.1 the magnified error function is e_(t) = e_ e(t) where e(t) solves the initial value problem + , e(0) = p(1) 1 [e n _ + (1 + a.) e q-1 1 q-2 e' = f (t, y(t)) e + pUT 1 ♦ +1 (t) + (1+« 1+ ... +a ) 6 Q ], assuming u. = y(t. ) + h P e. + 0(h P+1 ), i = 0(l)q - 1. Example 2. For Milne's method d = h' n r.-i ..(5) 1 2880 + 0(h°) , j»(t, u; 0) = ^(t, u; 0) = 1/3 2/3 2/3 1/3 1/3 2/3 2/3 1/3 f(t, u X ) f(t, u 2 ) f (t, u" 1 ) •J f (t, u ) G(t) = f (t, y(t)) 1/3 2/3 2/3 1/3 30 The magnified error function is e(t) where e' = f (t, y(t)) - y This system can be decoupled: e + 2880 1/3 2/3 2/3 1/3 1 1 1 -1 where (. 1 )--f y (t.x(t)) e 1 + *£$-, (e 2 )' = -i f (t, y(t)) e 2 . Example 3. For the modified linear multistep method equivalent to the fourth order Adams -Moult on method E = 1/2 5/6 1/2 1/6 1/2 5/6 1/2 1/6 1/U8' (I - S + E) -1 d = h — n k -1/U8 *™ * * 5 l(t, u; 0) = -1/80 17/720 3/8 1 -1/2U y (5) + h 5 n 1 1/8 -1/2U 0100 -11/2U 1 1/2U 0101 1/128 1/1+8 -1/1152 (k) f (t , y ) y y n n n + 0(h 6 ) , W + 1 1 u 2 , 1 3,1 Ux f(t, — u +u + 7f u + q- u ), 31 jjjt, u; 0) = 3/8 1 -1/21+ yt.K-» 2 -K + K>r|i Hi- G(t) = 3/8 1 -1/2U The magnified error function is f y (t,y(t)) [I i | I]. 1 1 e(t) + 1/U8 -1/1+8 y (u) (t) where and e(0) = I 1 + 1 2 + I y 3 + 1 J» e^o; " 2 Y l| + 6 Y l+ 2 Y l+ 6 Y U e« = f (t, y(t)) e + -^- y ( ^(t) + r4 f (t, y(t)) y U) (t) (M, 180 1+8 V This can be rewritten (e +-^y (1 °(t))' = f (t, y(t)) (e +J ^y {k) (t)) + ^y (5) (t 1+8 720 hence e(t) + rs- y . (t) satisfies the differential equation for the magnified error function of the fourth order Adams -Moulton method. 32 REFERENCES [l] B. Chartres and R. Stepleman, A general theory of convergence for numerical methods, SLAM J. Numer. Anal., 9 (1972), pp. U76-U92 . [2] C. W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice-Hall, Englewood Cliffs, N.J,, 1971. [3] P. Henrici, Discrete Variable Methods for Ordinary Differential Equations, John Wiley, New York, 1962. [h] P. Lancaster, Theory of Matrices , Academic Press, New York, 19&9. [5] A. Nordsieck, On the numerical integration of ordinary differential equations, Math. Comp., 16 (1962), pp. 22-1+9. [6] M. Spijker, On the structure of error estimates for finite difference methods, Numer. Math., 18 (l97l), pp. 73-100. [7] H. J. Stetter, Analysis of Discretization Methods for Ordinary Differential Equations, Springer-Verlag, New York, 1973. Form AEC-427 (6/68) AECM 3201 U.S. ATOMIC ENERGY COMMISSION U WIVE RSI TY-TYPE CONTRACTOR'S RECOMMENDATION FOR DISPOSITION OF SCIENTIFIC AND TECHNICAL DOCUMENT ( S*a Instructions on Rtwrm Side ) 1. AEC REPORT NO. COO-2383-0017 2. TITLE ANALYSIS OF FIXED-STEPSIZE METHODS 3. TYPE OF DOCUMENT (Check one): 1X1 a. Scientific and technical report I I b. Conference paper not to be published in a journal: Title of conference Date of conference Exact location of conference. Sponsoring organization □ c. Other (Specify) 4. RECOMMENDED ANNOUNCEMENT AND DISTRIBUTION (Check one): fX~l a. AEC's normal announcement and distribution procedures may be followed. I I b. Make available only within AEC and to AEC contractors and other U.S. Government agencies and their contractors. I | c. Make no announcement or distribution. 5. REASON FOR RECOMMENDED RESTRICTIONS: 6. SUBMITTED BY: NAME AND POSITION (Please print or type) C. W. Gear Professor and Principal Investigator Organization Signature wi>ce fj^ iW s- Date February, 1975 FOR AEC USE ONLY 7. AEC CONTRACT ADMINISTRATOR'S COMMENTS, IF ANY, ON ABOVE ANNOUNCEMENT AND DISTRIBUTION RECOMMENDATION: 8. PATENT CLEARANCE: □ a. AEC patent clearance has been granted by responsible AEC patent group. I~~l b. Report has been sent to responsible AEC patent group for clearance. f~~l c. Patent clearance not required. BIBLIOGRAPHIC DATA SHEET 1. Report No. UIUCDCS-R-75-702 3. Recipient's Accession No. 4. Title and Subtitle ANALYSIS OF FIXED-STEPSIZE METHODS 5. Report Date February 1975 |7. Author(s) Robert Skeel 8. Performing Organization Rept. No - UIUCDCS-R-75-702 |9. Performing Organization Name and Address Department of Computer Science University of Illinois at Urbana-Champaign Urbana, Illinois 6l801 10. Project/Task/Work Unit No. 11. Contract/Grant No. US AEC AT (ll-l) 2383 12. Sponsoring Organization Name and Address US Atomic Energy Commission Chicago Operations Office 9800 South Cass Avenue Argonne. Illinois 60^39 13. Type of Report & Period Covered 14. 15. Supplementary Notes 16. Abstracts The unified theory developed by Henrici [3] for one-step methods is extended to the more general case where the order of the system of difference equations can exceed the order of the system of ordinary differential equations. The analysis is applicable to every fixed-stepsize fixed-formula method known to the author. For many of these methods the concept of consistency is inadequate. A more appropriate concept, termed quasi-consistency, is introduced. 17. Key Words and Document Analysis. 17a. Descriptors 7b. Identifiers /Open-Ended Terms 7c. COSATI Fie Id /Group 18. Availability Statement Unlimited ORM NTIS-35 (10-70) 19. Security Class (This Report) UNCLASSIFIED 20. Security Class (This Page UNCLASSIFIED 21. No. of Pages 22. Price USCOMM-DC 40329-P71 ^ •gfr I ty H linos! h33 RIB MkflHUH Dull mm fflmn Mil m m flB^ aura Hm HI ^^B^H HK9 rISE Wffl m H mm ram ■ HffiSBl Hi ffinSH HH| iillii HffimBH liiiiiiSsiiHi ^n HHm HnltfllMMl lev, HI HK81 :i^.:E&l*«