UNIVERSITY OF ILLINQ.3 LIBRARY AT URBANA-CHAMPAIGN - '• k At ljYuCDCS-R-77-899 UILU-ENG 77 1757 C00-2383-0046 A SEARCH FOR BETTER LINEAR MULTISTEP METHODS FOR STIFF PROBLEMS by Antony King-Yin Kong December 1977 Digitized by the Internet Archive in 2013 http://archive.org/details/searchforbetterl899kong 5 1 v. o^ in UIUCDCS-R-77-899 A SEARCH FOR BETTER LINEAR MULTISTEP METHODS FOR STIFF PROBLEMS* by Antony King-Yin Kong December 1977 DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN URBANA, ILLINOIS 61801 * Supported in part by the US Energy Research and Development Administration under contract US ERDA/EY-76-S-02-2383, and submitted in partial fulfillment of the requirements of the Graduate College for the degree of Doctor of Philosophy. 111 ACKNOWLEDGEMENTS The author wishes to express his gratitude to Professor Robert D. Skeel not only for his ideas and guidance, but also for his many valuable suggestions and criticisms during the course of this project. He is also grateful to Professor Charles W. Gear, Professor Daniel S. Watanabe, and other members of the department for ideas and insights which contributed much to this work. The author is also indebted to Barbara Armstrong for her speedy and accurate typing of this thesis. The work was supported in part by the Energy Research and Development Agency under grant US ERDA/ EY-76-S-02-2383. Vlll TABLE OF CONTENTS CHAPTER Page I INTRODUCTION 1 1.1 STIFF EQUATIONS 1 1.2 PROPERTIES OF NUMERICAL METHODS FOR STIFF SYSTEMS 3 1.3 LINEAR MULTISTEP METHODS 6 1.4 OBJECTIVE 13 II MEASURE OF ACCURACY 16 2.1 INTRODUCTION 16 2.2 ASYMPTOTIC GLOBAL ERROR 16 2.3 GLOBAL ERROR BOUND 17 2.4 CONCLUDING REMARKS 22 III FROM A -STABILITY TO A(a)-STABILITY 24 3.1 INTRODUCTION 24 3.2 CHARACTERIZATION AND PROPERTIES OF A Q -STABLE METHODS. ... 25 3.3 SPECIAL CASES k = 1,2 28 3.3.1 k = 1 28 3.3.2 k = 2 29 3.4 SOME LOWER BOUNDS 31 3.5 HIGHER ORDER A(a)-STABLE METHODS 36 3.6 CONCLUDING REMARKS 45 IV MINIMAX 47 4.1 INTRODUCTION 47 4.2 STATEMENT OF THE PROBLEM 47 4.3 A FEASIBLE DESCENT ALGORITHM 49 4.4 CONVERGENCE 56 IX 4.5 MODIFICATIONS 63 4.6 CONCLUDING REMARKS 66 V ACCURACY VS A(a)-STABILITY 68 5.1 DESCRIPTION OF THE PROBLEM 68 5.2 MINIMAX FORMULATION OF THE PROBLEM 70 5.3 NUMERICAL RESULTS 73 5.4 CONCLUSION 87 BIBLIOGRAPHY 89 APPENDICES A POLYNOMIALS HAVING ROOTS IN THE OPEN LEFT HALF COMPLEX PLANE . . 92 B LIST OF A, A*(a), (rG) 1/k , r VALUES FOR METHODS FOUND USING THE M-ALGORITHM 96 VITA 101 CHAPTER I INTRODUCTION 1.1 STIFF EQUATIONS Systems of ordinary differential equations which arise in the description of physical problems often have greatly differing time con- stants and are very stable (Bjurel et al., 1970). They are known as stiff equations. A model stiff equation is given by y'(t) = Ay(t) y(t Q ) = n where yeR , A is an mxm real matrix with all its eigenvalues A-,, x 2 » lying in the open left half complex plane, and the ratio ••*. max Re A. i,j | Re A is large. All solutions y(t) = exp[(t-t n )A]n approach zero as t increases, hence the system is asymptotically stable about y = 0. The component of y(t) corresponding to an eigenvalue a. whose real part is large in magnitu- de will soon become negligible relative to the dominant component, i.e., the one corresponding to the eigenvalue whose real part is the smallest in magnitude, which could still be large enough to be of interest. For exam- ple, we consider the following equation when m = 2: ~y\{t)~ yj(t) ' A-1 -A-1 y } (t) -y,(t o r n i -A-1 A-1 ' _i y 2 (t) _y 2 (t )_ _ n 2_ where A << 0. The true solution is -t At y-j(t) = c^e + c 2 e -t At y 2 (t) = c^e - c 2 e c-| = (n-| + n 2 )/2 Co = (rin - n 9 )/2 The component exp(xt) is negligible after t >_ 1 , and thus we would like to be able to use a larger stepsize h from that time on when integrating the system numerically, since we need only follow the dominant component exp(-t). However, numerical stability of methods which are not specifical- ly designed for stiff equations, such as the one-step Euler's method, re- quires that the value, for the m-dimensional case, max | hx j | l such that a perturba- tion in the equations defining y by a fixed amount produces a bounded change in the numerical solution when any he(0,h fi ] is used as stepsize. A more precise definition of stability is given in Stetter (1973, p. 9). Both the concepts of convergence and stability, necessary for any method to be useful, are concerned with the limiting process as h ■> 0. In practice, we are more concerned with the effect of error accu- mulation -- not only roundoff error as introduced earlier, but also discre- tization error from approximating the true solution -- when some finite stepsize h is used. We thus introduce the concept of absolute stability for a given fixed stepsize. Since the concept is problem dependent, we de- fine absolute stability for the class of differential equation y'(t) = xy(t) y(t Q ) = n (1.2) where x is a complex constant. Definition The region of absolute stability 4 of a method is the set of all hx such that when the method is applied with stepsize h > to equation (1.2), the numerical solution y + as n + ». The method is then absolutely stable for stepsize h and for the equation (1.2) if hxe,4. The following terminology has been used to descri- be methods according to the extent of their region of absolute stability: Definition (Cryer, 1973) A method is A n -stable if the negative real axis {z : Im z = 0, Re z < 0} is in A. Definition (Widlund, 1967) A method is A(a)-stable, ae(0,V2), if the wedge S = (z : |Arg(-z) | < a, z f 0} is in A. The a angle of the largest such wedge is called the angle of absolute stability associated with the method. A method is A(0)-stable if it is A(a)-stable for some ae(0,-rr/2), and is A(7r/2)-stable if it is A(a)-stable for all ae(0,ir/2). Definition (Dahlquist, 1963) A method is A-stable if its region of absolute stability contains the open left half complex plane. From the definition, it is obvious that A-stability is equivalent to A(tt/2) -stability. Moreover, A-stability => A(a)-stability => A(0)-stability * A Q -stability Methods which are A n -stable shall be loosely called stiff methods. The angle of absolute stability is only one of a number of parameters which have been proposed for measuring the extent of the region of absolute stability A (cf. Gear, 1971, p. 213, Odeh and Liniger, 1971, Gupta, 1976). But it is probably the best such measure, especially for methods with automatic steps ize selection. When such methods are applied to the equation y'(t) = xy(t) + g(t) x complex the integration starts out with hx near the origin, and as the integration proceeds, the stepsize h increases and the value hx moves away from the origin. However if hx approaches the boundary of A, this would be detected by the error estimator, and any further movement of hx would be prevented. The implication of this is that not all of A is "used", but only that por- tion which can be reached from the origin by rays lying entirely inside A. The most "desirable" portion, where any h > can be used, is described by the angle of absolute stability. This parameter, a, also serves as a good indicator of the problem class for which the method is suitable, namely, those problems for which the large eigenvalues of the Jacobian lie inside the wedge S a . Hence, in an ODE solver, « can be an extra parameter which is used to identify among a family of methods of order k the A(a)-stable method that should be used. The value of a can be either supplied by the user or chosen automatically by some detecting device in the code. It would then be desirable to have A(a)-stable methods of any order for any oe(0,ir/2]. Since the true solution of (1.2) for any A lying in the open left half plane converges to zero as t ->■ °°, we would prefer methods whose region of absolute stability contains the open left half plane, i.e., A-stable, so that the numerical solution y n also converges to zero as n ■* °° for any stepsize h > 0. On the other hand, if A is not in the open left half plane, y(t) will not converge to zero, and in fact will diverge if Re A > 0. However, if hAGA, y n converges to zero and consequently invalid solu- tions are produced (Lindberg, 1974). Hence, methods whose region of abso- lute stability is equal to the open left half plane are desirable. We end this section with another definition: Definition (Odeh and Liniger, 1971) A method is A^-s table if its region of absolute stability con- tains a neighborhood of infinity. A^-stability is desirable especially for those problems which have a slowly varying component and a rapidly decaying component, so that transient that decays to zero rapidly in the true solution would not be present in the numerical solution as slowly dampened oscillatory com- ponents. 1.3 LINEAR MULTISTEP METHODS We restrict our study to the class of linear k-step formulas j! ¥n + j " h jf 6 / is the stepsize assumed to be fixed throughout the range of in- tegration. When k > 1, the first step of integration requires numerical values for y »y-i ^k-T ^ nce we are on ^ given n as an initial value, the k values have to be computed from r\. A common technique is to use a one-step method such as the Runge-Kutta method (Gear, 1971, Chapter 2) to compute y Q ,y-| ,... ,y k -. . The procedure used to compute y n ,y-,,.. . ,y. -, from n is called a starting procedure. The k-step formula (1.3) together with a starting procedure constitutes a k-step method. The starting procedure is said to be consistent with (1.1) if for j = 0,l,...,k-l, Vj = y.(n.h) -> n as h -> 0. At each step of the integration, in the case when e =(= 0, i.e., when the formula is implicit, we have to solve a system of nonlinear equa- tions for y . . For stiff equations, this is usually accomplished by some kind of modified Newton iteration. Results are usually satisfactory if the stepsize h is small enough and if the initial estimate is sufficiently ac- curate (Lambert, 1973, p. 239). The initial estimate could be obtained by using an explicit formula. 8 Unless explicitly specified, the following assumptions about the linear k-step formula (1.3) are made: k (A) z 3- = 1» a normalization. (1.4) j=0 J (B) The formula (1.3) satisfies the root condition, i.e., all the roots of the polynomial k j p(C) = l a.C (1.5) j=0 J lie inside the unit circle Ul l 1 and those on |c| = 1 are sim- ple. Occasionally, we require that (1.3) satisfies the strict root condition, i.e., all roots of p(c) lie inside |c| < 1 except for a root at ? = 1. (C) The order of the formula (1.3) is k, where formula (1.3) is said to be of order p if c = c, = ... = c =0in the following Taylor series k £ Dv(t+jh) - he y'(t+jh)] = j=0 J J = c Q y(t) + Cl hy'(t) + ... + c p h p y (p) (t) + +Cp+1 hP +1 y (p+1) (t)+... (1.6) It is said to be of exact order p if c n = c, = ... = c =0 and 1 p c ., 4 0, in which case, we call c ., the error constant of the p+1 p+1 formula. Formulas of order at least 1 are said to be consistent. (D) The polynomials p(c) and a(c) have no common factors, where p(c) is as defined in (1.5) and k i o(c) = s e,c J (1.7) j=0 J Stetter (1973, pp. 186-187) shows that if the polynomials p and a of a linear multistep formula possess one or more common factors, there exists an equivalent formula with a smaller step number in the sense that they generate the same numerical solution from given starting values. We also assume that the starting procedure is consistent. The k-step method is said to be consistent if both the starting procedure and the k- step formula are consistent. It is known that a k-step method is stable if and only if the k-step formula satisfies the root condition (Stetter, 1973, p. 208). It has also been proved that a k-step method is convergent if and only if it is stable and consistent (Stetter, 1973, pp.13, 211-213 -- note that some of the details are missing). Many authors require that the formula (1.3) together with an ar- bitrary consistent starting procedure be convergent (e.g., Henrici, 1962, p. 218). However, Stetter (1973, p. 78) indicates that, except in pathologi- cal cases, this requirement is no stronger than ours. It should also be remarked that the assumptions about constant stepsize and order are just idealizations that make our analysis tractable. In practice, both the stepsize and order will vary. Also, in practice, formulas that satisfy the strict root condition are more desirable. We call those formulas strongly stable. We denote a k-step method by (p,a) where p is as defined in (1.5) and a in (1.7), and use the symbol L[k] to represent the class of (p,a) which satisfy conditions (A)-(D). It is easily seen that, under assump- 10 tions (A) and (C), we have a (1) = 1 c ■ p(1) = C 1 = p'(l) - a(l) = and, in general, for q = 2,3,... Ik 1 c = - £ (j q a 1 - qj q " Bi) q q! j-1 J J We remark that in order for a method (p,a) in L[k] to be Aq- stable, it is necessary that (Cryer, 1973) order <_ k except the trapezoidal rule Furthermore, the region of absolute stability A associated with a method (p,o) can be described by A = (yeC : Uj(y)| < 1, j = l,2,...,k> where £,-(y), j = 1,2 k, are the k roots of p(0-u a (c), depending on the value of yeC (Lambert, 1973, p. 222). The stability condition requires that all the roots of a po- lynomial lie inside the unit circle. Since there exists a reasonably sim- ple algebraic condition for all the roots to lie in the left half complex plane, we wish to transform the unit circle to the left half plane. This can be done by using the bilinear transformation z+1 C = — (1.8) z-1 which maps the left half plane in complex z-space 1-1 onto the unit circle in complex c-space. Apply the transformation (1.8) to the polynomials p(c) 11 and a(0, and define the corresponding polynomials in z by k .z+1 r(z) = z a.z J = (z-1) k p( — j=0 J z-1 k , .z+1 s(z) = z b.z J = (z-l) k a( ) J-0 J z-1 It is easily seen that a k = p(D = b k = o(l) = 1 Furthermore, a method is of order k if, and only if, the following holds: r(z) z+1 2 k+1 c 1 log = _p + O(-j-rp) asz.o, s(z) z-1 z k+1 z k+2 (obtained from the definition of order by setting y(t) = exp(t) in (1.6)). Therefore, requiring that the method be of order k imposes k additional linear conditions on the a.'s and the b.'s: J J k b 1 a- = 2 s — - j = 0,1 k-1 (1.9) J 1-J+l 1-J i-j odd Hence, the method is uniquely parameterized by the k parameters b ,b-|,... ,b k _-| (b k = 1). Occasionally, we denote such parameterization by (bQ,b-|,... ,b k _-j), or by the polynomial s(z), or by (r,s) in analogy to (p.tr). The parameterizations (p,a) and (r,s) are related by the fol- lowing linear transformation (Duff in, 1969): a = Ma b = MB where a = [a ,a.j,... ,a k ] , b = [b^b^... ,b k ] , a = [a Q ,a 1 a k ] , 6 = 12 [Bq,B-j ,... ,6 k ] , and M = [M. .] is a (k+l)x(k+l) matrix whose elements are given by k-i k M. = (-D k \) i = 0,1,. ...k M kj =1 j = 0,1,... ,k with M" 1 = 2~ k M. We note that the polynomial r(z) of a stable method has roots all in the left half complex plane, with at most simple roots on the imaginary axis, and likewise for the polynomial s(z) of an A(0)-stable method (Widlund, 1967) -- the condition is relaxed for A n -stable methods whose s(z) could have at most double roots on the imaginary axis (Cryer, 1973). The case when all the roots of r(z) ( s(z) ) lie in the open left half plane corresponds to the method being strongly stable ( A^-stable ). The parameterization (r,s) has several advantages over (p,a). The restriction that (r,s) be of order k imposes a simple relation on the a.'s and the b.'s (cf. (1.9)). Moreover, the error constant c k+ i can be expressed simply as a linear combination of the even b-'s (cf. (2.3)). Also, the Hurwitz criterion (cf. Appendix A) for polynomials with roots in the open left half plane is "simpler" than the Schur criterion (Marden, 1966, p. 198) for polynomials having roots inside the unit circle. However, the a-'s and the b^'s have no natural interpretation like the a.'s and the B-'s have, namely, parameters of the linear difference equation defining the numerical solution y n< We finally remark that several versions of the parameterization (r,s) have been used by different authors. The one described in this sec- 13 tion is used in Liniger (1975), Genin (1973), and Jeltsch (1976). Dahlquist (1963), Widlund (1967), and Cryer (1973) use the same bilinear transformation (1.8) but define the polynomials r(z) and s(z) differently: z-1 v z+1 r(z) = (— ) P(— ) 2 z-1 z-1 k z+1 s(z) = (— )o( ) 2 z-1 Henrici (1962, p. 230) and Gear (1971, p. 195) use a different bilinear transformation and hence different r(z), s(z): 1+z c = 1-z 1-z k 1+z r(z) = (— ) k P (— ) 2 1-z 1-z k 1+z S(z) = (— ) k a(— ) 2 1-z 1.4 OBJECTIVE As pointed out in section 1.2, in order to be able to solve gen- eral stiff systems of the form (1.1), we would very much prefer to use methods which are A-stable, so that they would be suitable for all stiff problems. However, Dahlquist (1963) proves that the maximum order for A- stable linear multistep methods is only 2. To allow for methods having or- der greater than 2, we have to relax our stability requirement and consider A(0)-stable methods. Cryer (1973) shows that there exist A -stable methods of arbitrary order k by explicitly constructing a class of Ag-stable methods in L[k], which are later shown by Jeltsch (1976) to be in fact A(0)-stable. However, the error constant, which is an asymptotic measure 14 of the accuracy of the method (cf. Chapter II), of Oyer's methods becomes astronomical in magnitude when k is large, say when k = 10. Consequently, we are compelled to find more accurate A(0)-stable methods. Work has been done in seeking stiff methods with as high order as possible. Widlund (19G7) shows that A(a)-stable methods exist for any ae(0,iT/2) for the cases when k = 3,4. Using an interactive computer gra- phics program, Dill and Gear (1971) find stiffly stable methods (definition in Gear (1971, p. 213) -- note that stiff stability implies A(O)-stability) of orders 7 and 8. They only consider the class of methods most of whose B.'s being zero. However, Skeel (1977) shows that eyery k-step method is a (k+l)-value method. The same observation is mentioned in Wallace and Gupta (1973) and is proved, for the special case when the method is of order k+1 , in Osborne (1966). Hence insisting that most of the $.'s be zero would not reduce the number of previous values needed to be stored. Jain and Srivas- tava (1970) consider the polynomial o(c) = C k_r (c-c) r 0 0, what is the limitation on the angle of absolute stability for methods having error constant of mag- nitude C ? (iv) For any as(0,7T/2), do there exist methods which are A(a)-stable ? If they exist, how accurate are they ? Chapter II is devoted to a discussion on measures of accuracy. The answers for the above questions for the cases when k = 1 ,2 are given in Chapter III, where questions (i) and (iv) for a general k are also con- sidered. An algorithm for solving minimax problems numerically is introdu- ced in Chapter IV. Questions (ii) and (iii) for a general k are formulated as a minimax problem in Chapter V and some of the numerical results are listed. 16 CHAPTER II MEASURE OF ACCURACY 2.1 INTRODUCTION As stated in section 1.4, we are investigating limitations on the accuracy of stiff methods (definition in section 1.2). The results of the investigation are meaningful only if a reasonable measure of accuracy is used. For all practical purposes, we would say that numerical method A is more accurate than method B if the numerical results from using method A are closer approximations to the true solution than those from using method B for the same stepsize. Consequently, the order of the method alone is too crude as a measure of accuracy; rather a more refined measure of global error is needed. Since the actual global error accumulated will depend on the differential equation to be solved, parameters which depend only on the numerical method in either a global error bound or in an expression for the asymptotic global error will be suitable as a measure of how accurate the method is. Two candidates for the measure of accuracy are introduced in sec- tion 2.2 and 2.3, respectively. To simplify expressions for global error, we assume, for the rest of this chapter, that the starting values are ex- act, and that the roundoff errors are negligible compared with the discre- tization errors. 2.2 ASYMPTOTIC GLOBAL ERROR We first consider a measure of accuracy which indicates the mag- nitude of the asymptotic global error, i.e., the magnitude of the most dom- inant term as the stepsize h becomes arbitrarily small (Gear, 1971, pp. 75-6, 204-5): 17 y n - y(t n ) = h k c k+1 | n K(t n ,T)y (k+1) (T)dx + 0(h k+1 ) h - (2.1) where K(t,t) satisfies the homogeneous ordinary differential equation 3K df — (f.x) = — (t,y(t))K(t,x) K(t.t) = I 9t ay Thus, one measure of accuracy is the parameter c. + -., which is the only method dependent part of the leading term of the asymptotic error expan- sion. The error constant can be expressed in terms of the a.'s and the 5.'s as J C k+1 = 2 [j k+1 a - (k+l)j k 3.] (2.2) k ' (k+1)! j=l J J and in terms of the b.'s as (Genin, 1973) c k+l " 1 k b, E — ^- 2 k j=0 j+1 j even = 1 f 1 " 2 k+1 i-l S(Z: (2.3) 2.3 GLOBAL ERROR BOUND Under the assumption that the stepsize h used when applying (p, a) to solve system (1.1) is small enough such that h|%.0 J a Q C +a-,C +...+a k j = J k B = E j-0 J .k G = ' ' r k = | |G(s) |ds (2.6) J n G(s) is the influence function associated with (p,a) and is defined by G(s) = — E [a.(j-s) k - ke.(j-s)^" 1 ] (2.7) k! j=0 J J where, for any ve(-«,»), Note that a method of order > k is exact for polynomials of degree k and hence G(s) = for s < 0. Clearly, by definition, G(s) = for s > k. It is easily shown using (1.6) that G(s) is the kernel of the operator L h (Henrici, 1962, p. 248): L h y(t) = h k+1 f G(s)y (k+1) (t+sh)ds (2.8) J— oo where l_ h is the linear difference operator defined by 19 k k L h y(t) = z 0i y(t+lh) - h s s^' (t+1h) (2.9) n 1=0 i=0 Henrici (1962, p. 242) proves that r is finite. A refinement of the above bound can be derived by using the quan- tity k k-1 e_r+s,c K '+...+B. °° A , f - max Iy^I. ° ' c i = E V ( 2 - 10 ) j>0 J a C K +a 1 C K '' + ...+a |< j=0 J instead of rB, as follows: We start with the difference equation satisfied by the numerical solution, k k z Vi-k+i ' h z i=0 n J K n i=0 The difference of the above equation and (2.9), the latter evaluated at t = ,L °i y j-k+i " h ,L 6 i f(t j-k + i^j-k + i) " ° J i k t. k , yields k i= a i e j-k+i ■ h i= Bi§ J- k+i = " L h y(t J-k ) j ^ k where e = y - y(t ) and e = f(t ,y ) - f(t ,y(t )) for all v > 0. In V V V V V V V v — particular, e = e = for £ v £ k-1. Multiplying the above equation by y . and summing from j = k to j = n, n-J aw fa Y n-j fa a i e j-k + i - h fa Vj fa »1*J-kH ■ " fa 'n-jV'Vk' The first term is, by (2.5), simply e , whereas the second term can be shown, using (2.10), to be 20 h A 9 "-j'j Hence, we have n-k E j=k "- J J j=0 As in Henrici (1962, p. 247), L.y(t) can be uniformly bounded by e n - h ,l Vj 5 j - .:„ ^v(Vk-j) t 2 - 11 ' I L h y ^^ I - hk+lQY t o - l < t+kh - T From (2.11), by the triangle inequality, n_1 - k+1 |e n | < hfL e |e.| + h|Y IM e n l + n rGY(n-k+l) or hrL n-1 h k+1 rGY le I < e le.l * (n-k+1) n _ l-h|Y |L j=k J l-h| Yo |L if h| Y() |L < 1. We shall make use of the discrete Bellman inequality as stated in the following lemma (Babuska, Prager, Vitasek, 1966, p. 57): Lemma 2.1 and Let a , \i> ,9 be sequences defined for v = 0,1 ,... ,n, e > 0, V V v v — v-1 < i/> >t + E 9 J>^ v = 0,1 ,... ,n •j=0 Then 21 n-1 i *nl*n + s n Vy n ( 1+e s> Applying the above lemma with = |e +k _-j | , and k+1 h rGY hrL \\) = v, e = l-h|Y lL l-h|Y IL we obtain the following global error bound: (Vt^L exp[ ■ ] - 1 |y n - y(t n )| 1 h K rGY ^ (2.12) rL The bound (2.12) is better than Henrici's bound (2.4) since f <_ rB The inequality is usually strict. For example, for the 2-step Adams Bash- forth method (Gear, 1971, p. 109), f = 3/2 rB = 2 From the bound (2.12), one possible measure of accuracy is the quantity rG. There is no nice expression for rG in terms of the b.'s. We remark that the bound (2.12) can be further improved by using the fact that, from (2.8), n-k k+1 k n-k I E Y i Lh^n-k-i)! - h Y I I z Y.GCs+jJlds j=0 J n K J J k-n j=0 J If we define 1 r k n-k x= max I | z Y.:G(s+j)|ds k max z ( )(l-e) v = ., , j>0 v=0 k-2 e K_l Hence, the measure rG (or x) is preferable in that sense. 24 CHAPTER III FROM A Q -STABILITY TO A(a)-STABILITY 3.1 INTRODUCTION As pointed out in section 1.4, the reason for considering A Q - stability and A(a)-stability, which are weaker than A-stability, is to al- low for methods of order greater than 2. In proving the existence of A Q -stable methods in L[k] (definition in section 1.3) for arbitrary k, Cryer (1973) shows that the following methods, as parameterized by the po- lynomial s(z) (cf. section 1.3), s(z) = (z + d) k k+1 are A Q -stable for d _> 2 . Jeltsch (1976) later proves that these methods are in fact A(0)-stable. Hence A(0)-stable k-step methods of order k exist for arbitrary k. The error constant of Cryer 1 s method is given by (cf. (2.3)) 1 k k d k " J d k c k+l = " T E ( ^ < ' ~k ' 2 J-0 j j+1 2 K j even k 2 and thus |c. + -,| is greater than 2 . Although, as remarked by Cryer, the bound for d could be strengthened, yet the essential point is that c^.-, is of order 0(d ), and that as k increases so must d. It is the impractica- bility of Cryer's methods that motivates our investigation of limitations on the accuracy of A Q -stable methods. Some characterizations and properties of A Q -stable methods are listed in section 3.2. In section 3.3, we discuss in detail the accuracy of A Q -stable methods in l[1] and l[2]. Some lower bounds on |c^ + i| for Zb AQ-stable methods in L[k], for arbitrary k, are derived in section 3.4. The existence of A(a)-stable methods of arbitrary order k for any ae(0,-rr/2) is established in section 3.5. 3.2 CHARACTERIZATION AND PROPERTIES OF A Q -STABLE METHODS This section contains a survey of known results concerning A - stability which will be referred to in subsequent sections and also in Chapter V. The parameterization (r,s) for methods in L[k] (cf. section 1.3) is used throughout. We characterize A -stable methods in terms of some algebraic con- ditions on the polynomials r(z) and s(z). Restrictions upon the coeffi- cients b.'s are also given. All methods are assumed to be in -L[k], For any complex number z, Arg z denotes that value of arg z which is in the in- terval (-tt.it]. We first consider A n -stability. The following theorems are taken from Cryer (1973): Theorem 3.1 The following statements are equivalent: (i) (r,s) is A Q -stable. (ii) For all qe(-°°,0), the roots of r(z)-qs(z) lie in the interior of the left half complex plane, (iii) For Re z > 0, r(z)/s(z) is regular, and does not take values in (-<»,0). For z = iw, i = /^T, we(-°°,°°), r(iw)s(-iw) does not take values in (-°°,0). Theorem 3.2 If (r,s) is A n -stable, then the zeros of r(z) and s(z) lie in the 26 closed left half plane, and those on the imaginary axis are at most double zeros. Theorem 3.3 If (r,s) is Ag-stable, then b 10, b] 1 and, if k > 3, bj > for j = 2,3,...,k-l A characterization of A(a)-stable methods is given by Widlund (1967): Theorem 3.4 The following statements are equivalent: (i ) (r,s) is A(a)-stable. (ii) For all qeS a = , the roots of r(z)-qs(z) lie in the interior of the left half plane, (iii) For Re z > 0, r(z)/s(z) is regular, and takes its values in the com- plement of S . a If, in addition to A(a)-stability, we assume that (r,s) is A^- stable, then r(z)/s(z) is regular on Re z > 0. Thus, by continuity, A(a)-stability implies that r(iw)/s(iw) ^ S for all w€( -«,<») (3.1) On the other hand, the locus r(iw)/s(iw), we(-°°,°°), of an A^-stable method divides the complex q-plane into several components, one of which contains a neighborhood of °°, denoted by V. Obviously, if (3.1) is satisfied, V 27 contains S . Since r(iw)/s(iw), wet- 00 , 00 ), is the locus of q in the complex plane for which a root of r(z)-qs(z) is purely imaginary, all roots of r(z)-qs(z) for any q in the interior of V have to lie in the interior of the left half complex plane, a condition satisfied by r(z)-qs(z) at q = °°. From theorem 3.4, this means that the method is A(a)-stable. Hence, for methods which are A^-stable, (3.1) is necessary and sufficient for A(a)- stability. It can easily be shown that (3.1) is equivalent to the following condition: Re r(iw)/s(iw) - < cot a (3.2) | Im r(iw)/s(iw) | for all we(-°°,°°) such that r(iw) j 0. Condition (3.1) or (3.2) will be used in subsequent sections and in Chapter V to prove A(°0-stability when the method is A^-stable. We conclude this section by giving another necessary and suffi- cient condition for A(a)-stability if the method is known to be Ag-stable instead of A^-stable. The condition is a direct consequence of a theorem by Jeltsch (1976), rewritten using the parameterization (r,s) instead of (p,a): Theorem 3.5 The conditions (i)-(iv) are necessary and sufficient for a method in L[k] to be A(0)-stable: (i) (r,s) is A -stable. (ii) The roots of s(z) on the imaginary axis are simple, (iii) Let z = iw be a purely imaginary root of r(z), then 28 s(iw) Re > r'(iw) (iv) Let z = iw be a purely imaginary root of s(z), then r(iw) Re > s' (iw) If the method (r,s) is Ag-stable, then it is A(a)-stable if and only if conditions (ii)-(iv) in the above theorem hold, and (3.2) holds for all we(-oo,a>) such that r(iw) j 0. 3.3 SPECIAL CASES k = 1 ,2 The limitations of the accuracy of A^-stable methods in £[1] and L[2] are discussed in this section. A thorough analysis is possible since only a few parameters are involved. 3.3.1 k = 1 The polynomials r(z), s(z), and the error constant c 2 of a method in£[l], parameterized by (b fi ), are r(z) = 2 s(z) = z + b Q c 2 = - b /2 Since r(iw) 2b Q 2w = -s — p " 1_ p — ? f° r a "ll we(-°°,°°) s(iw) w +b Q w +b o the method (b Q ) is A-stable if and only if b Q > 0. If b Q = 0, we have the trapezoidal rule, which is in fact of order 2 since c ? = 0. Its region of 29 absolute stability is the entire open left half plane. Any method corresponding to b fi > is A^-stable. Hence, A-stable methods exist for any |c«| 1 0» and, moreover, those A-stable methods are A^-stable if |c«| > 0. The polynomials p(c) and a(c) corresponding to the method (b fl ) are p(c) = c - 1 l+b„ 2 aM . ^ + ^o It is obvious that y. = 1 for all j >^ 0, and that y Q = (l+b«)/2, and y. = 1 for all j > 1 , so that r = 1 f = max{(l+b )/2,l} Furthermore, the influence function is given by 1-b, '0 G(s) =< 2 and it is easy to show that - s <_ s <_ 1 otherwise 'I** G = x =< < b Q < 1 1 0, b 1 > Thus there exists no A Q -stable method in l[2] with |cJ < 1/12. The po- lynomial r(z) + qs(z) = qz 2 + (qb^z + (qb Q +2b 1 ) has roots in the open left half plane for any qe(0,°°) if b Q , b-, satisfy b Z. °» b i 1 °» b o +b l > ° ( 3 - 4 ) Thus, from theorem 3.1, any method (b Q ,b, ) which satisfies (3.4) is A Q - stable. If b Q = 0, b, > 0, r(z) and s(z) have a common factor and the corresponding method reduces to the one-step trapezoidal rule. Note that it is the only Ag-stable method whose order (2) is greater than its step number (1) (Cryer, 1973). From (3.3), r(iw) 2 Dl b Re s(iw) ( b n" w ) +b i w 31 is nonnegative for any we(-o°,°°) if b Q ^ 0, and b, >_ 0. Hence for any |cJ > 1/12, we can construct A-stable methods in l[2] by choosing b = 4 | c | - 1/3, b-, >_ (cf. (3.2) and the remarks after theorem 3.5). If b. > 0, the method is A -stable. oo The polynomials p(c) and a(c) corresponding to (b n ,b.) are b,+l 9 b,-l b,+l b,-l p(0 = J — C - b,c + - 1 - ■ -M? - l)(c - — ) '1 b Q+ b 1+ l 1-b b Q -b 1+ l o(c) = £ + ? + yi so that and V 1 i+1 Y = i . (-! — )J +I for all j > J b-,+1 n-b, b — L s(s - 1 +-=-) G(s) =/ 4 1+b, b r l J ( S . 2 )(s - 1 + --*-) 4 1+b, ^0 < s < 1 1 < s < 2 otherwise For fixed b^ >_ 0, if we choose <_ b, <_ l+b Q , then G(s) is of one sign on [0,2] and hence G = \c^\ , and if we choose 1 <_ b,, then r <_ 1. In either case, the method is A-stable, and is also A -stable if b n > and b n > 0. » oo o 1 3.4 SOME LOWER BOUNDS As pointed out in section 2.2, the error constant c. . associated with a method in z,[k], parameterized by (b^.b-, .... ,b. , ), can be expressed as a linear combination of the even b-'s (cf. (2.3)): J 32 1 k b. c k+1 = - -f 2 _i k ' 2 J-0 j+1 j even In this section, we derive some lower bound on |c^ + -. | for methods inz,[k] which are AQ-stable. We first consider the case when k is even. A necessary condition for the method (bQ,b-| ^v-\) t0 ^ e ^Q-stable is that b > 0, b 1 >. b- > j = 2,3 k-1 (cf. theorem 3.3). Since, by normalization, b k = 1, a lower bound for lc k+1 l is 1 I c k+l I = " c k+l — k K+l K+l 2 K (|<+1) the inequality being strict if k > 4. The case when k = 2 has been dis- cussed in section 3.3.2, where it was shown that the lower bound is at- tained by the trapezoidal rule. A similar argument when applied to the case when k is odd will lead to the lower bound l c k+l' = " c k+l - ° with strict inequality when k > 3. Again, the lower bound is attained by the trapezoidal rule in the case when k = 1 (cf. section 3.3.1). It will be shown in the next section that for any c^ < 0, there exist methods in l[3] which are Ag-stable and have error constant c^. For the general case when k l 5, k odd, a lower bound for Ic^-jl can be derived from the fol- lowing lemma: 33 Lemma 3.1 Let k > 5, k odd. If (dq,d-|,.. . ,b k _-j) is an A Q -stable method in L[k], then ,. b 4 b k-l, b k . 1 (b 2 + r + ... + _) ><; n 9 1 k = 5 k > 7 Proof Suppose the contrary, i.e., n k-l b a, b k_1 v=2 v-1 2 v even = vvi_ j k = 5 k > 7 (3.5) U roots of From theorem 3.1, A n -stability implies that, for any qe(0,°°), the qr(z) k , s(z) + = E H.(q)z J J-0 where H k (q) = 1 a- H.(q) = b. + qJ- j = 0,l,...,k-l lie in the interior of the left half complex plane. Hence, it is necessary that the following inequality be satisfied for all q > (cf. Appendix A): H k . 1 (q)H 1 (q) - H Q (q) > In terms of powers of q, the left hand side can be written as 34 %z + [ !i!ki.(!o. bi ) ]q + (bk _ ibi . bo) Since b > for v >_ 2 (cf. theorem 3.3), the coefficient of q is negative from (3.5). In order that the above quadratic (in q) be positive for all q > 0, it is necessary that [ !!^l..(i. bl , ] 2 <2l ;Cb |t . 1 h 1 -V or 2aibo + £!!fci . 1,2 . 2( !i^i . 1 )( !° . 1, + 1 u 2 k 2 k 2 ' k + (— - b, - -) 2 - 2a, b. ,b, < '1 'i u k-ri (3.6) which implies that <7- b i-;;> 2 - 2a i b k-i b i< since the first three terms in (3.6) are nonnegative. However, from (1.9), i. 2 b k . 2 b 3 (k-iy (^ - b - -) 2 ' k " 3(k-2) ~ 12(k-2) b i 1 ( -t 9 4 -b k > 7 v' > 2a 1 b k _ 1 b 1 making use of the fact that, by theorem 3.2, s(z) has no root in the open right half plane, and hence (cf. Appendix A) 4b k . 2 b 3 > (k-n'b, We thus have a contradiction. Q.E.D. 35 Theorem 3.6 Let k > 5, k odd. If (b Q ,b, b k_"|) is an A Q -stable method in £[k], then 1 C k+1 I > k K-M 2 K 3|< Proof Suppose that k-1 b 1 E <_ — v=2 v+1 ~ 3k v even Since b v > for v > 2 (cf. theorem 3.3), b. , k-1 b 1 J^i< e ^< — k v=2 v+1 " 3k v even and thus 3b k-i < ] Therefore, by lemma 3.1, 1 k-1 b 1 k-1 b k-1 b _ < k , E -v. < . s _ ± < E _v 3k k_1 v=2 v-1 3 v=2 v-1 v=2 v+1 v even v even v even We arrive at a contradiction. Since bQ >_ (cf. theorem 3.3), 1 k-1 b v 1 c k+ il>-jr £ > -c — Q.E.D. K ' 2 K v=2 v+1 2 K 3k v even It should be remarked that the above lower bound can be strengthened. The important point is that |c k+ J is bounded away from zero 36 if the method is A.-stable. Moreover, the infimum of |c. ,,| for A n -stable ' k+1 ' 2 methods lies between l/(2 k 3k) and 2 k +k-1 . 3.5 HIGHER ORDER A(a). STABLE METHODS In this section, we demonstrate how, for given ae(0,TT/2), A(a)-stable methods in L[k] can be constructed, and by doing so, prove that A(a)-stable methods of arbitrary order k exist for any ae(0,Tr/2). The ac- curacy of such methods is also discussed. A statement about the accuracy of Aq-s table methods in L [3] is given as a corollary. Since the cases when k = 1,2 have been analyzed in section 3.3, we assume throughout that k > 3. Dill and Gear (1971), when searching for stiffly stable methods of order greater than 6 (cf. section 1.4), make an interesting observation that there is a greater likelihood of stiff stability for methods whose p(0 polynomial has multiple roots near 1. Consider the method represented by p(o = u-n k . o(o = te+sKc-n^Va+fi) The common factor (c-1) makes the method of order k automatically. The method "almost" satisfies the root condition and has the same stability re- gion as P(C) -C-l. cr(c) = (c + 6)/(l + 5) which is A-stable if 6e(-l,l]. Now by choosing k-1 roots of a(c) to be slightly less than 1 instead of 1 , we may be able to get an A(a)-stable method in L[k] with a close to tt/2. To this end, consider the two parameter family of methods parameterized by the following polynomial: s(z) = (z +d)(z + D) k " ] (3.7) 37 where < d < D Since b ,b-|,... ,b k are coefficients of the polynomial s(z), using elemen- tary algebra, we can show that, for each j = 0,1,..., k, k , k-1 d k-1 b,- = D K " J [( ) + -( )] 3 M D j = D k " J Bj (3.8) where B,- is the quantity inside the brackets and is bounded by k B.- 1 ( ) (3.9) J j Since the method (r,s) is of order k, the coefficients of the po- lynomial r(z) are given by (cf. (1.9)) k b a, = 2 z — £ j = 0,1 J V=j+1 V-j v-j odd and hence k-1 i k s (z) 1 \*- J L ""i^ fc L j=0 J m=l m m odd where k-m s m (z) = i b.,z J m v -_ m+j (3.10) are monic polynomials of degree k-m, 1 <_ m <_ k. Define the monic polynomi- als u(z) and u-,(z) of degrees k-1 and k-2, respectively, by u(z) = s(z)/(z + d) (3.11) 38 U] (z) = (u(z) - u(0))/z Since s-j(z) = u(z) + du-j(z), the rational function r(z)/s(z) for z = iw, ^(- O0 » 0O )» is then given by r(iw) 2 ui(iw) k s (iw) = [1 + d- + E — ] (3.12) s(iw) iw+d u(iw) m=3 mu(iw) m odd If we are able to choose d and D so that the quantity u^iw) k s m (iw) d + i u(iw) m=3 mu(iw) m odd is sufficiently small, for any we(-°° t °°), with respect to the constant 1, then r(iw) m 2 = for all we(-°°,°°) s(iw) iw+d in which case we may be able to get A(a)-stability for a close to tt/2. And that is the main idea behind the discussion that follows. We shall find a uniform upper bound in terms of d and D for the following "undesired" por- tion of the function r(iw)/s(iw) over -» < w = Dy < °°: Ui(iDy) k s m (iDy) d + Z u(iDy) m=3 mu(iDy) m odd Then by appropriately choosing d and D, it will be shown that we are able to get A( a )-stability for any given aG(0,ir/2). Note that 2/(iw+d), we(-°°,°°), is the locus of a circle in the complex plane centered at 1/d with radius 1/d. An expression for |u(iDy)| is easily obtained from (3.7) and (3.11) by substituting z = iDy and factoring out the factor D, 39 k-1 | u ( i Dy ) | = D k_1 (y 2 + 1) 2 (3.13) We next find an upper bound for |s (1Dy)|, <_ m < k. From (3.8), (3.9), and (3.10), |s m (iDy)| 2 * \7 D k - m -\ +j (iDy) J '| 2 < D 2(k - m) P k . m (y 2 ) j ** 2 2 where Pi c _ m (y ) is a monic polynomial of degree k-m in y , independent of the parameters d and D. Combining (3.13) and the above bound, we have, for any m > 1 , s m (iDy) M' u(iDy) D m_1 uniformly in ye(-°°,°°), in which M' is some positive constant. It follows that, if D >_ 1, . k s m (iDy) M' k 1 M I I -* 1 < -» E - 1 tt<- T (3.14) m=3 mu(iDy) D c m=3 mD D^ m odd m odd where M is a constant independent of d and D (e.g. M = M'k). By a similar argument, we can prove that there exists a constant K > such that u^iDy) < K u(iDy) ~D uniformly in y6(-°°, 00 ). Together with (3.14), we obtain the following bound for the "undesired" portion of r(iDy)/s(iDy) (cf. (3.12)) when D >_ 1: Ui(iDy) k s m (iDy) dDK + M \d- + i — 1< ^— (3.15) u(iDy) m=3 mu(iDy) D m odd We now proceed to show how, for any given ae(0,ir/2), the 40 parameters d and D can be chosen so that the method is A(a)-stable. Assume that oe(0,tt/2) be given. To get A(a)-stability t we want r(iDy) | Arg 1 <_ tt - a s(iDy) for any ye(-°°,°°) such that r(iDy) + (cf. (3.1)). Consider an arbitrary ye(-°°,oo). if we choose d and D _> 1 such that dDK + M 2~ < sin(- - a) (3.16) then r(iDy) 2 ih(iDy) k s (iDy) |Arg 1 < |Arg 1 + | Arg [1 + d-^ + z — ]| s(iDy) iDy+d u(iDy) m=3 mu(iDy) m odd t dDK + M <_ - + sin ( « — ) <_ tt - a 2 T Therefore, from (3.1), the method (r,s), if convergent, is A(a)-stable pro- vided (3.16) holds. It remains to show that the polynomial r(z) , when (3.16) is satisfied, has roots in the left half plane. We shall need the following theorem (Levinson and Redheffer, 1970, p. 218): Theorem 3.7 (Rouche's theorem) Let f(c) and g(c) be analytic in a simply connected domain con- taining a Jordan contour J. Let |f(c)| > |g(c)| on J. Then f(c) and f(0 + g(0 have the same number of zeros inside J. Since k s, U) r(z) = 2[u(z) + d Ul (z) + z -5 ] m=3 m m odd 41 from (3.15) and (3.16), r(iw) - u(iw)| < |u(iw)| (3.17) " 2 for any we(-°°,°°). Applying theorem 3.7 with 5+1 1 c+1 ?+l f(c) = u(— ), g(0 = -r(— ) - u(— ) C-l 2 c-1 C-l and J being the unit circle, we conclude, from (3.17), that the polynomials u(z), r(z)/2 have the same number of zeros in the interior of the left half z-plane. Hence all the roots of r(z) lie in the open left half plane, and the corresponding method (r,s) is convergent (cf. section 1.3) -- in fact, strongly stable. We can hence obtain, for any ae(0,Tr/2), methods in L[k] which are strongly stable, A^-stable, and A(a)-stable by choosing the parameters d and D > 1 such that (3.16) is satisfied (Slight modification needed if D < 1). It should be remarked that d can be made arbitrarily small, and in fact zero, though the method will not be A^-stable when d = 0. From (3.8), the magnitude of the error constant c^ + -, is given by 1 k k-1 D k " j k-1 k-1 D k " J_1 2 k J-2 j-l j+l j=0 j j+l j even j even If we choose d such that d = 0(D~ ), then from (3.16), for A(a)-stability, D = 0(- ) as a + - o- - a 2 k-2 Hence, the magnitude of the error constant, being of order 0(D ) for lar- 42 ge D, is of order i k - 2 •TT 2 0([- ] ') as a + 1 -a 2 The results in this section can be summarized by the following theorem: Theorem 3.8 For arbitrary k and any og(0,tt/2), there exist A(a)-stable k-step methods of order k which are A^-stable and strongly stable. Corollary 3.8.1 The method corresponding to s(z) = z(z + D) k_1 (3.19) is A -stable if D 2 > 2 k " 3 + 1. Proof To prove the corollary, all we have to do is to obtain value for M' (cf. (3.14)) in terms of k when k >_ 3 (the case when k <_ 2 is trivial -- see section 3.3). We shall make use of the following identities, the veri- fication of which is straightforward: k-1 k-1 . * ^ (k-lh I ( ) = 2 + 2 c cos (3.20a) j=0 j 4 j=0 mod 4 k-1 k-1 v o ^ (k-l)Tr Y. ( ) = 2 J + 2 c sin (3.20b) j=l J 4 j=l mod 4 43 k-1 k-1 k o t (k-1 )tt z ( ) = 2 - 2 * cos (3.20c) j=2 j 4 j=2 mod 4 k-1 k-1 . , ^ (k-1 )tt z ( ) = 2 k ~ J - 2 * sin (3.20d) j=3 j 4 j=3 mod 4 Since for any <_ j <_ k-1 , k-1 |y| J i (y 2 + 1) 2 for all ye(-oo,c») we have, for m > 3, m odd, from (3.8), (3.10), (3.13), and (3.20), ^(m-D^m^V r , k " m , k " ] , k " m , M U1 2 D ; | 1 < [max{ z ( ), z ( )}] + u(iDy) j=0 m+j-1 j=2 m+j-1 j=0 mod 4 j=2 mod 4 k-m k-1 k-m k-1 ? + [max{ z ( ), z ( )}] j=l m+j-1 j=3 m+j-1 j=l mod 4 j=3 mod 4 k - 4 k4 < 2[2 k " 3 + 2 2 ] 2 < [2 2 ] 2 Therefore, if D _> 1 , k-^ k- 3 k s (iDy) 2 2 k 1 2 2 2 k " 3 m v J | m=3 mu(iDy) ~ D 2 m=3 mD m " 3 " 3(D 2 -1) D 2 -l m odd m odd uniformly in ye(-°°,°°). We note that the above upper bound is different from that in (3.14). In fact (3.14) is satisfied if we choose M = k2 k " 3/2 . Using a similar argument as in deriving (3.16), we can show that the method represented by (3.19) is A -stable if 2 k-3 D -1 or 44 ? k-3 D^ 1 2 K J + 1 Q.E.D. Remark From (3.18), the magnitude of the error constant of the method 2 k-3 represented by (3.19) when D = 2 + 1 is of order 2 0((k-l)2 k /2 ) as k ■> » Corollary 3.8.2 For any ae(o,7r/2) and any C > 0, there exist A(a)-stable 3-step methods of order 3 with error constant -C. Proof Assume that ae(0,fT/2) and C > are fixed. Consider the class of third order methods represented by s(z) = z(z 2 + 24Cz + R 2 ) where R ^max{l,12C}. The magnitude of the error constant of the method is C. It is obvious that 1 r(z) = 2[u(z) + -] 3 where, from (3.11 ), u(z) = z 2 + 24Cz + R 2 Hence the root condition is satisfied. The locus r(iw)/s(iw) can then be expressed as r(iw) 2 1 = -[1 + ] s(iw) iw 3u(iw) We first find a uniform upper bound for l/(3u(iw)) over -°° < w = Ry < oo. If R > 12,/ZC, then 45 24C — 12C — |u(iRy)| ■ R 2 [(l-y 2 ) 2 + ( ) 2 y 2 ] 2 > 24RC[1 - ( ) 2 ] 2 > 12»^RC R R Thus 1 1 3u(iRy) ~R(36^TC) uniformly in ye(-°°,°°). Using a similar argument as that in (3.16), the method is A(a)-stable if R satisfies 1 < sin(- - a) R(36^C) 2 in addition to the assumption that R >_ max{l ,12»/2~C}. Q.E.D. Remark 2 Note that the roots of r(z) have magnitude R +1/3. This implies that as a -> tt/2, the polynomial p(c)/(c-l) has roots close to 1. The corollary simply demonstrates the inadequacy of the error constant as a measure of accuracy for k = 3. 3.6 CONCLUDING REMARKS As remarked in section 1.4, an ODE solver could have for each or- der k a family of formulas parameterized by a. We have to know what a we want before we can select the right formula. While the order k is chosen automatically, the angle a would be either supplied by the user or also determined by the code. As an illustration, for specified values of k and a, we could choose among the class of methods of the form s(z) = (z +d)(z + D) k - ] the method whose parameters d and D satisfy (3.16) (the values for M and K can be found using a similar argument as in the proof of corollary 3.8.1). 46 The method will then be of order k and A(a)-stable. We note that the above class of methods can be parameterized by the two parameters d and D, or by the parameter D alone, in which case d can be a fixed constant or a func- tion of D (d has to be of order o(D) so that (3.16) can be satisfied when a ■*■ tt/2). However, the selection of the formulas for various a might best be done through some table residing in the code. The simplest such ar- rangement would limit the order depending on a. 47 CHAPTER IV MINIMAX 4.1 INTRODUCTION In the next chapter, we shall be concerned with investigating how large the angle of absolute stability can be among methods having error constants of the same magnitude. Since an analytical solution is difficult to obtain, we shall solve it using numerical methods. If we restrict our attention to A OT -stable methods only, then from (3.2), a method (r,s) is A(a)-stable if Re r(iw)/s(iw) sup - <_ cot a -« 3 (x) = max f(x,y) (4.2) yeY over the closed (not necessarily bounded) set X. The set X is called the feasible region or constraint set and any xeX is a feasible point. The approach that is used in solving most constrained optimiza- tion problem goes as follows: first, Lagrange multiplier theorems descri- bing the necessary conditions that an optimal solution must satisfy are formulated, then an iterative search is devised to locate points which sa- tisfy the necessary conditions. Those points are usually called stationary points. To this end, we require that the functions f(x,y) and g(x,u) be continuously Frechet differentiate (Luenberger, 1969, p. 172) with respect to x on X'xY and R n xU, respectively. It can then be proved that the func- tion 4> is Gateaux differentiate (Luenberger, 1969, p. 171) at each xeX', the Gateaux differential of at x in the direction d being given by (Dem'yanov and Malozemov, 1974, p. 188) 9cJ> 3f — (x) = max (— (x,y),d) (4.3) 3d y^R(x) 3x where ( , ) is the scalar product of two vectors and R(x) = (ysY | (x) = f(x,y)} (4.4) A stronger statement is in p. 191 of the same reference: 49 af (x+hd) = 4>(x) + h max ( — (x,y),d) + o(h;d) (4.5) yeR(x) 3x where o(h;d)/h -> as h -> uniformly in d, ||d|| = 1 (| is the Eu- clidean norm). 4.3 A FEASIBLE DESCENT ALGORITHM As pointed out in the previous section, we first state a theorem describing stationary points of on X (cf. (4.1) and (4.2)) and then pro- pose a feasible descent algorithm to find a stationary point. The fol- lowing notation will be needed, in which xeR n , yeR m , uelR p : Q(x) = {usU | g(x,u) = 0} 3f H(x) = {— (x,y) | yeR(x)} 9X H'(x) = {— (x,u) I ueQ(x)} 3X L(x) = conv [H(x) u H'(x)] s s COnv S = { l a.Z. Z.GS, a.>0, l0, l(x) to achieve its minimum on X at a point x*eX is that conv H(x*) n cone [- H'(x*)] 4 An equivalent necessary condition, but computationally easier to test, can be derived using the same idea as in Dem'yemov and Malozemov (1974, pp. 146-148 in which U is assumed to be a finite set only): Corollary 4.1 The necessary condition in theorem 4.1 is equivalent to 0eL(x*) (4.6) Definition A point x*eX which satisfies (4.6) is called a stationary point of on X. We proceed to describe a feasible descent algorithm for finding stationary points of <\> on X. For any feasible point x,.a nonzero vector deR n is said to be a feasible direction on X at x if there exists a scalar h > such that x+hdeX for all he[0,h] (Zoutendijk, 1966). A necessary condition for a nonzero vector d to be a feasible direction at xeX can be shown to be 39 (— (x,u),d) <^ for all ueQ(x) (4.7) 8X We note that some authors call da feasible direction if (4.7) holds. It 51 can be proved that, when the assumptions of theorem 4.1 are satisfied, the closure of the set of feasible directions at x is equal to the set of dejR n satisfying (4.7) (Dem'yanov and Malozemov, 1974, p. 199). We call a nonzero vector d a feasible descent direction for on X at x if there exists a po- sitive scalar h* such that, for any he(0,h*], x+hdeX and (x+hd) < (x) (Zoutendijk, 1966). The feasible descent algorithm works as follows. A feasible starting point x Q is given. Suppose that we have obtained points Xg.x-. ,... ,x eX. If x is a stationary point, the algorithm terminates. Otherwise, a new point x + -,eX is determined such that (x +1 ) < (x ), or it is concluded that no stationary point exists. The latter may occur sin- ce X is not bounded. The determination of x +1 consists of two parts: a feasible descent direction d at x is found, then a step is made in the direction d by choosing a step size h so that the new point x +1 = x +h d is in X and satisfies 4>(x , 1 ) < (x ). If the two conditions can V V V V+ I V be satisfied by any h > 0, or if |x | becomes unbounded as v increases, then we conclude that no stationary point exists. The feasible descent directions can be determined using the bin- ding functions and the active constraints. Since the function <(>(x) is not Fr6chet differentiate, such strategy could result in convergence to a non- stationary point (example in Dem'yanov and Malozemov, 1974, pp. 76-82 where X = R ). Such a phenomenon is known as zigzagging. It occurs when the sequence of steps becomes progressively very small not because a stationary point is in the vicinity but because new binding functions or new active constraints are encountered. One remedy is to take into consideration not just the binding functions and the active constraints but also the "nearly" 52 binding functions and the "nearly" active constraints. This also serves to prevent binding functions and active constraints from "disappearing" due to roundoff errors. Let e,y be nonnegative numbers and define the sets R e (x) = (yeY | *(x) - f(x,y) < e} (4.8) Q y (x) = (ueU | - u ig(x,u) < 0} 9f H e (x) = {— (x,y) | yeR (x)} 3X 3g H'(x) = {— (x.u) I ueQ (x)} M sx y L £y (x) = conv [H £ (x) u IT(x)] (4.9) Note that R (x) = R(x), Q Q (x) = Q(x), H Q (x) = H(x), Hj(x) = H'(x), and Lqq(x) = L(x). Moreover, the sets H (x) and H'(x) are compact and hence L (x) is a compact convex set. Let z (x) be the point of L (x) nearest ey ey ey the origin, i.e. , I |z (x) I I = min I |z I £y zgL (x) ey x ' If z (x ) = 0, or equivalently, OgL (x ), then we call x an (e,y)- quasistationary point of 4> on X, in which case e and y are reduced and the search continues with the new e and y. Note that (O,0)-quasistationarity is equivalent to stationarity. If z (x ) \ 0, define z (x ) d = - E1LJ! (4.10) ■12 (xjll :y x \>' I I ei The vector d is a feasible descent direction as shown by the following 1 emma : 53 Lemma 4.1 For any e > and u > 0, if x v is not (e,y)-quasistationary, then d as defined in (4.10) is a feasible descent direction for on X at x . Furthermore, Proof 8 — (x ) < - I |z (x ) 3d V £y V Since z (x ) is the point of L (x ) nearest the origin, it can be shown that (Dem'yanov a nd Malozemov, 1974, p. 252) IU ep (x v )|| 2 i(z,z Evl (x v )) for all zeLjxJ Thus, from (4.10), we have (z,d ) < - | |z (x )| | < for all zeL (x ,) (4.11) \ » v / _. II £ yV V / I I £ y\ y/ \ / In particular, 3f (— (x v ,y),d v ) i- ||z (x v )|| for all yeR(x v ) 3X (— (x ,u),d ) l- llz (x )|| for all ueQ(x ) 3X Using (4.3) and (4.5), and a similar argument for the function max g(x,u) ueu * it is then obvious that d v is a feasible descent direction. Q.E.D. After determining the feasible descent direction d , we next want to take a step in that direction. The stepsize h can be found by the fol- lowing step selection rule, which is a modification of the Armijo rule (Po- lak, 1971, p. 36): 54 Let y > 0, 3^(0,1), and a€(0,l) be fixed constants. The values y = 1, Be(o.5,0.8), and a = 0.5 are recommended by Polak (1971, p. 36). We start with an initial value for h given by h v = maxW.-jj — > If h v leads to an infeasible point, it is reduced by the factor 3 until the feasibility constraint is satisfied. Note that since d is a feasible direction, at most a finite number of reductions is needed. There are then two possibilities. If 4>(x +h d ) < <}>(x ) - ah I |z (x )| I (4.12) then we define x , , = x + h d . Otherwise, the stepsize h is reduced successively by the factor 3 until the above inequality is satisfied. The h corresponding to the smallest (x +h d ) obtained in this iteration is the stepsize used in defining x .-, , i.e. , x .-, = x + h d . Since d is a r 3 v+ 1 ' v+lvvv V feasible descent direction, the number of reductions must be finite. In case llx v +h v d v | is large, it is conceivable that no stationary point ex- ists, hence the algorithm terminates. The constant y ensures that the initial stepsize is not too small, whereas the choice h^.-j/B enables us to try a possible larger step. We remark that the original Armijo rule uses y as the initial trial and chooses as the stepsize the first h such that x +h d eX and V V V V <|>(x +h d ) < (x ) + ah (x ) V v v — V V- . V V Since (x) is not Frechet differentiable, we replace 3(x )/3d by - |z (x )| which is, from lemma 4.1, greater than the former quantity. 55 Note that we can choose as the stepsize the first h v such that x v +h v d v eX and (4.12) is satisfied. The convergence theorem proved in the next sec- tion still holds, and in fact holds for any line minimization rule which gives a better point, i.e., smaller cf>(x + ,), than the rule just mentioned. Since it is more desirable to have a larger decrease along the descent direction, we choose as the stepsize the one corresponding to the smallest 1 }, iv^\ kill, and (pjj k >_ 1 } be the strictly decreasing sequences of positive numbers used by the M-algorithm. We use the superscript k, i.e., x , to denote the point re- turned by the (e,y,p)-algorithm with e = e. , y = y. , p = p.. In other words, each x , k >_ 1 , satisfies l|ze k p k (x k )||

_ 1 } 56 to a stationary point is discussed in the next section. We finally remark that an approximation for the vector z _..(x) can be obtained as follows. Discretize the sets R £ (x) and Q u (x) and define the matrix D whose columns consist of vectors in H £ (x) and H'(x) evaluated at the discretized points (cf. (4.9)). Then z (x) can be approximated by Da* where a* is the optimal solution of the following convex quadratic program- ming problem: • • • T n T n minimize a D Da subject to a. ^ 0, l' = l,...,q a = [a-j^ a ] a,+ap+. . .+a = 1 The above problem can be solved by a number of methods, e.g., the "first symmetric variant of the simplex method" (Van De Panne, 1975, pp. 270-280) which is shown to converge to the optimal solution in a finite number of iterations. 4.4 CONVERGENCE In this section, we prove that the M-algorithm described in the previous section converges to a stationary point if it exists. The proof is in two stages: the convergence of the (e,y,p)-algorithm for positive e,y, p and then the convergence of the M-algorithm to a stationary point as e. -»■ 0, ijl + 0, and p^ ■*■ 0. The existence of stationary points is assured if the sequence (x ± v <_ v^, k >. 1 ) generated by the M-algorithm is bounded or the set X(x°) = {x£X | +(x) l*(x°)} is bounded, where x ex is a given starting point for the M-algorithm. For 57 simplicity, we assume that the latter holds. Lemma 4.2 For any e > 0, there exists 6 = 6(e) > such that R(x') c R £ (x) for all | |x'-x| | < 6, x',xex(x°). Proof The assertion follows directly from the definition of R(x) and R (x) (cf. (4.4) and (4.8)) and from the uniform continuity of (x) on X(x°) and f(x,y) on X(x°)xY. Q.E.D. Lemma 4.3 For any e > 0, n > 0, there exists h = h (e,n) > such that (x+hd) - <{>(x) 3f max (— (x,y),d) < n h y^R e (x) 9x Ox ___, _^ m n . _, lL ., lLJ ^, Proof for all he(0,h ], xeX(x°), and deK n such that x+hd£X(x°) and ||d|| = 1. From the uniform continuity of 3f/9x on X(x )xY, there exists po- sitive h' = h'(n) such that 9f f(x+hd,y) = f(x,y) + h(— (x,y),d) + hA(x,y,hd) 9x where |A(x,y,hd)| < n for any x£X(x°), yey, he(o,h'], and deK n such that x+hdeX(x°) and | |d| | = 1. By lemma 4.2, there exists positive h" = h"(e) such that R(x+hd) £R £ (x) for any xeX(x°), he(0,h"], and deR n such that x+hdeX(x°) and ||d|| = 1. 58 Let h^ = minlh'.h"}. Then for any xeX(x°), he(0,h ], and deR n such that x+hdeX(x°) and ||d|| = 1, <}>(x+hd) = max f(x+hd,y) <_ max f(x+hd,y) yeR(x+hd) y£R (x) 3f < max f(x,y) + h max ( — (x,y),d) + hn yeR £ (x) y e R £ (x) 3x 3f = and n > 0, there exists h = h (y , n ) > such that \p (x+hd) - i/j (x) 8g max ( — (x,u),d) < n h ueO (x) 3x for all he(0,h ], xeX(x°), and deff n such that x+hdeX(x°) and ||d|| = 1, where ip(x) = max g(x,u) ueU Proof Similar to that of lemma 4.3. Theorem 4.2 For any e > and y > 0, if (x v >_ 0} is the sequence of points generated by the (e,y,0)-algorithm with initial value x n eX(x ), then |z (x )| I + as v ■> », v > 1 ' ey x v ' ' ' ' — Proof Let n > be given. Define (cf. lemma 4.3, 4.4) 59 h* = mirUh^e.O-aJnJ.h^y.rO.Y} Since U(x ) \j > 0} is a monotone sequence bounded below by inf <}>(x) = inf n <}>(x) > - °° xeX xex(x ) there exists a positive integer N = N(e,y,n) such that *(x ,) - (x l+1 ) < agh*n for all v >_ N 'v+1 Let h be the stepsize last tried by the step selection rule (not neces- v sarily the stepsize chosen). Then x +h d eX(x ) and (x ) - (x .,) > <|>(x ) - (x +h d ) > ah llz (x )| TX v y YX v+l'— TV v' rv v v v 7 — v 1 ' ey ■ v ' ' Consider v > N. If h > 3h*, then — v — *(x J - *(x .,) a3h*n z (x ) < v - ^- < (x +[h /3]d ) - *(x ) > - a z (x ) i i £1J \ v / h ,/e v Since, from (4.11), ( — (x,,y),d ) < - |z (x )| for all yeR (x ) 3X using lemma 4.3, we obtain z (x )| < n ey v v" 1 Case 2. > *(x v ) 60 or ♦ (x v +[fi v /B]d v ) - »(x y ) o h v /B Using (4.11) and lemma 4.4, we have z (x )| < n Q.E.D. ey v ' ' One immediate corollary of the above theorem is that the (e,y,p)-algorithm terminates after a finite number of iterations if e ,y iP are positive. Note that the above theorem holds if the new point is de- fined by n , i.e. , x , , = x + h d , as remarked in the previous section. V V+l V V V To prove that the Pl-algorithm yields stationary points of on X, k we have to show that for small enough e. ,y. ,p. , if x is close to an k (e y )-quasistationary point, then x is also close to a stationary point. We shall need the following lemmata: Lemma 4.5 For any xeX and any n > 0, there exist positive 6 = 6(n,x), e = e(n,x), and y = y(n,x) such that the following hold for any ||x'-x| < 6, x'eX: (i) For any ee[0,F], eyery element of R (x 1 ) is within an n-neighborhood of an element of R(x). (ii) For any ye[0,TT], e^ery element of Q (x 1 ) is within an n-neighborhood of an element of Q(x). Proof Suppose there exist a positive number n', a positive sequence ie^ i 2l 0) converging to zero, and a sequence of points (x -eX i >. 0) con- 61 verging to x such that for any i > 0, there exists y.eR (x. ) which is not l £ i 1 in an n '-neighborhood of any veR(x), i.e., ||y i - v| | > n' for all vGR(x) (4.13) Since y.eY for all i * 0, there is a subsequence {y^ iel>, I c {0,1,...}, converging to y£Y. For each lei, (x.j) - fU^) < e n . By continuity, as i -> °°, iei, 4>(x) - f(x,y) < which implies that yeR(x). It follows that for sufficiently large i, there exists yeR(x) such that ||y--y|| ll^oll/ 2 Proof From the hypotheses, it can be proved that every element of conv G' is within an |z Q | |/2-neighborhood of an element of conv G, i.e., for any z'Econv G', there exists zeconv G such that ||z' - z|| < ||z ||/2 62 It follows immediately that ll z 'll > ll z M- M Z M /2 = M z oH /2 Q - E - D - Lemma 4.7 For any xeX, there exist positive 6 = 6(x), e = e(x), and y = y~(x) such that for all ||x'-x| < <5, x'eX, ee[0,eTl, ye[0,vT], !|z Ep (x')|| > ||z 00 (x)||/2 Proof The proof for the case when ||zqq(x)| | = is trivial. Suppose that ||zqq(x)|| f 0, i.e., x is not a stationary point. From ths uniform continuity of the functions 3f/3x and 3g/3x in yeY and ueU, respectively, we can find n = n(x) >0 such that the following hold for any |x'-x| | < n, x'eX: 3f 3f l|z nn (x)|| II— (x.y) (x',y')|| < — M for all y.y'eY, I |y-y' 1 1 < n 3x 3x 2/n 39 39 ll z nn( x Hl ||_(x,u) - _(x\u')|| < ^ for all u.u'eU, | |u-u' | | < n 3x ax 2>/rr Applying lemma 4.5, we obtain positive numbers 6' = 6'(x), e = e(x), and y = y(x) such that the hypotheses of lemma 4.6 are satisfied with 3f 3g G = {— (x,y) | yeR(x)} u {— (x,u) | ueQ(x)} 3X 3X 3f 3g G' = {— (x'.y 1 ) | y'eR (x 1 )} u {— (x'.u 1 ) | u'gQ (x')> 3X 3X y for any x'eX, ||x'-x| < 6 = min{n,_ 1 } generated by the M-algorithm is a stationary point of on X. Proof k Let x* be a limit point of {x k >_ 1}. Then there is a subse- k quence {x keO, I c {0,1,...}, converging to x*. Using lemma 4.7, there exists a positive integer K such that for all k >_ K, kel, Since p.+0as k + °°, so must ||z 00 (x*)|| ■* Q.E.D. Theorem 4.3 does not guarantee the convergence of the sequence (x k >_ 1). If we further assume that the function has at most a finite number of stationary points, it could be shown that the generated sequence will converge to a stationary point of on X. 4.5 MODIFICATIONS Since X is described by a nonlinear inequality, the function g(x,u) can be defined arbitrarily in the sense that the same region is ob- tained if we replace g(x,u) by g(x,u) times a positive function. Subse- quently the choice of feasible descent directions is also arbitrary depen- ding on how g(x,u) is defined. Such an arbitrariness can also be seen from theorem 4.1. The vectors in H'(x*) can be adjusted by any positive number and the necessary condition for x* to be a global minimum still holds with H' (x*) replaced by 64 3g H'(x*;t) = {t(x*,u)— (x*,u) | ueQ(x*)} 3X where t(x*,u) > for all ueQ(x*). The convergence property of the M- algorithm will not be altered if we require that the numbers t(x,u) be po- sitive and be bounded above by a constant. One particular choice for t(x,u) is motivated by the following consideration in the two dimensional space: Suppose that in choosing a feasible descent direction we approxi- mate the constraint set by a disc and use a linear approximation for the objective function. The resulting problem is to minimize f(x) = c T x c.xeR 2 over the constraint set defined by X = {xett 2 | g(x) = (x - a) T (x - a) - r 2 £ 0}, aett 2 Note that with appropriate scaling of the decision variables, a convex con- straint set can be approximated by a disc. It will be seen that the radius of the disc is immaterial in the determination of feasible descent direc- tion. Suppose that we start with an initial point x lying on the boun- dary of X (see figure 4.1). By the linearity of the objective function, the optimal solution x* must lie on the boundary. Thus the vector x*-x gives a feasible descent direction which leads to the optimal point in one step. Geometrically, x*-x is the vector which bisects the angle formed by the vectors 3f(x )/ax and 3g(x )/3x. Using the Kuhn-Tucker theorem (Luen- berger, 1969, p. 249), we can show that x*-x is given by 65 f(x) = f(x) = f(x u ) Figure 4.1. A Two Dimensional Case ifV) £<*°> C* - X° = - 3f o. -{— (x u ) + 3f(x )/ax 1 1 3g 0, I c | 1 3g(x )/ ax | 3x x u )> Note that the direction of x*-x is independent of r, the radius of the disc. Thus it seems like it may be advantageous to scale 3g(x)/3x so that it has the same norm as the vector 3f(x)/3x. Based on this observation, we could pick t(x,u) so that the vectors t(x,u)3g(x,u)/3x, ueQ (x), have norm equal to the smallest norm of the vectors 3f(x,y)/3x, yeR (x). The above example also indicates that an appropriate scaling of the variable x is desirable. One common way is to scale x by a linear transformation matrix T, i.e., x = U 66 Using linear algebra, we can show that the above scaling is equivalent to the algorithm in which the next point is defined by x xl = x + h TT T d V+l V V V instead of x ,, = x + h d . We shall not go into details. V+l V V V 4.6 CONCLUDING REMARKS We have described an algorithm for finding only stationary points of (J) on X. A stationary point need not be a local minimum, it can also be a local maximum or a local saddle point. The first possibility can only occur at the starting point x if there exists a yeR(x ) such that 3f(x ,y)/3x = 0. The second possibility can be tested by exploring dif- ferent directions to see if can be further minimized. Some sufficient conditions for a point x* to be a local minimum are suggested by Dem'yanov and Malozemov (1974, pp. 125-126). Note that if the convexity or the Slater assumption (cf. theorem 4.1) fails to hold, the M-algorithm can still be used to locate local minima of on X provided that the above tests are performed on the generated limit points. If the function $ is not convex, there may be more than one local minimum, and in the best case, we shall only find one of them. A partial remedy is to apply the M-algorithm several times using different starting points. It should also be remarked that the formalism is general enough to include cases when X is defined by more than one inequality constraint. However, because of the Slater condition, problems involving equality con- straints can only be handled if the equality constraints can be solved to reduce the original problem to a new problem of smaller dimension having 67 only inequality constraints. 68 CHAPTER V ACCURACY VS A(a)-STABILITY 5.1 DESCRIPTION OF THE PROBLEM In this chapter, we are concerned with discovering the upper bound on the angle of absolute stability (definition in section 1.2) for methods in L[k] having a given error constant. For each A > 0, define A*(a) to be the supremum of all a for which there exists an A(a)-stable method having error constant -A , and be zero if no A(0)-stable method with that error constant exists. Note that the problem independent part of the leading term of the asymptotic error expansion (2.1) of such methods has magnitude (hA) . Possible values for A*(A) are shown in figure 5.1. Our A*(a) Figure 5.1. Possible values for A*(a) objective is to find A*(a) for A > 0. We expect that in all probability A*(a) is strictly increasing on some semi-infinite interval in (0,°°). From sections 3.3 and 3.5, the A*(a) values for A > when k = 1,2,3 are as shown in figures 5.2, 5.3, and 5.4, respectively. An analytical solution for the general case when k ^4 is diffi- cult. We thus find the A*(a) values numerically. A mathematical formula- A*(a) Figure 5.2. A*(a) for k = 1 69 A*(a) 1/2/3 Figure 5.3. A*(a) for k = 2 A*(a) * A Figure 5.4. A*(a) for k = 3 70 tion is discussed in the next section. We remark that the results from section 3.5 provide a "lower bound" for A*(a) over some semi-infinite interval in (0,°°). Moreover, A*(a) approaches the line a = tt/2 as A -> °° at least as fast as ( A -2k/(k-2)^ It alsQ follows from section 3.4 that A*(a) = for < a <_ l/[2(3k) 1A ]. 5.2 MINIMAX FORMULATION OF THE PROBLEM A mathematical formulation of the problem described in the pre- vious section is suggested by condition (3.2). To use that, we consider only methods which are A -stable. For reasons that will be apparent later, oo we further restrict our attention to methods that satisfy the strict root condition. Since convergent methods satisfy the root condition (cf. sec- tion 1.3), and from theorem 3.5, A(0)-stable methods are almost A^-stable save for possible simple roots on the imaginary axis, the above restric- tions lead to only a "slightly smaller" feasible region. We note that there exist methods which are not stable but whose locus r(iw)/s(iw), w^- 00 , 00 ). is outside the wedge S for some a > 0, e.g., the 5-step method represented by (b ,b 1 ,b 2 ,b 3 ,b 4 ) = (0,15, .14, 3, .007) From (3.2), an A^ -stable method satisfying the strict root condi- tion is A(a)-stable if and only if Re r(iw)/s(iw) - 1 cot a for all w> (5.1) | Irn r(iw)/s(iw)| We need only consider the semi-infinite interval we[0,°°) because the com- plex conjugate of r(iw)/s(iw) is r(-iw)/s(-iw) so that the rational func- 71 tion on the left is symmetric about w = 0. To indicate explicitly the fact that both r(z) and s(z) are actually functions of b = (b~,b..,...,b. ,) and z, we write them as r(b,z) and s(b,z), respectively. Since |s(b,1w)| > for any w 1 0, (5.1) is then equivalent to Re r(b,iw)s(b,-iw) sup - < cot a (5.2) wlO |lm r(b,iw)s(b,-iw) | The restriction that the method satisfies the strict root condition ensures that the numerator and denominator of the rational function on the left are not zero simultaneously; so does the A^-stability assumption. Suppose that a value of A (> 0) is given. We then solve for the infimum of the quantity on the left of the inequality (5.2) over the family of methods in L[k] which are A^-stable, satisfy the strict root condition, i/ and have error constant of magnitude a . Hopefully, the arccotangent of the infimum is the value of A*(a). We proceed to put the problem into the form described in section 4.2. As remarked in section 4.6, the equality constraint 1 k b k 2 k j=0 j+1 j even can be used to reduce the dimension of the problem by one. One way is to use b-,,b 2 ,... ,b. _, as the decision variables and set k } k b i b o = A " T E ~ U 2 K j=2 j+1 j even We denote the decision variables by x = (x,,x,,,... »x. _•,) where x. = b., j = 1,2 k-1 , and write r(x,z), s(x,z) for r(b,z), s(b,z), respectively, implying that the equality constraint has been eliminated. 72 k-1 The feasible region should be the set of xeB at which the po- lynomials (in z) r(x,z) and s(x,z) have roots only in the open left half plane. We could use the Hurwitz criterion (cf. Appendix A) to obtain a system of nonlinear inequalities which describes the feasible region. A simpler way is described as follows: Consider the region X = {xeB k " 1 | g(x,w,j) < for all we[0,»), j=l,2} (5.3) where g(x,w,l) = - |r(x,iw) | 2 g(x,w,2) = - |s(x,iw) | 2 It is obvious that both g(x,w,j) and 3g(x,w,j)/3x are continuous in B x[0,°°)x{l,2}. The set X contains all x such that r(x,z) and s(x,z) have no roots on the imaginary axis (in the complex z-space). Thus X con- sists of a number of disconnected open sets each of which is either entire- ly feasible or not feasible. If the initial trial x lies in a feasible component, so will all subsequent trials if the step selection rule does not take too large a step landing into an infeasible component. To prevent that, we only accept those points x such that r(x,z) and s(x,z) have roots all in the open left half z-plane and reject all other points. In other words, the feasible region is implicitly described by X and the above test, but only the function g(x,w,j) is needed in the determination of feasible descent directions. The objective function for the minimax formulation is Re r(x,iw)s(x,-iw) f(x,w) = |Im r(x,iw)s(x,-iw) | whose value at (x,w) is equal to - cot | Arg r(x,iw)/s(x,iw)| . Although the 73 negative part of f(x,w) is unbounded, it is never needed in determining the maximum of f(x,w) over we[0,°°). From the fact that as w •*■ °°, f(x,w) -> g(x,w,j) -> - » j - 1,2 it is not difficult to show that the results of Chapter IV still apply so long as the positive part of f(x,w) is bounded, which is true for those x representing methods which are A(0)-stable. Hence, if we can find an x eX which represents an A(0)-stable A^-stable method satisfying the strict root condition, then we can use the M-algorithm to find the infimum of the func- tion (x) = max f(x,w) we[0,°°) over the set X(x°) = {xeX | <|>(x) < <(>(x )} together with the test on r(x,z) and s(x,z) mentioned earlier. The M- algorithm will then either conclude that no stationary point exists or give a point x* such that the value of (x*) is close to the value of a local infimum of $ on X(x ). Note that since X(x ) is not closed, the infimum may not be attained by any point in X(x ). The initial trial x can be found by first finding an A(0)-stable method and then by varying some of the parameters b ,b,,...,b. _, such that the magnitude of the error constant is A without violating the feasibility of the new point in the x-space. 5.3 NUMERICAL RESULTS The M-algorithm described in section 4.3 is implemented to run on the Cyber 175 using double precision arithmetic. The numerical results 74 will be given after a brief discussion on the Fortran program. The evaluation of <}>(x) is in two steps. First, we test if the method represented by x is A(0)-stable. From (3.2), the method, being A^-stable (see constraints given in the previous section), is A(0)-stable if and only if the following condition is not satisfied by any positive real w of Im r(x,iw)s(x,-iw) = 0: Re r(x,iw)s(x,-iw) < The test involves finding the positive roots of an even polynomial of de- gree 2(k-l ) since r(x,iw)s(x,-iw) = z (-1 ) J C(2j)w 2j + iw E (-1 ) J+1 C(2j+l )w 2j j=0 j=C = u(x,w) + iv(x,w) where C(J) = E (-D V a b. for j = 0,l,...,2k-l v=0 v J_v Note that from (1.9), it can be proved that C(j) = for j >_ k, j even. If the method represented by x is A(0)-stable, then (x) is found by locating all the local maxima of f(x,w) on <_ w < °°. Again, we have to solve for the positive roots of an even polynomial of degree at most 3(k-l) since it can be shown that for w > such that Im r(x,iw)s(x,-iw) f 0, 3f v(x,w)[3u(x,w)/8w] - u(x,w)[3v(x,w)/aw] — (x,w) = - 3W |v(x,w) |v(x,w) where the numerator simplifies to L3(k-1)/2J v E (-l) v z (2v-4j+l)C(2j)C( v=0 j=0 The set of local maxima contains R(x) as a subset. (x.w.l >l) J A(2j)w 2j g(x,w,2) z (-D J B(2j)w 2j 75 Similarly, the set Q(x) is contained in the set of local maxima over w^[0,°°) of the functions g(x,w,j), j = 1,2, given by (cf. (5.3)) k-1 E J-o k z j=0 where for j = 0,1 ,.. . ,2k, A(j) = £ (-D V a a, v=0 v J_v B(j)= I (-l)\bj.v v=0 As remarked in section 4.3, the sets R (x) and Q (x) are discre- tized so that | |z (x)| can be approximated by the solution of a convex quadratic programming problem. The discretization contains, in addition to points in R(x) and Q(x), points of the form w±6 in R (x) for each weR(x) and (w± - y. Start with 6 = Jv. If g(x,w-6,j) >_ - y then (w±6,j) are accepted; else set 6 = 6[g(x,w,j) + y]/Cg(x,w,j) - g(x,w-6,j)] and repeat. Empirical results show that the convergence rate is faster if we include the additional test to check if (0,j) is in Q y (x). For given x and w, the vector 3f(x,w)/3x can be found using two 76 linear recurrence relations, the derivation of which is straightforward: 3f 9X (x.w) = - 1 T,R + T 9 I v =1 ,3,. . . where Im r(x,iw)s(x,-iw) 1 _ T 4 R v-l " T 3 (WI -1 + — ) vs2 . 4 -- x v , 1 £ v < k-1 , is the v-th component of x, T, = 2[Re s(x,iw) + F(Im s(x,iw))] T 2 = Im r(x,iw) + F(Re r(x,iw)) T 3 = Re r(x,iw) - F(Im r(x,iw)) T. = 2w[Im s(x,iw) - F(Re s(x,iw))] F = u(x,w) v(x,w) and R v , I v are defined by the recurrence relations R, = l K ' " " w R v-2 v = 3 5 I 1 = w l v = " w ! v-2 v = 3,5,... Similarly, the vector dg(x,w,j)/dx for given (x,w,j) can be ob- tained using the above recurrence relations: f- 4 Re r(x,iw) R v 3g - 4 w Im r(x,iw) R _, - 2 Im s(x,iw) I j = 1 v = 1 , 3 , . . . — (x,w,j) = / 3X V 1 ,2 Re s(x,iw) [wl , + v- I ;+l -] j = 1 v = 2,4, .. . j = 2 v = 1,3,... j = 2 v = 2,4,. . . The program starts with e, - u-i s p-i ■ .001 and decreases the parameters by a factor of 1/2 in each iteration, i.e., = M k /2 e k+l = e k /2 l k+l p k+l = p k /2 77 In case the points generated by the (e,y,p)-algorithm with e = e. , y = y. , p = p. converge to a nonstationary point, the program makes one more at- tempt by setting e. + , = y. + , = p. + , = .1 before an abnormal end is signal- led. Before we can use the M-algorithm to find A*(a), an initial feasible point b = b [a] has to be known. It can be supplied by the user as an input to the program. The program also has the capability of finding b [a] if it is given a point b [A 1 ] representing a method with error con- k i stant of magnitude (a') instead with A' 4 A. The program first applies the M-algorithm for the value A' to obtain a point b* = b*[A']. Then it tests if the vector b' whose v-th component is given by v = 1 ,3, . . . v = 0,2,. . . is feasible for the value A and represents an A(0)-stable method. Mote that from the definition of b', the corresponding method has error constant of magnitude (A) . The above definition for b' works satisfactorily for odd k and for even k when A > A 1 in the sense that in most cases, b' is a feasible point and does represent an A(0)-stable method. For the case when k is even and A' > a, we make use of an empirical observation that b* always tends to zero and hence the polynomial s(b*,z) must have at least one real (negative) root. Instead of fixing A and trying to find b 1 so that the corresponding method has error constant of magnitude (a) , we de- crease the largest negative real root of s(b*,z), say y, by a factor of .1 and define the new polynomial s(b',z) by 78 z - .ly s(b'.z) = s(b*,z) z - y It is straightforward to show that for v = l,2,...,k-l, k b ' = b J + .9y z y J ~ v_1 b* with bg = bg = and that the value A corresponding to b' is smaller than A'. The above definitions for b' are only heuristics for finding an initi- al feasible A(0)-stable method for the M-algorithm. In all cases, if the vector b' is not feasible or does not represent an A(0)-stable method, the program will try to proceed in smaller steps. The positive real roots of a given polynomial are found by first generating the Sturm sequence for the polynomial and its derivative in the interval [0,°°) and then using the ZEROIN routine (Shampine and Allen, 1973, pp. 244-246) which locates a root of a continuous function inside a given interval by a combination of bisection and secant rules. The convex qua- dratic programming problem is solved using the first symmetric variant of the simplex method (Van De Panne, 1975, pp. 270-280). Numerical results for A* (A) for a finite number of points on < A < oo for the cases when k = 4,5,6,7 are plotted in figures 5.5, 5.6, 5.7, 5.8, respectively. The isolated symbols + represent the (A, a ) values with 1/k A = |c k+ i I and a being the angle of absolute stability for the following methods: BDF (Backward Differentiation Formula), CHEB2, CHEB3, CHEB4 (Gup- ta, 1976), and It, L k (Gupta and Wallace, 1975). The BDF of order 6 and the other methods in Gupta (1976) lie outside the region described in the plots and hence are not included. The table following each figure contains the values of the parameters bQ»b-i ^k-1 ^ or t ^ ie metnoc * corresponding to 79 90° f 80 c 70 c 60 c 50 ( 40< 30° - 20 c CHEB4 10 c 0.35 0.45 0.55 0.65 Figure 5.5. A*(a) for k = 4 0.75 0.85 80 u u l u 2 u 3 .0 .0022 .4165 .6103 .0 .0041 .4526 .6241 .0 .0066 .4913 .6389 .0 .0094 .5284 .6530 .0 .0136 .5787 .6716 .0 .0198 .6411 .6953 .0 .0282 .7159 .7226 .0 .0383 .7963 .7513 .0 .0462 .8544 .7714 .0 .0597 .9458 .8025 .0 .0950 1.1544 .8700 .0 .1164 1.2675 .9048 .0 .1428 1.3980 .9434 .0 .2021 1.6630 1.0177 .0 .2864 1.9979 1.1049 .0 .4367 2.5206 1.2290 .0 .6673 3.2156 1.3773 .0 1.2635 4.7096 1.6522 .0 3.5655 9.0 2.2637 .0 13.2348 21.0 3.4392 Table 5.1 k = 4 81 90° 4- pur 80° LHt 70° • + L 5 60° of + 50° / + / BDF 40° - / + / CHEB2 30° 20° 10° n° 8 1 1 1 1 \ 1 0.4 0.5 0.6 0.7 0.8 Figure 5.6. A*(a) for k = 5 0.9 1.0 82 u .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 D l 1.366 1.366 1.366 1.366 1.366 1.366 1.675 3.081 6.553 10.660 15.241 20.226 31.206 37.131 49.725 70.236 82.374 113.365 146.699 182.088 d 2 .816 .817 .819 .834 .863 .981 1.998 3.359 6.171 9.036 11.929 14.837 20.687 23.622 29.508 38.363 43.292 55.141 67.011 78.895 D 3 3.238 3.238 3.238 3.238 3.238 3.238 3.568 4.709 6.733 8.518 10.141 11.648 14.417 15.707 18.144 21.527 23.296 27.296 31.025 34.544 .373 .373 .374 .381 .394 .448 .837 1.068 1.381 1.606 1.786 1.938 2.189 2.296 2.487 2.729 2.847 3.099 3.316 3.508 Table 5.2 k = 5 83 80 o . . 70° - 60 c 50° - 40° 30° 0.55 0.65 0.75 Figure 5.7. A*(a) for k = 6 CHEB2 0.85 84 u .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 u l .233 .298 .385 .500 .653 .749 .860 1.139 1.314 1.756 2.034 2.741 3.188 4.713 6.082 7.871 12.836 17.215 26.848 37.431 u 2 3.470 3.928 4.491 5.188 6.054 6.565 7.137 8.496 9.302 11.227 12.375 15.130 16.779 22.032 26.405 31.780 45.421 56.451 78.745 101.254 u 3 3.579 3.882 4.244 4.678 5.200 5.500 5.829 6.586 7.021 8.026 8.605 9.943 10.716 13.061 14.907 17.076 22.209 26.077 33.367 40.213 u 4 4.218 4.459 4.738 5.061 5.436 5.645 5.869 6.371 6.650 7.273 7.620 8.392 8.821 10.065 10.992 12.033 14.334 15.951 18.794 21.280 u 5 1.295 1.329 1.367 1.410 1.459 1.486 1.514 1.576 1.609 1.682 1.721 1.806 1.852 1.978 2.068 2.164 2.363 2.494 2.709 2.884 Table 5.3 k = 6 85 80° 70° 60° 50 c 0.7 0.8 0.9 Figure 5.8. A*(a) for k = 7 1.0 86 u .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 .0 u l 9.750 14.853 46.134 133.735 151.480 178.912 202.449 226.560 251.198 276.330 301.921 327.951 354.396 381.237 436.028 u 2 26.031 35.990 87.855 206.759 228.727 261.810 289.483 317.238 345.065 372.957 400.908 428.912 456.964 485.061 541.376 29.627 37.941 75.803 148.501 160.819 178.928 193.719 208.269 222.599 236.729 250.674 264.450 278.068 291.539 318.076 D 4 15.260 18.553 31.785 53.244 56.584 61.381 65.211 68.907 72.485 75.959 79.337 82.628 85.840 88.980 95.061 D 5 10.578 11.942 16.798 23.442 24.388 25.716 26.751 27.732 28.665 29.556 30.410 31.230 32.020 32.783 34.236 D 6 1.898 2.048 2.506 3.020 3.086 3.176 3.245 3.308 3.368 3.424 3.477 3.526 3.574 3.619 3.703 Table 5.4 k = 7 87 the symbol a in ascending order of A. The numerical values for (a,A*(a)) 1/k and the corresponding (rG) and r are listed in Appendix B. 5.4 CONCLUSION From the empirical results, it is observed that the parameter b Q corresponding to the method found by the M-algorithm is always zero (or very close to zero) which implies that a(c) has a root at e = -1. Moreo- ver, except for the root z, = 1 of p(0 and the root ? = -1 of a(s), the other roots of p(c) and a(s) are close together and they seem to be ap- proaching the point c = 1 as A -> °°. It is interesting to note that for all the methods found, the ratio (rG) to |c k+ -.| is not very large (< 2), and that r < 1.3. But each of the methods has a pair of roots of p(c) of magnitude close to 1. Thus, the global error accumulated when any of the methods is used is not that much affected by the positions of the roots of p(c) provided a fixed stepsize is used throughout. We could compare the BDF with the methods found by the ri- al gorithm as follows: Since the k-th order k-step BDF has error constant of magnitude l/(k+l), we consider the ratio H 2 " a BDF[k] \ - A*((l/(k+l)) 1/k ) where ^BDFfkl 1S tlie an 9^ e °f absolute stability associated with the BDF. For 4 £ k i 6 (the 7-th order 7-step BDF is not A(O)-stable), the above ra- tio is greater than 3. In other words, the method found by the M-algorithm having error constant of magnitude l/(k+l) is at least 3 times "closer to being A-stable" than the BDF. On the other hand, for 4 <_ k <_ 6, the ratio 0/(k+l)) 1/k 88 where a. is the value of a such that A*(a. ) = a RnF r k -|. Thus for the method corresponding to (a, ,A*(a.)) which has the same angle of absolute stability as the BDF, the error constant has magnitude less than (10/13) times that of the BDF. This implies that if, in an ODE solver, the stepsize is chosen to be proportional to the k-th root of the reciprocal of the magnitude of the error constant, then the method corresponding to (a.,A*(a.)) is 1.3 times more efficient than the BDF. However, numerical tests (based on DIFSUB (Gear, 1971) with appropriate modifications) show that for problems which can be successfully handled by the BDF, the method corresponding to (a, ,A*(a. )) is only about .7 times as efficient as the BDF. The reason may be that our theory did not take into account the local instabilities that could arise from step changing and order changing. This is especially likely since the p(c) and a(c) polynomials associated with the methods have roots which are much nearer the unit circle than those of the BDF. A final point is that we can see from the plots that some of the methods obtained by Gupta and Wallace (1975) and Gupta (1976) are close to the (a,A*(a)) values. 89 BIBLIOGRAPHY BabuSka, J.; Prager, M.; and Vitcisek, E., 1966, Numerical Processes in Dif- ferential Equations , SNTL — Publishers of Technical Literature, Prague. Bjurel , G; Dahlquist, G. ; Lindberg, B. ; Linde, S.; and Oden, L., 1970, "Survey of Stiff Ordinary Differential Equations," Report No. NA70.11, Dept of Information Processing, Royal Inst, of Technolo- gy, Stockholm, Sweden. Coddington, E. A., and Levinson, N., 1955, Theory of Ordinary Differential Equations . McGraw-Hill Book Co., New York. Cryer, C. W. , 1973, "A New Class of Highly-Stable Methods: A n -stable Methods," ITT, vol. 13, pp. 153-159. Dahlquist, G. G. , 1963, "A Special Stability Problem for Linear Multistep Methods," BTT, vol. 3, pp. 27-43. Dem'yanov, V. F., and Malozemov, V. N., 1974, Introduction to Minimax , John Wiley and Sons, Inc., New York. Dill, C, and Gear, C. W. , 1971, "A Graphical Search for Stiffly Stable Methods for Ordinary Differential Equations," Journal of the ACM , vol. 18, pp. 75-79. ' Duff in, R. J., 1969, "Algorithms for Classical Stability Problems," SI AM Review , vol. 11, pp. 196-213. Gantmacher, F. R. , 1959, Theory of Matrices , vol. II, Chelsea Pub. Co., New York. Gear, C. W., 1971, "Algorithm 407 - DIFSUB for solution of Ordinary Dif- ferential Equations," Communications of the ACM , vol. 14, pp. 176-179. Gear, C. W. , 1971, Numerical Initial Value Problems in Ordinary Differenti- al Equations , Prentice Hall, Inc., Englewood Cliffs, New Jersey. Genin, Y., 1973, "A New Approach to the Synthesis of Stiffly Stable Linear Multistep Formulas," IEEE Transactions on Circuit Theory , vol. CT-20, pp. 352-360. Gupta, G. K. , 1976, "Some New High-Order Multistep Formulae for Solving Stiff Equations," Mathematics of Computation , vol. 30, pp. 417-432. Gupta, G. K., and Wallace, C. S., 1975, "Some New Multistep Methods for Solving Ordinary Differential Equations," Mathematics of Computa- tion , vol. 29, pp. 489-500. 90 Henrici, P., 1962, Discrete Variable Methods in Ordinary Differential Equa- tions , John Wiley and Sons, Inc., New York. Jain, M. K. , and Srivastava, V. K. , 1970, "High Order Stiffly Stable Methods for Ordinary Differential Equations," Dept. of Computer Science Report No. 394, University of Illinois, Urbana, Illinois. Jeltsch, R., 1976, "Stiff Stability and Its Relation to A Q - and A(0)- Stability," SIAM Journal on Numerical Analysis , vol. 13, pp. 8-17. Krall, A. 11., 1967, Stability Techniques for Continuous Linear Systems , Gordon and Breach, Science Publishers Inc., New York. Lambert, J. D., 1973, Computational Methods in Ordinary Differential Equa- tions . John Wiley and Sons Ltd, London. Lehnigk, S. H., 1966, Stability Theorems for Linear Motions with an intro- duction to Liapunov's Direct Method , Prentice Hall, Inc., Englewood Cliffs, New Jersey. Lindberg, B., 1974, "On a Dangerous Property of Methods for Stiff Differen- tial Equations," BII, vol. 14, pp. 430-436. Liniger, W., 1975, "Connections Between Accuracy and Stability Properties of Linear Multistep Formulas," Communications of the ACM , vol. 18, pp. 53-56. Luenberger, D. G., 1969, Optimization by Vector Space Methods , John Wiley and Sons, Inc, New York. Levinson, N., and Redheffer, R. M., 1970, Complex Variables , Holden-Day, Inc. , San Francisco. Marden, M. , 1966, Geometry of Polynomials , American Mathematical Society, Providence, Rhode Island. Mitrinovic, D. S., 1970, Analytic Inequalities , Springer-Verlag, New York. Odeh, F., and Liniger, W. , 1971, "A Note on Unconditional Fixed-h Stability of Linear Multistep Formulae," Computing , vol. 7, pp. 240-253. Osborne, M. R. , 1966, "On Nordsieck's Method for the Numerical Solution of Ordinary Differential Equations," jSIT, vol. 6, pp. 51-57. Polak, E., 1971, Computational Methods in Optimization: A Unified Approach , Academic Press, New York. Shampine, L. F., and Allen, R. C, Jr., 1973, Numerical Computing: an in- troduction . W. B. Saunders Co., Philadelphia. Skeel, R. D., 1977, "Equivalent Forms of Multistep Formulas," In prepara- tion, Dept. of Computer Science Report, University of Illinois, 91 Urbana, Illinois. Stetter, H. J., 1973, Analysis of Discretization Methods for Ordinary Dif- ferential Equations , Springer-Verlag, New York. Van De Panne, C, 1975, Methods for Linear and Quadratic Programming , North-Holland Pub. Co., Amsterdam. Wallace, C. S., and Gupta, G. K. , 1973, "General Linear Multistep Methods To Solve Ordinary Differential Equations," Australian Computer Journal , vol. 5, pp. 62-69. Widlund, 0. B., 1967, "A Note on Unconditionally Stable Linear Multistep Methods," BIT, vol. 7, pp. 65-70. Zoutendijk, G., 1966, "Nonlinear Programming: A Numerical Survey," SI AM Journal on Control , vol. 4, pp. 194-210. 92 APPENDIX A POLYNOMIALS HAVING ROOTS IN THE OPEN LEFT HALF COMPLEX PLANE 93 The following is a list of equivalent conditions each of which is necessary and sufficient for a polynomial , x n n-1 p(z) = c n z + c n _-jz + ... + C-|Z + CQ c n > to be Hurwitzian, i.e., to have roots all in the open left half complex plane: (1) (Hurwitz criterion) The leading principal minors of the nxn matrix H whose (i,j)-th element is given by H-jj = c n _.j_2j i»j = l,2,...,n are positive (Krall, 1967, p. 48). A stronger statement (Lienard and Chi- part criterion) is in Gantmacher (1959, p. 221). (2) (Routh criterion) The elements in the first column of the tableau R are positive, where the elements in R are generated by the fol- lowing scheme (Krall, 1967, p. 52): Start with the first two rows whose elements are c n c n-2 c n-4 c n-6 ••• c n-l c n-3 c n-5 c n-7 ••• Every row after the first and second in the tableau is determined by the two preceding rows according to the following rule: If the elements of the two preceding rows are §0 9] §2 93 hg h-i hp ho ... then the elements q ,q-,,... of the current row are given by % q j = g j+l " T-Vl j = °' 1 ' — n 94 The tableau consists of the first n+1 rows generated by the scheme. (3) (L. Cremer-Leonhard separation criterion) The zeros {w v > and {w^} of the polynomials Re p(i/^w) and Im pCi^w)//^, respectively, sa- tisfy > w^ > w-| > W2 > W2 > ... (Lehnigk, 1966, p. 154). (4) The polynomial 2 n-1 , . n-2 n-3 c n-l z + ( c n-l c n-2- c n c n-3) z + c n-fn-f + + (c n „ 1 c n _ 4 -c n c n _ 5 )z n_4 +'... + n-2m+l / \ n-2m + c n-l c n-2m+l z + ^ c n-l c n-2m" c n c n-2m-V z + ••• is Hurwitzian (Krall, 1967, p. 47). We next give some necessary conditions for p(z) to be Hurwitzian. Suppose that p(z) is Hurwitzian. Then it is obvious that Cj > j = 0,1 ,...,n Using (4), it follows that c -, c .-cc . -, > i = 2,4,... n-1 n-j n n-j-1 J •"'»•.•■ Moreover, from (3), using the Dougall's inequality for symmetric functions (Mitrinovic, 1970, p. 97), it can be shown that the following inequalities are satisfied: 2 c 2 4 c 4 n-2 c n _ 2 n c n > . . . > > — n c Q n-2 c 2 4 c„. 4 " 2 c n . 2 2 c 3 4 c 5 n " 4 c n-3 n " 2 c n-l n-2 c, - n-4 c - ~ 4 c c _ 2 c Q 3 n-b n-J 95 when n is even, and 4 c 4 n-3 c n . 3 n-1 c„_ 1 > > • . • > > n-1 c. n-3c " 4 c c 2 c 2 n-5 n-3 2 c_ 4 c c n-3 c n-1 c 3 5 n-2 n > > • • • > n-1 c, n-3c, 4c. 2 c 3 n-4 n-2 when n is odd. APPENDIX B LIST OF A, A*(A), (rG) 1/k , f VALUES FOR METHODS FOUND USING THE M-ALGORITHM 96 97 A A*(A) (rG)"* r .3815 12.17 .5069 1.1388 .3848 16.45 .5086 1.1381 .3883 20.55 .5103 1.1374 .3916 24.06 .5118 1.1366 .3959 28.24 .5137 1.1353 .4010 32.87 .5158 1.1336 .4069 37.53 .5182 1.1314 .4130 41.78 .5205 1.1288 .4172 44.43 .5221 1.1268 .4236 48.08 .5243 1.1236 .4372 54.59 .5287 1.1162 .4441 57.34 .5307 1.1123 .4517 60.03 .5328 1.1078 .4660 64.32 .5364 1.0993 .4823 68.27 .5450 1.0896 .5050 72.47 .5692 1.0769 .5310 76.06 .5919 1.0636 .5767 80.32 .6403 1.0445 .6687 84.85 .7323 1.0204 .8190 87.77 .8888 1.0098 Table B.l k = 4 98 A A* (a) (rG) 1 ^ f .4045 .01 .6233 1.2074 .4046 .07 .6233 1.2073 .4048 .20 .6232 1.2072 .4063 1.22 .6228 1.2066 .4091 3.21 .6218 1.2052 .4197 10.47 .6184 1.2001 .4821 43.49 .6238 1.1649 .5296 57.68 ' .6602 1.1321 .5923 68.83 .7300 1.0937 .6361 73.72 .7776 1.0716 .6704 76.54 .8206 1.0606 .6988 78.40 .8517 1.0563 .7448 80.75 .9145 1.0484 .7641 81.54 .9385 1.0451 .7977 82.73 .9847 1.0393 .8394 83.91 1.0417 1.0326 .8594 84.38 1.0708 1.0296 .9010 85.23 1.1315 1.0256 .9361 85.81 1.1848 1.0231 .9666 86.25 1.2300 1.0210 Table B.2 k = 5 99 A A* (a) (rG) ,/K r .5677 40.52 .7072 1.1858 .5763 43.54 .7139 1.1807 .5858 46.59 .7234 1.1750 .5965 49.67 .7280 1.1687 .6085 52.72 .7412 1.1617 .6149 54.23 .7467 1.1579 .6217 55.72 .7548 1.1540 .6363 58.64 .7662 1.1457 .6442 60.06 .7694 1.1414 .6610 62.80 .7915 1.1324 .6700 64.12 .7986 1.1278 .6893 66.65 .8169 1.1182 .6995 67.85 .8291 1.1134 .7230 70.29 .8492 1.1031 .7473 72.42 .8740 1.0930 .7601 74.00 .8985 1.0852 .8106 76.66 .9420 1.0711 .8379 78.06 .9716 1.0632 .8821 79.91 1.0166 1.0522 .9174 81.12 1.0570 1.0448 Table B.3 k = 6 100 A A* (A) (rG) ,/N r .7131 54.28 .9039 1.1676 .7430 58.98 .9322 1.1524 .8343 68.85 1.0306 1.1122 .9351 75.26 1.1381 1.0785 .9479 75.87 1.1545 1.0749 .9653 76.64 1.1722 1.0703 .9786 77.19 1.1875 1.0669 .9908 77.67 1.2021 1.0640 1.0022 78.09 1.2134 1.0613 1.0129 78.47 1.2270 1.0589 1.0229 78.81 1.2374 1.0568 1.0324 79.12 1.2483 1.0557 1.0414 79.41 1.2595 1 . 0548 1.0499 79.66 1.2674 1.0540 1.0658 80.12 1.2865 1.0524 Table B.4 k = 7 101 VITA Antony King- Yin Kong was born on February 10, 1953 in Hong Kong. After finishing secondary school in his own country, he attended Iowa State University of Science and Technology from 1970 to 1973. He was elected to the Upper 2% Scholars and was awarded the Gertrude Herr Adamson Mathematics Awards. He received the degree of Bachelor of Science in Computer Science, Mathematics, and Statistics, with distinction, in 1973. In the same year, he was elected a Phi Kappa Phi Fellow to study in the Department of Compu- ter Science in the University of Illinois. He served as a research assis- tant from May of 1974 until September of 1977 and coauthored with Professor Robert D. Skeel in the paper titled "Blended Linear Multistep Methods." »rmAEC-427 (6/68) AECM 3201 U.S. ATOMIC ENERGY COMMISSION UNIVERSITY-TYPE CONTRACTOR'S RECOMMENDATION FOR DISPOSITION OF SCIENTIF C AND TECHNICAL DOCUMENT / See Instructions on Reverse Side ) AEC REPORT NO. C00-2383-0046 ' A SEARCH FOR BETTER LINEAR MULTISTEP METHODS FOR STIFF PROBLEMS TYPE OF DOCUMENT (Check one): |X| a. Scientific and technical report ^] 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) RECOMMENDED ANNOUNCEMENT AND DISTRIBUTION (Check one): IX1 a. AEC's normal announcement and distribution procedures may be followed. 71 b. Make available only within AEC and to AEC contractors and other U.S. Government agencies and their contractors. J c. Make no announcement or distrubution. REASON FOR RECOMMENDED RESTRICTIONS: SUBMITTED BY: NAME AND POSITION (Please print or type) C. W. Gear Professor and Principal Investigator Organization Department of Computer Science University of Illinois Urbana,/i,llinois 61801 FOR AEC USE ONLY AEC CONTRACT ADMINISTRATOR'S COMMENTSTlF ANY, ON ABOVE ANNOUNCEMENT AND DISTRIBUTION RECOMMENDATION: PATENT CLEARANCE: LJ a. AEC patent clearance has been granted by responsible AEC patent group. Lj b. Report has been sent to responsible AEC patent group for clearance. LJ c. Patent clearance not required. IBLIOGRAPHIC DATA HEET 1. Report No. UIUCDCS-R-77-899 itle anJ Subt itlc A SEARCH FOR BETTER LINEAR MULTISTEP METHODS FOR STIFF PROBLEMS 3. Recipient's Accession No. 5. Report Date December 1977 6. Author(s) Antony King-Yin Kong 8. Performing Organization Rept. No - UIUCDCS-R-77-899 Performing Organization Name and Address Department of Computer Science University of Illinois Urbana, Illinois 61801 10. Project/Task/Work Unit No. COQ-2383-0046 11. Contract/Grant No. US ERDA/EY-76-S-02-2383 !. Sponsoring Organization Name and Address US Energy Research and Development Administration 9800 South Cass Avenue Argonne, Illinois 60439 13. Type of Report & Period Covered thesis 14. i. Supplementary Notes i. Abstracts For arbitrary k > 1 and ae(0,Tr/2), A(a)-stable k-th order k-step formulas exist, so that in an- ODE solver, a can be an extra parameter used to identify among a family of methods of order k the A(a)-stable method that should be used for the particular problem. Two measures for assessing the accuracy of k-th order k-step formulas are proposed. The problem of finding the upper bound on the angle of absolute stability for the k-th order k-step formulas having the same accuracy (with respect to one of the measures) is considered. Analytical results are obtained for k = 1, 2, 3 whereas a numerical search is used for the cases when k = 4, 5, 6, 7. Key Words and Document Analysis. 17o. Descriptors linear multistep method stiff ordinary differential equations b. Identifiers/Open-Ended Terms |e. COSATI Field/Group • Availability Statement unlimited ''"M NTIS-35 I 10-70) 19. Security Class (This Report) UNCLASSIFIED 20. Security Class (This Page UNCLASSIFIED 21. No. of Pages 105 22. Price USCOMM-DC 40329-P7 1 CCD o 1Q7ft AUG 15