CENTRAL CIRCULATION BOOKSTACKS The person charging this material is re- sponsible for its renewal or its return to on n- ?7 ^ Whkh U Was borrowed on or before the Latest Dote stamped £. 3 . Carver, Chalk River Nuclear LaDoratories In search of a robust integration algorithm for general library use: some tests, results and recommendations .... 36-1 r. £. Zadunaisky, CNI£-Gbservatoric Nacional de Fisica Cosmica and G. Lafferriere, University of Buenos Aires On an iterative improvement of the approximate solution of seme ordinary differential equations 37-1 i. ¥, Chang and R. horris, University of Nebraska — Lincoln and G. Corliss, harcuette University f. portable automatic Taylor series (ATS) compiler for solving ordinary differential equations 3Q-1 r. Kaps, university of InnsbrucK Rosenbrock-type methods for the numerical solution of stiff systems 39-1 - . S. WatKins, Utah State University Determining realistic initial values for stiff systems cf ordinary differential eouations occurring in ionospheric physics 40- 1 . ■.. Paine, Australian National University A comparison of methods for computing Sturm-Liouville eigenvalues 41-1 iii WORKSHOP Technique? tor cifficult problems (e.g. oscillatory systems) no paper D. G. bettis, moderator, University of Texas at Austin F. Deuflhard, University of Heidelberg C. W. Gear, University of Illinois at Urbana-Champaism L. R. Petzold, Sandia Laboratories, Livermore lhe numerical solution of second order systems of ODEs 43-1 { / t . H. Enright, moderatcr, University of Toronto D. G. tettis, University of Texas at Austin P, Deuflhard, university of Heidelberg F. T. krogh, Jet Propulsion Laboratory Method cf Lines 44-1 G. o. Eyrne, moderator, University of Pittsburgh M. £. Carver, ChalK River Nuclear Laboratories A. C. hindmarsh, Lawrence Livermore Laboratory G. K. Leaf, Argonne National Laboratory M. Minkoff, Argonne National Laboratory in. L. Schiesser, Lehigh University R. F. Sincovec, Boeing Computer Services Implementation of implicit Runge-Kutta codes no paper F. h. Chipman, moderator, Acadia University T. A. tickart, Syracuse University J. C. Butcher, University of Auckland M. A. Epton , Boeing Computer Services S. P. Nrirset, University of Trondheim L. F. Shampine, Sandia Laboratories, Albuquerque Panel en testing 45-1 T. E. Hull, moderator, University of Toronto h. t. Carver, Chalk River Nuclear Laboratories k. ri. Enright, University of Toronto I. Gladwell, university of Manchester r. I. Krogh, Jet Propulsion Laboratory L. t . Shampine, Sandia Laboratories, Albuquercue A user interlace standard for ODE solvers 46-1 A. C. hindmarsh, Lawrence Livermore Laboratory An incomplete list of better E0R1RAN codes for IVPs in ODEs 47-1 R. u. Skeel, University of Illinois at Urbana-Champaign IV CONFERENCE OVERVIEW R. D. Skeel Program Committee Chairperson This meeting is another in a series of single-topic SIGNUM Meetings. It differs from most of its predecessors in two ways. First, it is devoted to a well established and well defined area of numerical mathematics; and although the scope of the meeting is quite narrow, it has a large vertical extent ranging from the very theoretical down to the very practical. Second, by publicizing the meeting well in advance and by obtaining the generous support of the National Science Foundation, it was possible to make the meeting truly international. The numerical solution of the initial value problem for ordinary differential equations is often cited as one of the notable successes of numerical mathematics. Nevertheless, there are evidently many who believe that much can still be done, and the significant advances being reported at this meeting confirm that belief. The planning of the program has been influenced to some extent by the following categorization of research efforts in numerical ODEs: (i) theory — the analysis of numerical methods and techniques, (ii) algorithms — the synthesis of numerical methods and techniques, (iii) software — implementation of methods and techniques as general purpose library routines, (iv) applications — design of algorithms and tailoring of software for specific areas of application. These various aspects of numerical ODEs might be visualized as a triangle with algorithms in the center, theory at the top corner, and software and applications at the two supporting bottom corners. A rough breakdown of symposium presentations is fifteen in theory, fourteen in algorithms, five in software, and five in applications. This preponderance of more theoretical papers is balanced by the more practical orientation of the workshop. The organization of the technical program is in a large part due to the efforts of T. E. Hull, J. D. Lambert, L. F. Shampine, and H. J. Stetter in reviewing the summaries for about thirty-four papers and I would like to thank them very much. Also appreciated is the sponsorship by Tom Hull of the OlympiODE. I would like to acknowledge the efforts of Bill Gear and Dan Watanabe in the overall organization of the conference as well as their advice on the planning of the program. It should also be noted that Bill organized the workshop. The assistance of Dennis Gannon and Mitch Roth with local arrangements is gratefully acknowledged. The University of Illinois Department of Computer Science has supported this conference in numerous ways. Barb Armstrong has done a great deal of secretarial work, and business manager Clyde Helm has been very helpful. Finally, we are very pleased to have so many of the top researchers from throughout the world participate in this meeting. The efforts and expenditures of all participants is much appreciated. SOME PROPERTIES OF LINEAR MULTISTEP AND ONE-LEG METHODS FOR ORDINARY DIFFERENTIAL EQUATIONS Germund Dahlquist Royal Institute of Technology S-IOOI4I* Stockholm, Sweden 1. I'hat is a one-leg method and what relation does it have to a linear multisteD method? This paper deals with the numerical integration of the initial value problem for the differential system, (1.1) dy/dt = f(y), y£f s . To simplify the writing we restrict the discussion to autonomous systems. The generalization of re- sults for non-linear systems to the non-auto- nomous case is straight-forward, see e.g. Dahlquist 0978). In the constant stepsize case, t n = t +nh, a one-leg method is written, k k (1-2) I «,r^ «hf( I 6,y n+ ,) , j=o J n+J j=0 < where Jb.= 1. It is "The one-leg twin" of the more familiar linear multistep method, {K3) J o °J P n + =h J o B J f( ^)- Some one-leg methods have been in use for a long time. One of the simplest examples is the implicit midpoint method, Vl-^^n+l+V)' which is the one-leg twin of the trapezoidal method, The backward differentiation methods are their own one-leg twins. In terms of the displacement operator E, Ey a y ... and the generating polynomials, n k k pU) = I a. ? j , a(c) = I B.? j , j=o J=0 (1.3*) oU)9 = ha(E)f(y ) n n the methods can be written, (1.2') p(E)y = hf (o(E)y ) (one-leg) n n' i (linear multistep) e assume that p,o are relatively prime. The normalization convention and the consist- ency conditions read, (1.1») 0(1) = 1, (normalization) (1.5) p(D=0, p'(D = 0(1) (consistency) For constant stepsize the following relations between the twins are well known, Dahlquist (1975a). a) LINEAR IDENTITY. They yield the same difference equation for a linear autonomous system. b) NON-LINEAR EQUIVALENCE (first part). If {y } satisfies the one-leg equation (1.2') and if y n = o(E)y n , then {y } satisfies the linear multistep equation (l.3'). n In order to prove the second statement , one multiplies (1.2') by the operator a(E) and makes use of aommutativity , o(E)p(E) = p(E)a(E). Dahlquist ( 1975a) introduced one-leg methods, since they gave neater formulations in the analy- sis of numerical stability for non-linear systems. There are, however, other reasons for the consider- ation of these methods. For example, they give an efficient implementation of the linear multistep twins. Note that one needs k vectors only, v n> v n+1 '• • • '^n+k-1 » in order to compute y n+k from (1.2). Then y n = a(E)y n is computed but not stored. The most straightforward implementation of (1.3) needs 2k values from the past (except for "sparse" methods like Adams or BDF), if one does not re- compute k old values of f(y). However, Skeel (1978) has recently found other implementations of a general linear multistep method with the same advantage as the one-leg implementation. A one-leg method may also be considered as a numerical method in its own right. Nevanlinna and Liniger (1979) point out that it is advantegeous compared to a linear multistep method in variable step computations. On the other hand, the latter has, in a sense, usually higher accuracy than its one-leg twin, see Section 2. In the variable stepsize case, let h m =t m -t m -i. We temporarily write a one-leg method in the form, where the aj n ,Bi n , depend on the stepsizes • ' n n+k' We fi- 11 ^' by dimensional ana- ' h n+2- "n+1 lysis, that it is sufficient to consider methods of the form, (1.6) P n (E) ^n = h f (0 (E)y ) , in ^ n r\ J n a h > n where the coefficients of p n -!> J = 2,3,. depend on the homogeneous of the first degree in 3 h n+j /h n+j ,k, only, and h n is 1-1 h „ h ., . The consistency conditions read + ' ' - + ' • n+K p n (l) =0, P n (E)t n =h n . The method (1.6) is normalized if, .8) P ;(D = a n (D = 1 • The condition p ' ( 1 ) = a ( 1 ) is introduced for for- mal convenience n only. It means no restriction, since it will just give a natural scaling of & n> Note that in general, p n (E) and ojE) don't rte. Therefore, the short proof of non-linear equivalence given above, does not work for arbi- trary changes of stepsize. In fact, it _ follows from an argument of Nevanlinna and Liniger (1 979, p.6U) that there is no non-linear equivalence for arbitrary stepsize sequences. However, for geometric stepsize sequences, i.e. if (1.9) h n = h • x n , the polynomials p n ,o n do not depend on n (although they depend on x ) . In this case p n (E) and o n (E) commute. The following results are proved very much like the special case in Dahlquist (1975a). PROPOSITION 1. Linear identity for arbitrary step- size sequences. For a linear autonomous system the one-leg method (1.6) gives the same difference equation as the linear multistep method, J. 10) P n (E)y n = h n a n (E)f(y n ) . Non-linear equivalence for geometric stepsize sequences defined by (1.9). Assume that {y n } satisfies the one-leg difference equation, (1.11) p(E)y n = v n , and let y^ = o(E)y n . Then (1.12) p(E)y n = a(E)v n . If v =h f(y ) we can write (1.12) in either of n n n the forms , (1.13) p(E)? n = a(E)(h n f(y n )) , (1.13') p(E)y n = h n a(tE)f(y n ) . Conversely, assume that for n*0, y n satisfies (1.12). Let P m , ^ be polynomials of degree less than k, such that for some m, 0■'" ■ 2. The local error Let y(t) be an exact solution of the differential system ( 1 . 1 ) , and define the local error of the one-leg method (1.2') by, (2.1) e = h^p (E)y(t )-f(a (E)y(t )). n n n n n " The displacement operator E operates on a seque space, Ey m Here (Ey) = y m+1 , (not on a function space" J Q a jn y m + j P n (E) ^ and similarly for a n (E). Introduce the operators L ip(t ) = p (E)tp(t )-h a (E)«p'(t L , S , n n [2.2 n n n n' S cp(t ) = a (E)jp(t )- 3 only. The relations (2.8) were first obtained by Liniger (1978), for constant stepsize, after that Warming and Dahlquist had realized that in general the order will not exceed 2 for one-leg methods. In view of the non-linear equivalence (for constant stepsize) it may seem paradoxical that these methods do not have the same "order of However, Prop.1 can only accuracy prove that a(E)y follow that y be used to 0(h v ). It does not y(o(E)t ) n „/ n rv\ ■ • -y(t ) = O(hP). Similar observations n n have been made for other numerical methods, see Butcher (1973). 3. Stability and contractivity Consider the application of the method (p,o) with constant step h to the test equation, y' = Xy, and set Y p = (y n »y n+l .....y 5 ^. 1 ) T . The difference equation can then be written in vector-matrix (3-D n+1 A(Xh)Y where A(Xh) is a k*k companion matrix. We recall that the stability set is S = jxheC: {A n (Xh)}" is bounded!. The contractivity set K, Nevanlinna and Liniger (1978), with respect to a given norm 11*11 in C k is {XhCC: HA(Xh)ll < 1} . *) Obviously KcS. We use terms like A-contractivity etc. in obvious analogy with A-stability etc. We define the G~norm of a vector, lv||2 = v gv, v e C. u ■ A corollary is that A-stability implies A-con- tractivity for some G-norm (a property previously called G-stability ) . Unfortunately, it also follows that for a method that is not A-stable, one cannot find a G such that Kq contains both and °°. For other norms K is usually not circularly bounded. When a one-leg method is applied to the equation y' = X(t)y, the difference is still of the form (3.1), and if X(t)h€K for all t, then llA(Xt n )ll< 1 , Vn. For a linear multistep method, the difference equation becomes more complicated. The matrix A depends on the values of X at several points. This example indicates both why one-leg methods may be easier to analyse than linear mul- tistep methods in more general situations, and why contractivity is easier to generalize than plain stability. For G-~norms , there is a straight-forward gener- alization to nonlinear systems with a certain monotonicity property related to Kq, Dahlquist ( 1975b, 1978). The theory is given for constant stepsize only, but the extension to arbitrary step- size sequences offers no difficulty, provided that there is a G that can be used for all occuring (p ,a ). We shall see in Section k that there exist one-leg 2-step methods with this property. Methods, which do not posess this property, may be open for an analysis, based on estimation of the rate of change of G. For maximum norms , Nevanlinna and Liniger (1978,1979) have obtained very interesting results . In non-linear problems they assume a kind of diag- onal dominance for the Jacobian. Their analysis is normally restricted to a smaller class of methods than the G-norm theory, but in some respects the maximum norm theory is stronger, e.g. in connexion with Ag-contractivity. See also the end of Section h. 4. A-contractivity of two-step methods for arbitrary stepsize sequences Let P n U) a n U) = a 2n ?2 + a 1n ? + a 0n' 2-cC "In"' "On Assume that the method is A-stable for each n,i.e. 5 >1 Rea (r.)/p ( C ) > 0, Vn . n n For brevity we drop n in the subscripts of the coefficients. It must be remembered, however, that they depend on h ,p/ n + -,• LEMMA 1. The consistent, normalized A-stable two- step methods can be expressed in terms of three non-negative parameters, a,b,c, 2a 2 = c + 1 (1*.1) 2d! = -2c 2a = c - 1 1*62 = 1 + b + (a+c) , l*6i = 2(l-b) UBq = 1 + b - (a+c). 1-3 . V = 1 26!, a+ c = 2(8? - 6 ). Formation = .-+ -1). Then oU)/pU) can be ex- ■ ■ i ■ z+c L , 1 , C > . = " + 0/(5-0 and match coefficients, the normalization (1.8). We iow use the notations and results of 2 to express second order accuracy. It is enient to put t£ = t n+ i and express stepsize rms of the parameter, a, defined thus, :. = o(t 2 + t )/2 . a > means increase of stepsize. We substitute this and (U.1) into (2.7). It is :'icient to consider q= 1 and q=2. We obtain r some computation: t a(1-o)c * D. The method defined by able and second-order accurate, iff — + — =1-c 2 , a>0, b>0, c>0. a 1-a = the condition reads a=0. If a = 1 and c>0 the condition reads b=0. If c = then (l-b)(T 2 +T ) = 2a, x 2 - t = 2. ■ In order to facilitate the comparison with a more familiar parameterization, we mention the following relations between (h n+ 2,h n+1 ) and (a,h n ), which follow from (U.2) and (2.7), ■ n+2 + 1 = T 2 h n -t„ h n 1 - a + a/c )h (1 a/c )h It can be shown, Dahlquist (1979), that the only possible G-matrices for the pair (p n ,o n ) de- fined by (H.1) are of the form (apart from a posi- tive multiplier), (It. 5) G = g 1+g-2c 1-g 1-g 1+g+2c ac) 2 < Uabc triable stepsize method is characterized by the triplet (a(o), b(a), c(a)). We want G to be the 'or all a, apart from a scalar multiplier. .lows that g and c are independent of a, hence by (b.1), p n has to be the same all the time. jne that c * 0. If we require b(a) to be continu- en, by ( U . 3 ) a(a)/a is continuous at >,ce a(a)>0 everywhere, we conclude that = a' =0, and hence a(0) = 0, b(0) = 1 - c 2 . ,g=1, and hence the only possi- '2(1-c) 2(1+c) ■ ]■ < c < 1 raa earlier obtained by Nevanlinna I by the It is easily verified thai that. . with g= 1. (I4.9) (1 - c 2 -b(ct) - = Ua(o)h a) ■ , and we summarize our n PROPOSITION 4. For any constant c, < c < 1 , equations (^.8) and (*t.l) define a second-order accurate one-leg met- , ;h is A-contract i with respect to the G-norm (k.J). The parame + a and h are given by (h.i> 1. There is a geometric idea behind the method (I4.8). In the plane with cartesian coor-i (ac,b), eqn. ( U . 3 ) defines a one-parameter family of straight lines, the envelope of which is given by (h.8). Then (h.9) shows that this envelope is exactly the parabolic boundary of the domain de- fined by the second inequality (lt.6), for g= '. This coincidence, which we cannot yet "explain", indicates that ( U . 8 ) gives a complete account of methods that satisfy the restrictions formulated in this section. (Alternative methods are perhaps obtained by dropping the continuity requirement for b(oe).) It is remarkable that the BDF method, for which c = 2 when a = , is excluded. Nevanlinna (1979) has somewhat earlier than we found the same family (although with a very different parametrization) and shown that it is A-contractive with respect to the maximumnorm , for arbitrary stepsize sequences. His result is more flexible than ours, since it allows p n (?) to be changed during the computations, which we cannot do. On the other hand our result covers (at the time of writing) a different class of non-linear problems. The analysis can easily be extended to the case where the stability region is assumed to include a circle tangent to the imaginary axis at the origin. We have done some calculations concerning the unique stable method with k = 2, p= 3 and some explicit methods with p = k = 2 , but we must say like Fermat: the paper isn't large enough. REFERENCES, Butcher J. C. (1973), The order of numerical methods for ordinary differential equations, Math.Comp.27, 793-806. Dahlquist G.(l975a), On stability and error analysis for stiff non-linear problems, Dept.Comp.Sci. , Roy . Inst . of Techn. , Stockholm, Report TRITA-NA-7508. Dahlquist G. (1975b), Error analysis for a class of methods for stiff non-linear initial value problems, Numerical Analysis, Dundee, Springer Lecture Notes in Mathematics, 506, 60-7^. Dahlquist G. ( 1 978 ) , G-stability is equivalent to A-stability, BIT 18, 38U-U01 . Dahlquist G.(1979), Some contractivity questions for one-leg- and linear multistep methods, to appear as Dept.Comp.Sci., Roy . Inst . of Techn. , Stockholm, Report TRITA-NA-7905. Linicjer W.(1978), personal communication. Nevanlinna 0. and Liniger 11.(1978,79), Contractive methods for stiff differential equations , Part I, BIT 18, h'/(-h-{h; Part II, BIT 19, 53-72. Nevanlinna 0.(1979), personal communication. Skeel R.D.(i I , mivalent forms of multistep formulas, Dept.Comp.Sci., Univ. of Illinois at R, P or1 UIUCDCS-R-78-9 1 *0. 1-4 D-STABILITY M.van Veldhuizen, Wlskundlg Seminarium Vrije Universiteit , De Boelelaan 1081. 1007MC AriSTERDAM In this paper a new stability criterion for discretization methods for stiff ordinary differential equations is proposed. This new stability concept takes into account the variation in time of the eigen- vectors of the Jacobian matrix of the system. As such the new criterion complements the existing criteria; it does not replace them. The new criterion is applied to multistep methods, a generalized Runge-Kutta method. and Galerkin methods ( with numerical quadrature) .We refer to the new stability criterion as D-stability because of the relationship with the decomposition error. In the definition of D-stability we consider linear homogeneous systems of the type (1) x' * A(t)x, where fait) is a matrix and x' ■ dx/dt. The matrix Att) may be obtained as the Jacobian matrix [evaluated along a smooth solution) related to a non-linear system. We restrict ourselves to discretization methods which discretize the system (1) by the recursion (2) y . , = G(t. ,h)y. *i+1 l *i wt-ere h is the stepsize.t =t +ih, and where G(t.,h) is a matrix. If the discretization method is a one-step l i method, then y is the approximation to x(t,). If the discretization is a multistep method, then the recursion (2) is the formal one-step recursion corresponding to the multistep recursion and the meaning of y.has to be changed accordingly. For stiff systems it is very difficult to obtain stability for practical values of the stepsize h. In order to obtain valuable information about discretization methods for stiff systems for realistic values of the stepsize G. G.Dahlquist restricted himself to the scalar model equation x' = Ax by introducing his concept of A-stability ,cf .[2]. This concept has been modified in various ways, but in all modifications the time- dependence of the eigenvectors of the Jacobian matrix is neglected. However, the time-dependence of the eigenvectors of the Jacobian matrix cannot always be neglected in the analysis of stiff discretized problems. The trapezoidal rule provides an example. For this method the source of the difficulties is the bad decomposition error. As such the difficulty was known to Dahlquist and Lindbergh] and Van Veldhuizen [7] .H'^wever, no simple criterion for investigating such phenomena in other methods seems to exist. In this paper we propose the concept of D-stability. We define it in such a way that a D-stable method does not suffer from a bad decomposition error for a large class of stiff problems. Definition . The discretization method resulting in the recursion (2) is called D-stable if for all t in the interval, all h £ ^°' n ^ and a11 stiff systems of a class 3 ( see below) ||G(t,h)|| < M, M a constant independent of h and the system in the class o< Consequently , D-stability is a weak concept as only the transition in one step should be bounded uniformly 2-1 ir the stepsize h. On the other hand, D-stability is a strong concept by requiring this uniformly in a class of stiff systems. Obviously, the concept of D-stability complements, but not replaces, the concept Df --stability and related ones. Now we come to the description of the class o of stiff systems. This class jj s-ould be sufficiently large because otherwise the concept op D-stability becomes a triviality. On the other hand, we aim at a simple criterion and so we must be restrictive in the choice of the class O . The class J as described below seems adequate. Class . The class S consists of all linear systems of the type x* = A[t)x satisfying the conditions: (51] x(t) e £ 2 . (52) A(t) = T(t)A(t)T(tr 1 for all t. Here A(t) is a diagonal matrix, A(t) = diag{ A(t)/e, ptt] }. Also ReC Alt)) < \ < 0,and e £(0,e Q ]. - 1 (53) X, u, T and T depend smoothly on t and e, and the derivatives from order zero up to a certain order ( sufficiently high ) are bounded uniformly in t and e e ( , £_]. Obviously , the systems in the class o are stiff. In many cases it is realistic to consider a subset of the class J . E.g. one might consider systems with real eigenvalues only. Also, it is possible to restrict the class o to the set of stiff systems with a particular type of coupling (strong, weak] between the smooth component and the transient one. In this survey we do not go into these details. The stiffness is introduced in the class O by the eigenvalue X(t)/e. We must have e e (0,£ ] in the class O because otherwise the definition of D-stability becomes a trivial one. We now describe some of the results;we look at three classes of discretization methods. Multistep Methods . A multistep method defined by the generating polynomials p and a of degree £ k i6 D-stable if and only if cr(z) » B n z . This result once more justifies the choice of the backward differentiation methods in the code by Gear[5]. However, the one-leg method (see Dahlquist[3] ) corresponding to the multistep method is D-stable if the method is A -stable. This result emphasizes the importance of one-leg methods as introduced by Dahlquist[3] . leralized Runge-Kutta Methods . Verwer[8] introduces the concept of internal S-stability. Roughly speaking, a generalized Runge-Kutta method is internally S-stable if it has good damping properties at each stage when applied to the model x' = Xx. It turns out that the internally S-stable methods of Verwer[8] are not D-stable. This shows that D-stability is not a consequence of good damping properties. lethods . Here we restrict ourselves to Galerkin methods with approximate quadrature as considered among others ) by Axelsson[l] and Weiss[9]. With suitable Gaussian quadrature or Radau quadrature these methods reduce to collocation methods, while the corresponding method with Lobatto quadrature is not a collocation method. We have the following result: all A-9table collocation methods are O-stablejthe method 1 Lobatto points ia not D-stable. This result is important in view of the use of such methods in 2-2 REFERENCES [1] Axelsson.C: A Class of A-stable Methods. Nordisk Tidskr. Inf ormationsbehandling (BIT) 9,185-199(1969) [2] Dahlquist.G.G.: A Special Stability Problem for Linear Multistep Methods .Nordisk Tidskr. Informations- oehandling (BIT) 3.27-43(1963) [3] Dahlquist.G.G.: Error Analysis for a class of Methods for Stiff Non-linear Initial Value Problems. Lecture Notes in Mathematics. vol . 506. pp. 60-72. Berlin-Heidelberg-New York:Springer 1976 ["♦] Dahlquist.G.G. .Lindberg.B. : On some Implicit One-Step Methods for Stiff Differential Equations. TRITA-NA-7302,0ept of Information Processing , Royal Institute of Technology .Stockholm. [5] Gear.C.W.: Algorithm 407. DIFSUB for solution of ordinary differential equations ,C. A. CM. , 14, 1 65-190 (1971) [6] Russell ,R.D. .Christiansen, J . : Adaptive Mesh Selection Strategies for Solving Boundary Value Problems. Siam J .Numer.Anal. (to appear) [7] Veldhuizen.M. van: Convergence of One-Step Discretizations for Stiff Differential Equations (thesis). Mathematical Institute, University of Utrecht .Netherlands (1973) [8] Verwer.J.G. : S-stability properties for generalized Runge-Kutta Methods. Numer. Math. 27, 359-370(1977) [9] Weiss. R.: The Application of Implicit Runge-Kutta methods to Boundary Value Problems. Math.Comput.28. 449-464(1974) 2-3 Stability of explicit time discretizations for solving initial value problems Rolf Jeltsch, Institute of Mathematics, Ruhr-University Bochum, D-4630 Bochum, Fed. Rep. of Germany Olavi itevanlinna, Department of Mathematics, Oulu University, SF-90101 Oulu 10, Finland stability regions of explicit "linear" time dis- cretization methods for solving initial value problems are treated. A discretization is called "linear" if tne right nand side of the differen- tial equation, its total derivatives and iterates of these occur only linearly. Tnus the class of methods considered includes kunge-Kutta-methods, linear multistep methods, multistep methods using higher derivatives, predictor-corrector methods and cyclic methods. If an integration metnod needs m function evaluation per time step, then we scale the stability region by dividing by m, and show that for any two methods, satis- fying some reasonable conditions, the boundaries of the scaled stability regions necessarily intersect. Using this comparison result bounds for the size of the stability regions for three different purposes are given: I) for "general" nonlinear ordinary differential systems 2) for systems obtained from parabolic problems and 3) for systems obtained from hyperbolic problems. First we consider the size of disks one can have inside the stability region with origin on t their boundary. In the scaled sense the largest disk is obtained with the Luler's method and for any smaller disk there exist explicit linear k-step methods for order p = k which are stable in such a disk. A-stable methods with high order. For hyperbolic problems one is interested in methods with long stability intervals along the imaginary axis. It is shown that in the scaled sense, the leap-frog method has the largest inter- val and that for any shorter interval kS {2,3,4} there are explicit linear k-step methods which areJ stable on that interval and of order p = k. however for k = 1 mod 4 an explicit linear k-step method with p = k has no nontrivial interval of stability] on the imaginary axis. By the previous result about the largest stability:! disk with the origin on its boundary it follows that the radius of such a disk is bounded by m fori predictor-corrector methods, where the corrector is iterated m - 1 times and the function is evaluated m times. It is shown that this disk can j be obtained only if one uses different correctors j at each iteration step. When one, as is customary,: iterates a fixed corrector, then as m tends to infinity the maximal radius does not tend to in- finity, and in fact it is shown that in the limit ; the maximal radius is 2. For parabolic problems one is interested in methods which are stable for long intervals of the form | -a,o] , a > o . It is shown that for explicit linear multistep methods of order p = k lac length of the interval is always bounded by 2, but already for p ■ k - 1 one can construct methods which are stable in bounded sectorial sets for arbitrarily big radius and any opening less than n . hence, although an explicit method cannot be A-stable we can have approximately 3-1 STABILITY PROPERTIES FOR A GENERAL CLASS OF METHODS FOR ORDINARY DIFFERENTIAL EQUATIONS J.C. Butcher, University of Auckland, Auckland, New Zealand Abstract The concepts of G-stability for linear multistep methods and B-stability for Runge-Kutta methods are combined in a unified approach to non-linear stability for a class of methods general enough to include these as special cases. Recent work on these developments, by Kevin Burrage and the author, is reviewed and, to exemplify the new ideas and techniques, an algebraic stability analysis is presented for the backward differentiation methods of orders 1, 2 and 3 • Introduction The property of A-stability, long accepted as a desirable condition for methods to possess if they are to be used with stiff problems, gives only limited guidance to the numerical behaviour if these problems are also non-linear. For this reason, there has grown up an interest in non-linear problems of the form y'(x) = f(y(x)), where y has values in an inner-product space X, which are such that < 0, for all u,v 6 X . For such a problem, an immediate consequence is that for two particular solutions y and z, By(x)-z(x) is a non-increasing function of x . It is interesting to know whether a particular numerical method exhibits similar behaviour in the results it generates. This is simply expressed for one step methods as the requirement that for two solution sequences . . . ,y .. ,y ,... and ...,z ,,z ,..., it holds that • n-1* n* ' H y n" z n" - ,ly n-l~ z n-l"' For' multistep methods, in which r vectors in X are used to represent the solution information at the end of a step, it is necessary to construct a norm on X using the tensor product of the inner- product on X that comes with the problem and a new inner-product on 1R that depends on the method. Bounds on the behaviour from step to step can then be expressed in terms of the norm 4-1 formed froa the tensor product. tf G La the matrix of the Lnner- producl i n IK , which i a character- Lei Lc of the method | the corresponding stability propci t \ i -. known as i itabilitj .Hid \>a^ firsl Used by G. Dahlquisl S . The corresponding propertj ol Runge-Kutta methods was i i i — t proposed i>\ the present author [4l .uui , i ii collaboration with K<-\ i n Burrage, waa further developed [1] and combined with tbe ideas of G-stabilit> to g i \ e the algebraic stability property for genera] linear methods [2]. This paper i-^ pai tly .i re\ Lew of the joint work by Di . Bui i age and the author, and partly a study j t'rom the viewpoint of algebraic Stability , of three special methods. Be methods, the backward differen- t Lai ion methods of orders 1, 2 and 3, are important components of many algorithms in present use for solving stiff problems. Part itioned general linear methods . Let \. , be an sxr matrix, A an B 1 . an sxs matrix and rvr matriXj B an rvs matrix. A partitioned general linear method is a general linear method (or (A,B) method) [3], such that V and B are each (s+r)x(s+r) matrices partitioned as follows A 1 2 A 22 B = B ll ° B 21 ° where, of course, the various matrices are not necessarily the same si/es. We will, for convenience, denote such a method by writing the non- trivial parts of its coefficient matrices together in the format 12 22 15 11 21 For such a met hod, when USed to solve the equation \'(\) i(\(\)) with step size Ii, the a+r vectors which represent al] the quantities computed in step n, saj (n) „(n) Y (n) (n) (n) (n ) ] '2 "••'\s ' y i ' y 2 '•••'', are related to t he previous step by (2) Y (n) = J a (H) (n-1) + j=l ^ h J b (ll) f(Y (n) ), i=l,2,...,s, j=l iJ J (n) ' (22) (n-1) y = V a y ' + i j=l ij j h rb (2l) f(Y (n) ), i=l,2,...,r, j=l ij J . (12) (22) .(11) .(21) where a v , a . b , b v are (3) 1J !J typical elements in ij A 12> A 22> B ll 21 respectively. Y (n\ Y (n) .(n) 12 s can be interpreted as s stages, as in a Runge-Kutta method. We see from (2) and (3) that no use in later steps is made of these quantities once the step is completed. On the other hand, (n) (n) (n) . +u y ,y ,...,y represent the r 12 r quantities needed at the end of a step to allow the computation of the next step. Although we have used the symbol y for one of these pieces of i information, it is not necessarily an approximation to the solution value near the n integration point x . Rather, it should be thought of as an approximation to u. v(x )+hv. y (x ) where u.,v. are numbers characteristic of the method. The special case r=l , with A equal to 1, gives s stage Runge-Kutta methods, whereas 4-2 12 22 11 21 *1 2 • °k Pi H Pk Po °1 o 2 .. • °k h p 2 Pk Po 1 . . . ... • 1 . • . . . . . 1 . . . 1 . . . 1 . . . • • • • 1 • • •- • • • • • . . . .1 gives the k-step method y n = °l y n-l + a 2 y n-2 + '" +a k y n-k + h (3 f(y n ) + 3 1 f(y n _ 1 )+...+3 k f(y n _ k )) , l • l v(n) (n) xn which Y v = y = y„ , 1 1 n (n) .(n) 'n-l : = y. n-k+l , 2 y k + l= hf ( y n>> y ( Zl= M ^n-1>- ( n ) LC/ y = hf(y 2k n-k+l k+2 ). A non-linear test problem . Let (X, <•,•>) denote a real pseudo inner- product space. We recall that an inner- product <', #> on X is a bilinear functional such that = for u,v 6 X and such that for non-zero u € X, > 0. In a pseudo inner- product this inequality is weakened to -ii,u - > . Let II • H be the corresponding pseudo norm defined by Hull = 1 ^ 1 for all u e X . Let f : X -. X be continuous and such that, lor all u 6 X 2 where 5 is a real number. A consequence of (4) is that any solution y to the differential equation y (x) = f(y(x)), is bounded in pseudo norm by a multiple of exp(i'x). To verify this, we need only differentiate Hexp(-?x)y(x)|| 2 = exp(-2?x) to obtain a zero upper bound on the derivative. This rather artificial test problem can be used to study the behaviour of the difference of two solutions, as in the G-stability theory, even when the differential equation under study is non-autonomous. If G is a non-zero non-negative symmetric rxr matrix, then we can define a pseudo inner-product on X using the elements g. . (i,j=l,2, ,r) of G, and an existing pseudo inner- product on X. Denoting the pseudo inner-product on X define i v = (v 1 ,v 2 ,... G by <-,-> G , define its value for u = (u, ,u„, > v r ) b y > u r )» i,J=l g . . (4) defined from <•,•>„, as the measure of r • • the size of a vector in X with a view to obtaining bounds on the behaviour of /(n)(n) (nK , . ,-,-> (y ,y ,...,y ), computed using (1), 12 r in terms of the behaviour at the previous step. We will attempt to find bounds that correspond in some way to bounds on the behaviour of the solution to the underlying problem. Algebraic stability and its general - izations . A partitioned general linear method (l) is said to be (k,/.) algebraically stable for k a positive number and I an arbitrary number, if there is a non-zero non-negative symmetri; matrix G, and a non-negative diagonal matrix D, such that the (symmetric) matrix 4-3 M I T 4 2 D ~ A 22 GB 2] " 2 * A 12 DA 12 • 2 '' A 12 DB 1] DA I2~ B 21 GA 22 DB 11 +B 11 D ~ B 21 GB 2] T T - 2tB ll DA 12 " MB 11 DB 13 i -^ non-negat Lve, It is a principal result in Burrage and Butcher [2], that if a method satisfies this condition and if a problem satisfies (4) with > = h » , then (> 1 ,y 2 (n)x ,y ) can be (pseudo-) estimated by 'G = ftfe are interested in obtaining for a method a (k,-?,) algebraic stability result in which k and i, are functionally related by to(t) = k ' . To the extent that tp is a good approximation to the exponential function, we obtain behaviour for the computed results that is a good model for the behaviour of the exact result. In particular, if a method is (1,0) algebraically stable it will be called, simply, algebraically stable. In addition to this non-linear generaliz- ation of A-stability, we are also interested in knowing whether, given k as an arbitrarily small positive number, (k,/, ) algebraic stability can be obtained for sufficiently small I (generalizing the "behaviour at infinity" aspect of L-stability ) , and whether we can obtain (l,i.) stability for sufficiently small I (generalizing part of the definition of stiff stability). Finally, Ln t h i s section, we note that, Lf Ik , Ls non-singular ( as Ls typicallj the case Ln methods suitable for s1 iii problems), .i congruence trans formal Lon can be app] i cd to M given by (S), to yield a symmetric matrix <>i i simpler form for many methods . Ii -1 V,, = V,., - ^oi B ll 12 a,icl ' ' '*' transforming matrix is I " B l' A 12 T then M = P MP is given by ;6) M = kG-A 22 GA 22 T „ -A 12 D A 22 GB 21 ■ DA 12~ B 21 GA 22 DB 11 +B 11 D - B 21 GB 21 - 2£B n DB ii Where M takes on a simpler form than M, we would use its non-negativity as a criterion for (k,/, ) algebraic stability. Application to backward differentiation methods . The backward differentiation methods of orders 1, 2 and 3 respectively, given by (7) y = y , + hf(y ), J n J n-1 J n ' (8) y n = 4 y n-l " Tn-2 + l hf(y n } ' (9) y n = TTY^!- ify n -2 + IT y n-3 + TT hf l , or, what is equivalent, k = 1/(1-2/,) , t<0 , kg 11 -1 kg ll- g 22 -1 2-g i:L -24 - 4 -4 '12 3 3^12 kg 12 •3 3 e 12 kg 1 3 22 1 3 !— i 3 9 fe H 9 7tl and kg ll" g 22 kg 12" g 23 kg 13 kg 12 -g 23 kg 22 -g 33 kg 23 kg 13 kg 23 kg 33 -18 _6 _9 _6 ll"ll g 12 ll"ll g 13 -2 11 JJ3 _6 ll"ll g 12 ll"ll g 13 22 11 11 121 8 11 — 72i 121' Not surprisingly, determination of which (k,/,) pairs give (k,/,) algebraic stability, is trivial for the method niven by (7) and progressively more complicated for (8) and (9)- Ih f;orem 1 . The first order backward differentiation method is ( k , J, ) il '^nbraically stable whenever k 1/2 = i/(2-(i+2t) 1/2 ) , t>q . Proof . We will show that, with t related to k as in (10) and (11), g ll ,g 12 ,g 22 may be cnosen so tnat M is of rank 1 with non-negative diagonal elements , Let a = kg 12> 3 = kg 22 so that M being of rank 1 would require that g 1T = a 2 /pk + p/k 2 , g 12 = a/k , grt^ = P/k (implying that G is positive if 3^0) and furthermore that (12) a + 3 (4+2 a/k) = , (13) (12-8-03 = 4a 2 /k + 43 2 'k + 1. Since the values of I given in (10) and (11 ) satisfy the inequality 12 - 8.1 > 0, any real solution to (12) and (13) will necessarily have 3>0 (because the right hand side of (13) is positive). Hence, it remains to show that, with the given value of t, (12) and (13) do in fact have a real solution pair. If 3 i s eliminated, the resulting equation for a can be written in the form (14) 2(a+k) 4 + (a+k) 2 (k-k 2 -2/,k 2 ) + (k 3 -k 4 +2*,k 4 ) = Pi 00 f . < hoosf ,-1/2 ; n If k < 1 , the value of t given by (10) is satisfied by a + k = 0. If 1 I, the value of t given by (11 ) is 4-5 ~»«n h aa bo sake (\A) •> perfect square .iiui the negativity of the coefficient Di ia+k)" guarantees the existence oi a > posit i \ <• root for ( a+k ) " . theorem Hie third order back- ward different i.it Lon method is ik,') algebraically stable for (k,t ) pel .u ed so I tut (15) 4 17 U li 4k -1^ -3)3/2 24 V k JJ : °<"S# do) t-^+*| ^ 2 (3K + 1), k>$ Proof i ^ s in the proof of theorem 2, we \> i 1 1 ->lu>\> that G may be chosen so that M is of rank 1 with non-negative diagonals. Let a = kg-i-i» P = kg,!) ■y = ku ■ , , so that •33 G = -2- + £— + -I k - k 2 v k^ °£ + _J i: k-y k a i L44ft-108k ) , J ' I l 44k ■ ■ 63k 3 )-y ■ < l 6k * 6+1 2k*) = 0, ( 21 ) (12965 2 -1944k6+720k 2 +144k ) , Z -16mk 4 Y+(l44k 3 5 2 +216k 4 fi+8lk ! '+16k 4 )=0 Eliminating y from (20) and (21 ) Leads to an equation that can bo writ ten as (22) (a6 4 +bfi 2 +c) 2 - df, 2 = . where a = 20736, b = 864k(8-3k), c « 192k 4 m - 63k 3 (l6+8lk), d = 65536k 4 (km+27) 2 - 2359296k 3 (l+9k) 2 . What we will show is that the value of n resulting from -l chosen by (IS) and (16) leads to a repeated non-negative 2 root for 6 in this equation. In the case k > 64/49, the choice of I implies that d=0 and that c = 9 k 5/2 (4-9k l/2 ) 2 (8-7k l/2 ) < . It turns out V at k < 64/49 would 2 not allow a positive root for 6 to exist if d=0 so we consider the condition that (22) will have repeated roots if d 7^ 0. This condition is and, writing m = 3 3— 1§ -t- , the conditions on the last row and column of M for it to be rank 1 are (17) 9k 2 Y + 3kap + 3Py - ok 2 = , H-) 9ky + 2k p - 6av = (10) k 3 + 9k 2 a 2 + 9kp 2 + 9y 2 - mk 3 y = , and we note that, since the values of I siven by (IS) and (16) imply that m>0, the positivity of y, if real solution trios (a,P,y) to (17), (18 ) and (19) are shown to exist, will be assured. It 3 is convenient to substitute a = -rk + 0 Since 27k 3 /2 + 54k 1//2 - 8(l+9k) = (3k l/2 -4)(3k l/2 -(2V"2))(3k 1/2 -(2-V"2)), (26) is clearly satisfied for k £ [(2i/2")/3,(2+/2)/3] and, since the left hand side of (26) vanishes when k=0 or k = 64/49, it will be sufficient to prove that the derivative of this expression is positive when 1 In k 6 (0,(2-V"2)/3) and negative when k l/2 £ ((2V2)/3, 8/7). The value of this derivative is (27) ^L k 2 + 27k 2 - 72 - -|(4-3k) and, when k 1 ' 2 € (0, ( 2-Jl ) /3 ) , this is bounded below by 81 2-£L + 27 _3_ 2^2 1 72 - |(4-0) 2 > because -^k 1 ' 2 + 27k" 1 / 2 is easily verified to be an increasing, and - -^•(4-3'<) ' a decreasing, function in the interval in question. For k 1 ' 2 6 ((2V2)/3,8/7), ^k l/2 + 27k" l/2 is an increasing function and we bound (27) by substituting k = 64/49 to obtain -27/8 < . A final remark: the backward differ- entiation method with order 3 is (1,- y^) algebraically stable and it is noted with Lnterest that the abscissa - yr exactly eea with what is obtained in a linear Pf— stability) anal ys i , . REFERENCES [l] Kevin Burrage and J.C. Butcher, Stability criteria for implicit Runge-Kutta methods, SIAM J. Numer. Anal, (to appear). [2] Kevin Burrage and J.C. Butcher, Non-linear stability of a general class of differential equation methods, Computational Mathematics Report No. l8, (1979), University of Auckland. [3] J.C. Butcher, On the convergence of numerical solutions to ordinary differential equations, Math. Comp 20 (1966), 1-10. [4] J.C. Butcher, A stability property of implicit Runge-Kutta methods, BIT 15 (1975), 358-361. [5] G. Dahlquist, On stability and error analysis for stiff non- linear problems, I., Depart- ment of Information Processing Report No. NA 75.08, (1975), Royal Institute of Technology, Stockholm. 4-7 A-STABILITY OF ONE- AND MULTISTEP METHODS nerhard WANNER Abstract. This work with S.P.N«$rsett and E. Hairer deals with a beautiful and reo- metric discussion of A-accertability pro- perties of stability functions arising f~om one- and multi-sten methods. We "ive many necessary or sufficient con- ditions for A-stability and nrove the con- jecture of Ehle (1969) as well as the con- jecture of Daniel and Moore (1970). Properties of Order Stars. Consider a rational stability function of a one-step method M . „, ,. P(z) de^(P)=a (nb. of zeros) nvz; " OTzT dee(Q)=a. (nb. of poles) with the order condition ,~ s z Bf N n + 1 n+2^-, p + 3\ (2) e -R(z)=e ..z- + e ._z +0(z f J ) . P+i r + ^ R is said to be A-acceptable (or the method is A-stable ) if|R(z)(;i UU111I1I mir.M'.i ■ 1 ■ -. * 1 : tIMUUII . 1 1 : 1 1 11.11111:1 1 1 • 1 1 1 : 1 ; 1 ■. 1 1 1 \ • UtltlUII III tin 1 1 •■.::*. 1 * 1 ■»«*H1 The form of these stability domains is clumsy and not very instructive. So the first treatments on A-acceptability have been loaded with difficult computations. In contrary, the set z (')) A ={z; R(z) } for the same rational approximation is sket- ched in Fip;.2 and pives much more insicht (like a hand under an X-ray skreen) . At every point where R(z) has a contact of order p with exp(z), the set A consists of p+1 regularly distributed finders, so we call it the order star of R. Fifr. Example of an order star, F1r;.1. Example of stability domain. The following properties are conseouences of simple complex analysis: [ Lemma 1 . On the imaginary axis, A and D are complementary. So we can have A-acceptability only if A does not meet ilR . Reason: Because of |e ^|=1 for y real. [ Lemma 2 . ^or z+0, A behaves like a star with p+1 "finders" of equal width w/(p+l), sepa- rated by p+1 similar "white finrers" of the complement of A. Reason: r or z+0, z€A is equivalent to |R(z)e- z | = h - e- z (e n + 1 z n + l +e n + 2 z^ + 2 +..)|>l 5-1 or Re(e p+lZ P +1 + (e p+2 -e D+1 )z^ +2 + ..)<0. Thus the sign of the error constant deter- mines the color of the finders. f Lemma 3. For z*«m, A tends to the negative half-plane with exactly two boundaries leroing to infinity. *~ ^ z Reason: elementary properties of e . Lemma 4. Each bounded finder of A contains at least one pole of R; each bounded white finger contains at least one zero. If several bounded fineers collapse, then they contain together at least as many poles or zeros as there are fingers. This lemma is less trivial, see flf . While Lemmas 1 and 3 are linked to special properties of the exponential function, Lemmas 2 and 4 remain true for any order stars of two functions R(z) and S(z). 2. A-acceptability of Pade approximations. For Pade approximations we have p = a+in and Lemmas 2 to h do not permit any other con- figuration as u black fingers (bounded) to the right and o white bounded fingers to the left. Because of Lemma 1 and 4, we have A-acceptability iff the imaginary axis passes behind the last bounded black finger (as indicated by dots in Fig. 2). This is equivalent to o22w-2, !R(-)|<1, and the coefficients of Q(z) have alternatine signs. Then R is A-acceptable. A similar version of this has been proved differently by Crouzeix-Ruamps (1977). Theorems of Ehle's Thesis are special cases. Theorem 11. Suppose p2.2w-3, |R(iy)|roof. At least [(p + 1 ) /2] fingers of A start in C . They cannot cross iff! and must there- fore be bounded. Since in total we have u> poles, it follows that [(p+1) /2] . a [ Theorem 7. If R is A-acceptable and p=2w, then for the error constant sign(e p+1 ) = ( -1 ) Proof. Look at the color of the fingers. [ theorem 8. If R is A-acceptable and p = 2a), "•" ivi Ns.)i&iin ■ Proof. Consider the order star R where R(z) is "compared" to the diagonal approximation R (z). The following are inverse results assuring A-acceptability: y real, and the coefficients of 0(z) have alternating signs. Then R is A-accentable . "Proof. Similar as for 9 or 10. see 121 . Following are the exact characterizations of rational approximations of order 2u-3 and 20)-^: All rational approximations of order 2u>-k can be written as ,(i) (5) R(z) = I i : o M (i) (Dz u - 1 /I i : o N^ i; (0)z- where N(t)=P u (t) + A 1 P u _ 1 (t) + .. + A k P w _ k (t), P (t) = (^ r ) 0, ii) A a +A 3 <0, iii) A 3 >0. 3 is T heorem 13. R(z) of (5) of order 2ui-4 A.+A,<0 is A-accentable iff i) 1+A +A,>0, ii) 2 "4- ' x " ; "1 3" iii) (2u-5)A 1 A i( + (2ui-l)A 3 >0, iv) A^O. The proofs involve Theorems 9 and 11 above and some continuity arguments of the poles of R(z) as functions of A 1 ,A 2 ,A^ ,A^ . See /6/ for details. The following theorem gives an order bound for approximations with multiple poles or zeros: Theorem l'l. Suppose that R possesses only a ( different zeros. and u o different poles. Then p\ . Extension to multi-step methods. Eor a k-step method the rational function (1) is replaced by a characteristic polynom- (1') 0(z,R)=Q o (z)R k +Q 1 (z)R k " 1 + ..+Q k (z)=0 (0)*0, deg(Q k )=a, deg(Q )=w. R(z) becomes an algebraic function to be con- sidered on its Riemann surface M. For linear multistep methods we have a F • "*l r * r ** •£.«*«•« Order function R^z), the order condition can now mqS! w-a-^oo V theorems ' RTT 1 * he *rit1 • - (1978), Wtf? z R #_x _ p + 1 + nf 7 n+2 1 /2/ H. Wanner, E.Hairer.R. P. Nrfrsett, When n \ KZ > ■ e p + l z ; ' I-stabillty implies A-stability, The above Lemma 2 now carries over to the PI"' 18 (1978), p 5a3 corresponding sheet of the Riemann surface, ,,. _ „ .„ ,, *,^ t i •■ . -, tl _ . when we now consider the order star as a /V -"ai er Hncondl onallv stable methods subset of M. Also for Lemmas 1,3, and U ^"^^'r ^/Th" \l L e °" at l0nS ' . , , j mw to anpear in Numer. M ath. ,1979 exist natural extensions. Thus: '' 7 ' 7 ^v,^„o~c £ 7 P or,^ 11 r^^moir, ^r,..^ i „ fv,. / '' / S . P . Fdr se t t , . Wanner , The real-nole Theorems o,7,c, and 11 remain true in the _ ... . - ' . . ' . multl-ster case sandwich for rational approximations and oscillation enuations, appears in PTT Theorems 6 and 8 have been known as the 1979. "conjecture of Daniel and Moore", the special .... „ _ „. .. _ ,. n case u = l is Dahlouisfs classical Theorem /5/ S.P.^rsettO. Wanner Perturbed collo- that p = 2 is the highest attainable order of cation and Runre-Kutta methods, submitted an A-stable linear multistep method. zo lumer - atn - Theorem 9 needs the following modification: /6/ Y .Burrapre ,C .Wanner , Characterization of Theorem 9'. Suppose that p>2u>-l and _ _ , . _ R(ly) ! to he writt en. jble La extends the result of Jeltsch (PIT 1978) to u>l. Theorems 5 and 1*J do not extend as they are to the multistep case, since the Proofs use the fact that C Is simply connected, which is no longer true for arbitrary Riemann surfaces . 5. Extension to second order enuations. If a stability test for methods for y"=f(t,y) is made at y"+By=0 with P real, the characteristic polynomial becomes (1") Q(z,R)=Q o (z 2 )R k +..+Q k (z 2 )=0, where now R=l becomes a double zero of Q(0,P)=0. But since there appears z2 in (1") only, there is no monodromy around this point on the corresponding Riemann sheets and these two zeros continue as holo- morohlc functions in the neighbourhood of the origin. One of these two solutions then satisfies the order condition (2') and there is a^ain a beautiful star on the corresponding sheet of the Riemann sur- face. So Theorems 6, 7, 8, and 11 remain true also in this case, The special case u=l was proved by Dahlquist (BIT 1978). Theorem 6 has the following corollary: Corollary 6 ' . If the solution of 0(z,R)=0 Tsee Lambert -V/at son , J.Inst. Math. Applies. 1976: this is eouivalent with the fact that all roots of 0(iy,R)=0, y real, satisfy |R|=1), then p<2 u . See Hairer /3/ for details. Hairer has con- structed unconditionally stable methods of higher order for this type of enuations. A-stable and P-stable Runee-Kutta methods 5-3 ORDER AND STABILITY OF COLLOCATION METHODS WITH REAL EIGENVALUES. S.P. N0RSETT INSTITUTE OF NUMERICAL MATHEMATICS UNIVERSITY OF TRONDHEIM, TRONDHEIM-NTH The collocation method for obtaining an approximate solution to the initial value problem (1) y'=f(t,y), t€(a,b), y(a)=y Q , amounts to finding a polynomial u of degree at most m such that i) u(t ) = y ii) u'(t + c. h) = f (t +c. h,u(t +c. h)) , l i ' i i = 1 ,' • • ,m , iii) y l = u(tj = t Q +h) where c. , i = 1,««*,m are m distinct real numbers. By using a result of Alekseev and Grobner on the nonlinear Variation-of-constants formula for (1) a very short and nice proof can be given for the already known result that the collocation method possesses the same order as the corresponding quadrature formulae . When the collocation scheme is applied to the test equation y' = Ay, z = Ah, we get •j i = R(z)y Q = (P(z)/Q(z))y where ( m) P(z) = N(1)z +---+!r ' (1 ) , /.) = M(0)z m +---+!I (m) (0) and N(t) = (t-c. ) • • • (t-c )/m!. R(z) is an approximation to exp(z) of order s>m method has order . (N is called the N-polynomial to P. and is related to the C-polynomial p of N0rsett by N(t) = (-1 ) m p(1-t). ) It is well known that the collocation method is equivalent to a special subclass of Runge-Kutta methods with RK-matrix {a u = /£. (tMt} where i , i = 1 I. (l) J n (t-c, )/ ji (c. -c, ) k i k k*j k*j By writing Q(z) = ( 1 -y, z ) • • • ( 1 -y z) we have that y , ••*,Y are the eigenvalues of 1 w A and that if y. is a simple eigenvalue of A, a corresponding eigenvector x. of A is (x H' ' nj I J = N' (c. )y. +• .,( m) , . m + N (c, )y. These results are, as pointed out by Butcher, essential in connection with the practical implementations of implicit RK- methods. In particular we are interested in RK-methods where the eigenvalues are real and hopefully equal. Two questions arise : first which collocation points do we need in such a case and secondly what is the maximum order of the corre- sponding method. An upper bound on the order can be given by considering the order of rational approximations R(z) to exp(z) where R has real poles. For example, we can show by using the theory of order stars developed 6-1 by V. -.he following interesting result : Let =P k ( ), P k €TT k , Q. € Tij .as r < j complex different sn the order s satisfies s < if Q has no real zero k+r+1 if J has at least one real zero . This result was earlier obtained for the case r = by N0rsett & Wolfbrandt in a rather more technical way. The proof is now short and is a counting of bounded dual fingers versus unbounded dual fingers . The A-stability of these methods, and also of all one-step methods giving y, = R(z)y when applied to the test equation y 1 = Ay, reduces to the question of A-acceptability of R(z) = P fc (z)/Q. (z), i.e. |R(z)| < 1 for z £ C". For the optimal order approximation, s = k+j , R is an entry in the Pade-table. By intro- ducing the concept of order stars (see the talk by G. Wanner) one easily proves that the diagonal and the first two sub- diagonal entries are the only A- acceptable approximations. When Q has only real zeros y. = Y> I = 1,«»« t k, k = j, the approximations have the form (of order at least k) R(z ) = Q (-1 ) m L( k - m) (1/ Y )(Yz) m ) / (1-yz)" , approximations is as the following table shows k 12 3 5 V 1112 When y is a solution to L (1/y) = 0, Q £ v and L-acceptable approximations mav appear. Some results in that direc- tion will be discussed as well as the I-acceptability of the special approxi- mations R(z) P „ (z)/(1-yV ) k 2 k suitable in connection with integrating discretized hyperbolic equations. Some of the questions rised by Baker and Bramble who first introduced these functions will be dealt with. For further details see the following papers , |1| G. Wanner, E. Hairer, S.P. N^rsett; "Order stars and stability theorems", BIT 1 8 (1978), 475-489. |2| S.P. N^rsett, G. Wanner; "The real-pole sandwich for rational approximations and oscillation equations" . To appear BIT. |3| S.P. N0rsett, G. Wanner; "Perturbed collocation and Runge- Kutta methods". Submitted to Num. math. y e R . Optimal order k+1 is obtained for y a solution to L! (1/y) = 0. Let the k different Y - values be denoted by y with, < 1 /y < 1/y < < 1/Y... for each 1 1 '2 k of these Y-values we obtain the restricted-Pade-approximations , and an interesting property is that the order star corresponding to y = Y contains just one bounded finger of multiplicity k-v+1. Based on that result the A- acceptability of the restricted-Pade- 6-2 SAFE POINT METHODS FOR SEPARABLY STIFF SYSTEMS J D Lambert University of Dundee §1 Introduction The solution of the 2*2 linear constant coefficient system y' = Ay , y(x ) = y n , where A has real eigenvalues X. << X~ < and eigen- vectors CpCj is y(x) = k. (exp(X.x))c. + K 2 (exp(X 2 x)c 2 represented in the following diagram. X, real, *, « t=2 ^., m Vt < ° (1.1) P 1 [x,(K 1 +E ] )(exp(X ] x))c ] + (K 2 +e 2 )(exp(X 2 x))c 2 ] is a numerical approximation to the exact solution point P . Even if P f is close to P , the high gradient at Pj invokes instability if an explicit method is used. Correcting P. in the c.- direction was proposed by Lambert flj and developed by Alfeld and Lambert [2] into the CDS technique (Correction in the Dominant Space). The new proposal is to correct P. in the direction of the tangent to the local integral curve through Pj . For all K j , this tangent passes through S[x-I /Xj , (< 2 +e 2 ) (I-A-/X. ) (exp(X 2 x))c_] , which we call the safe point , since it invokes no high gradients. Consider the m * m linear system y' ■ Ay + g(x) , where A has eigenvectors {c f } and eigenvalues {X } such that If (1.1) is satisfied we say that the system, and the matrix A , are simply separably stiff . The exact solution at x = x is represented by the point P[x",k. (exp(X.x)c.+^(x) ] , and the numerical solution by P. Cx, (<• +e. ) (exp(X .x) )c.+tKx)+<5 (x) ] . The tangent at P. to the local integral curve now contains the safe point S[x-1 /X. ,ij)(x)+6(x)- 0J>' (x)+6' (x))/X. J , which, note, no longer lies in the subdominant space. Since 1/X. is small we could accept the solution at S as an approximation to y(x) ; this invokes a persistent error (cf CDS) of n'(X)/X, . More advantageously, we can accept the solution at S as an approximation to y(x-l/X.) , which invokes a persistent error of n"(x)/xf . §2 One-step Safe Point Methods Problem y' = A(x)y + g(x) =:f(x,y) , y(x Q ) = y Q , y,g £ K , eigenvalues {X (x)} satisfying (1.1) V x e [x ,x N ] . Procedure (1) Starting from x , take forward step, length n n i by conventional explicit method, the basic method . (2) Compute X. by the power method. (Note (i) the stiff er the better! (ii) estimates for initial iterate available from previous step; (iii) eigenvectors not needed (cf CDS).) (3) Proceed to safe point, and accept as an approxi- mi I mation to y(x +h +5 .,) , where 'n+1 ■1A,(VV This gives an irregular grid of data points, hence the motivation for one-step basic methods, which will yield numerical solutions on the grid (x J*n- X + jWV. )} • 7-1 Algorithm (SPM ) , z , ■ v + h £,(x ,v ;h ) n+1 n n f n'-'n* n '*l = V. + W ( W z n+l (i) (ii) )(iii) (2.1) Local TE p the order of (i) 0(h P+1 ) ♦ 0(d +1 ) , n n+ 1 Note: It is of course possible to start with the correction step (iii); this is equivalent to applying (2.1) to a perturbed form of the problem. §3 Other Safe Points Regard (2.1 (iii)) as an Euler step with £ as discretization parameter (cf singular perturba- tion?) . We could replace this by RK (any q-stage explicit Runge-Kutta method of order q , q S 4) . Clearly, from Fig 1, S can no longer be the safe point. It turns out that RK. has no safe point, but Y/X RK.. has a safe point given by ■1.596071638 . Using this in place of (ii) and replacing (iii) by an RK.. method, 3 (2.1) yields a method, which we label (SPM), , with local TE = 0(h p+ ) + 0(£ .) ; the persistent error can now be safely ignored, even for systems of modest stiffness. In the sequel we consider only correction steps of the form given in (2.1), but in every case replacement by a correction based on RK, is possible. §4 Stability and Sensitivity of (SPM ) . Definition Let A be simply separably stiff and have distinct eigenvalues {X | t =1,2, . . . , m} . Then a SPM is said to be dominantly stable if all solu- tions of the difference equation in {y } which result from applying the method to y 1 = Ay tend to zero as n ■* °° for all step lengths h such that Vt e *B ■ C = 2 « 3 ' , m , where J£ B is the region of absolute stability of the basic method. Theorem 3 (SPM) | and (SPM). are dominantly stable. Sensitivity Applying (2.1(i)) to y 1 = Ay yields z , ■ P(h A)y , where P(-) is a polynomial of n+ 1 n n degree p Dominant stability depends on the when t = 1 . If fact that (l+£A )P(h X ) ~ n ' we denote by X. the value of X. computed by the power method, then stability is retained only if |(l-A t /A )P(h X )| < 1 , which leads to the ,? , .. ,r ,. . . -d an excessive T n 1 ,-x,: requirement (X.-X.)/A ~ (hX.)" precision. Analogous difficulties beset CDS when the basic method is one-step. §5 Linear Multistep SPM In order to obtain a solution on a regular grid, we propose to correct by one forward (Euler) step to the safe point and then back by a further (Euler) step. This process retains the same safe point (eg Fig 1). We now use as basic method the K-step explicit LMM (p*,cr*) . Algorithm (SPM ). z n+k = (E k -p*(E))y n +ha*(E)f(x n ,y n ) (i) 'n+k ■!A,(x n +k) (ii) Vk = Z n + k + W^WW (iii) W = ? n + k " W (X n + k + WVk } (iv) • (5.1) Local TE = 0(h p+l ) + 0(£ 2 . ) n+k 3 (SPM), is similarly defined, replacing (iii) and (iv) by RK 3 steps. (SPM), is absolutely stable if the polynomial 7r(r;hA t ,£A t ):= r k -Q t [ r k - p *( r ) +h X t a*(r) ] (5.2) is Schur for t = 1 , 2, ..., m , where 2 Q t = 1 - (£X t ) . This is clearly so for t = 1 , while for t = 2, .... m , it = p* - hA a* + 0(X /Xj ) 2 and we expect a stability region little different from that of (p*,CT*) . For (p*,a*) = 4th order Adams-Bashf orth (5.2) is Schur in the following region (where we assume I I' It follows that for t = 2, . . . , m the interval of absolute stability is only slightly less than that of AB, . Let t = 1 in diagram, and for Q, ~ 2 read 1 - (X./X.) . Then stability is retained if (QphX.Q.) lies in the indicated parallelogram, which implies the much improved sensitivity requirement h|X.-X.| < .077 . §6 Predictor-corrector SPM A predictor (p*,a*) and corrector (p,a) in PECE mode usually has better absolute stability than (p*,0 ) . Applying the correction process (5.1 (iii) (iv)) after the PECE cycle - indicated by PECE£ - raises the sensitivity problems encountered 7-2 with (SPM), . However, for the mode VE^CEE, , (5.2) is (5.2) is replaced by tt = (l-Q t )(l+hX t 6 k Q t )r k + Q t [p-hX t a+hX t B k Q t (p*-hX t a*)] . An order of magnitude argument suggests a good com- promise between sensitivity and subdominant stabi- lity is got by setting to zero terms in hX Q , which forces the choice a = B.r , ie (p,a) is a backward differentiation formula (BDF) . Further, it is found to be advantageous to choose p = p . The PL pair is then completely specified by (p,a) = BDF of order k ; p* = p ; (p*,cr*) of order k . §7 Extensions The extension to the general separably stiff linear problem (ie more than one dominant eigen- value) is straightforward; but computing estimates for the dominant eigenvalues becomes more difficult. Nonlinear systems may also be tackled by evaluating X. = X.(x,y) , dominant eigenvalue of the Jacobian, at appropriate approximate values for y . Success will depend on the eigensystem of the Jacobian varying slowly with x . References [1] Lambert, J D "The numerical integration of a special class of stiff differential systems", Proc 5th Manitoba Conf on Num Math, Oct 1975. [2] Alfeld, P and Lambert, J D "Correction in the dominant space: a numerical technique for a certain class of stiff initial value problems", Math Comp Vol 31, No 140, 922-938, Oct 1977. 7-3 LB1.-8882 THE NEED FOR IMPROVED NUMERICAL ALGORITHMS IN MOLECULAR SCATTERING* Lowell D. Thomas National Resource for Computation in Chemistry Lawrence Berkeley Laboratory University of California Berkeley, California 9A720 USA Many problems in theoretical chemistry re- quire an understanding at the microscopic level of molecular scattering. That is, an understand- ing of what internal and external changes take place as two molecules approach one another, col- lide and then recede from one another. From this knowledge at the microscopic level one can then calculate many properties of macroscopic chemical systems which cannot be measured presently in the laboratory nor observed directly in our environ- ment. The applications of molecular scattering information are many and far-reaching. In laser physics this information is needed in the search for new chemical systems with better lasing pro- perties and for performance improvement of exist- ing ones. The discovery within the past ten years of a rich variety of molecules and mole- cular clouds in interstellar space has created a need in astrophysics for molecular scattering data. Many interstellar molecules cannot be made in earth-bound laboratories, leaving computation as the only method for the determination of their properties. The study of our own atmos- phere involves many molecular scattering pro- cesses, including, for example, the formation and destruction of ions in the upper atmosphere and the creation of air pollutants. Many energy- related problems also require understanding of molecular scattering events. To understand combustion processes requires knowledge of the associated reactive scattering events and in nuclear reactors there are many inelastic and reactive collisions involved in the final dis- position of the generated energy. This list- ing is just a brief mention of a few of the applications. Quantum mechanical approaches to mole- cular scattering can be reduced to the solu- tion of a system of second-order, ordinary differential equations of the form, *This research was supported in part by the National Resource for Computation in Chemistry under a grant from the National Science Founda- tion and the Basic Energy Sciences Division of the United States Department of Energy under Contract No. W-7405-ENG-48. y"(r) + k y(r) = V(r)y(r) k 2 >0 y(0) = lim y(r) = sin(kr+n) + cos(kr+n)R n constant (1) (2) (3) (4) (5) Here V(r) is called the potential energy function or most often simply the potential. Another impor- tant class of equations have non-local potentials, i.e. , y"(r) + k 2 y(r) =/ V(r ,r ' )y ( ')dr • (6) However, the discussion here will be limited to the local case. Eq. (1) will generally be a matrix equation with V(r) a real, symmetric matrix and k diagonal. In this case the sin and cos functions in Eq. (4) are to be regarded as diag- onal matrices and R is real symmetric. The matrix R is the desired quantity. From it all micro- and macroscopic scattering information can be derived. This problem is most commonly solved by start- ing with Eq. (3), assuming an arbitrary, non- singular matrix for y'(0) and then numerically integrating outward to some large r such that V(r)=0. the numerical solution then has the form y(r) = sin(kr+n)A + cos(kr+n)B (7) and the desired solution is found by multiplying from the right by A -1 . That is R = BA (8) The single equation problem illustrated in the figure below is representative of the nature of the problem. For small r where V(r)>k 2 , the solution has exponential behavior. For large r where V(r)< 1/2 and is second-order accurate when 6 = 1/2. For this example, algorithm (9) becomes (I 6AtJ )(u u ) Atf(u n ) This noniterative ODE method was suggested indepen- dently by Lomax [9] and Liniger and Willoughby [8]. This scheme is equivalent to solving (4) by Newton's method and stopping the iteration after one step. (c) Approximate factorization . For PDE applica- tions it is convenient to split (3) as (10) dt ' f(u) fj(u) + f 2 (u) The Jacobian is then 3u 3f, 3f, (11) J and the algorithm (9) becomes (12) [I-u)At(J n+p + J2 +P )]p(E)u n = Atf(u) . An approximate factorization of (12) yields (13) (I-uAtJ^ +P )(I-uAtJ2 +P )p(E)u n = Atf(G) . It is shown in [13] that this approximate factoriza- tion does not upset either the second-order temporal accuracy, or the unconditional stability of the method when applied to linear test equation du/dt = Xiu + A 2 u (ReX!,A 2 < 0). The A-stability is preserved because we have chosen p(E)u n as the unknown variable. For simplicity, we illustrate the PDE applica- tion by considering only the hyperbolic equa- tion (1). The algorithm for the complete equa- tion (2) will be given in the final manuscript. Comparing (10) and (1), we identify the differen- tial operators (14) fj(u) 3F(u) 3x f,(u) 3G(u) 3y The Jacobians (11) are replaced by the operators 3 (15) '1 3x 3y B where A = 3F/3u, B = 3G/3u are Jacobian matrices. (For the equations of gas dynamics in two spatial dimensions, these are 4><4 matrices.) Inserting (14) and (15) in (13) yields* (16) (i + ooAt £ A n4 )(l +uAt ^B n+P )p(E)u n = » t pF(u) + 5G(u) ] L 3x 3y J An obvious computational sequence to implement (16) as an ADI method is *In (16), notation of the form 3 denotes (I + uAt ■£- B) p(E)u dy p(E)u" + uAt ^ (Bp(E)u 11 ) 9-2 (17a) (i + uAt ^A n+P Jp(E)u* I a* 3y J (17b) (l + -.At ^ B n+P ) D (E)u n = o(E)u* , (17c) u where o(E)u* is a dummy temporal difference. The spatial derivatives appearing in (17) are to be approximated by appropriate finite-difference quo- tients. For three-point spatial difference approx- imations, the x- and y-operators on the left side of (17a, b) each require the solution of a block- tridiagonal system of equations with each block having dimensions m * m, where m is the number of components of the vector u. There is an obvious extension of (17) to three-spatial dimensions. In practice, two-step A-stable methods are of primary interest. The class of all two-step methods which are at least second-order accurate can be written as (18) (1 + Ou n+2 (1 + 2?)u n+1 + Cu" = At [ef n+2 + (e - 29 +|) f n+1 - ( 5 - 6 + f ) f n ] , where 9,£ are arbitrary real numbers. By applying the theory of positive real functions [4], one finds that (18) is A-stable if and only if 5 < 29 - 1 , S>--f MIXED SPATIAL DERIVATIV ES. The appearance of mixed spatial derivatives 3 z /3x3y on the right-hand side of (2) precludes an efficient reduction of a multidimensional problem into a product of one- dimensional operators. This difficulty can be avoided by treating the mixed derivatives with an explicit LMM which is "compatible" with an A-stable LMM. We demonstrate the rather remarkable result that this can be done and still retain unconditional stability. The only penalty is that the parameter space for A-stability is slightly reduced. NUMERICAL COMPUTATIONS . Finally, we give a brief survey [1,11,12] of recent two- and three-dimensional flow computations made using ADI-LMMs of the form (17). REFERENCES 1. R. M. Beam and R. F. Warming, An Lmplicit fac- tored scheme for the compressible Navier-Stokes equations, AIAA J. 16 (1978), 393-402. 2. G. Dahlquist, A special stability problem for linear multistep methods, BIT 3 (1963), 27-43. 3. G. Dahlquist, Error analysis for a class of methods for stiff non-linear initial value problems, Lecture Notes in Mathematics 506 , Springer-Verlag, Berlin (1976), 60-74. 4. G. Dahlquist, Positive functions and some applications to stability questions for numer- ical methods, preprint, MRC Symposium, Madison, Wise, April 1978. 5. J. Douglas, On the numerical integration of 8. u xx + u yy Indust . Appl i t by implicit methods, J. Soc. Math. 3 (1955), 42-65. J. Douglas and J. E. Gunn, A general formula- tion of alternating direction methods, Numer. Math. 6 (1964), 428-453. A. C. Hindmarsh, GEARB: Solution of ordinary differential equations having banded Jacobian, Lawrence Livermore Laboratory Report UCID-30059, Rev. 1, 1975. W. Liniger and R. A. Willoughby, Efficient integration methods for stiff systems of ordi- nary differential equations, SIAM J. Numer. Anal. 1_ (1970), 47-66. 9. H. Lomax, Stable implicit and explicit numeri- cal methods for integrating quasi-linear dif- ferential equations with parasitic-stiff and parasitic-saddle eigenvalues, NASA TN D-4703, 1968. 10. D. W. Peaceman and H. H. Rachford, The numeri- cal solution of parabolic and elliptic differ- ential equations, J. Soc. Indust. Appl. Math. 3 (1955), 28-41. 11. T. H. Pulliam and J. L. Steger, On implicit finite-difference simulations of three dimen- sional flow, AIAA Paper 78-10, 1978. 12. J. L. Steger, Implicit finite difference simu- lation of flow about arbitrary two-dimensional geometries, AIAA J. 16 (1978), 679-686. 13. R. F. Warming and R. M. Beam, An extension of A-stability to alternating direction implicit methods, submitted to BIT. Also NASA TM 78537, 1978. 9-3 TOLERANCE PROPORTIONALITY IN ODE -CODES Hans J. Stetter Tech. Univ. Vienna 1 . Meaning of Proportionality Most routines for the numerical solution of (1 .1 ) v' (t) = f(t,y(t)), y(0) = y , t £ [D,T], y(t) £]R S , request the specification of an accuracy level through a tolerance parameter 6 € IR. The choice of 6 determines the computation for a given problem (1 .1 ) so that all computational quantities are de- pendent upon 6 . Originally, the numerical solution of [1 .1 ] by a given code is only defined on a "grid" G(6) = {t ; n = 0(1 )N} c: [0,T] ; by a suitable interpolation n (see section 2) we may generate an approximation n(t,6} defined for all t€ [0,T], with global error e(t,6) :=n(t,6) -y(t). Tolerance proportionality means that (1.2) e(t,6) «*<5 • vCt), with a function v: [0,T]-»1R independent of 6, holds for an interval of 6-values including those commonly used. (1 .2) permits an assessment of e from the results of two computations with differ- ent tolerances 6 . A vaguer form of proportionality occurs in the Toronto Assessment Programs ([1], [2]). Here the efficiency of a code in finding an approxi- mate solution of a given accuracy for a given problem (1 .1 ) at a given point t is measured on the basis of 2 . Local equivalent of proportionality The numerical solution of (1 .1 ) by a one-pass forward step method cannot be controlled with a global strategy. But (1 .2), with v unspecified, is not a global requirement, at least not for small 6: Within first order perturbation theory it is equi- valent to ( cf . , e .g . , [3] ) the local requirement that n(t,<5) satisfies (2.1) n'(t,6)*«f(t,n(t,6)) +6 • y(t), with y : [0.T] -»1R independent of 6. (2.1 ) can be controlled on a step-by-step basis. For specified values n . and ri 6 IR , there n-1 n s is a unique vector u £ IR such that the solution of z'(t) = f(t,z(t) ) +u , z(t , n n-1 'n-1' satisfies z(t ) = n . We assume that o ur continuous n n approximation n(t,6) has been determined in this fashion; u may be called its "backward error" in n (t „,t ), cf. [4]. (2.1) requests n-1 n (2.2) u (6) «6 • y(t) in [t ,,t ]. n n-1 n It is well-known (see [3] or [4]) that u and the "local error per unit step" L are re- lated by (2.3) U (1 +0(h )); n n thus proportionality requires L to be proportional to 6. n (1.3) lle(t.6)ll«6 «v(t), where e and v(t) are determined by a logarithmic least squares fit. In the following, "proportionality" will refer to the stricter relation (1 .2) unless spe- cified. We will regard various effects in cur- rent ODE-codes for their impact on the achieve- ment of proportionality. For a code "of order p", with L = 0(h ), this 1/p n n implies h = 0(6 p ) . Thus, even if we manage to control the computation properly w.r.t. leading powers of h, we must expect (1 +0(6 )) -factors in the proportionality relation (1.2) if we re- place the ra-sign by equality. Therefore we need not bother about (1 + 0( 6) ) -factors which arise from linearization w.r.t. 6 in (2.1) or from the evaluation of quantities along y(t) in place of n(t,6). We assume that the actual control of the computation is effected through a vector L € ]R n which is computed in each step: On the basis of 10-1 L and 5, an acceptance procedure determines wheth- er to accept the current step or not, and a step - . and order) control procedure calculates the stepsize (and order) to be used in the next Cor re- peated) step. 3. Local Error Estimator Effects In a one-step method of order p, we may com- pute L from any procedure Lit „ ,n„_ , -h_jf ) which n satisfies n-1' 'n-T V (3.1) L(t,y(t).h,f) = h P (p(t+h)M+0(h)) xaptur Lffacts We assume that the acceptance procedure deci- des on the basis (4.1) ^ $6 accept HL n ll > 6 reject where II. ..II may include a weighting, e.g. by the current solution values. (The type of weighting and the choice of ]R -norm are supposed to be independent of 6.) If we have a local error estimation of type (3.1), we have if we only achieve (3.2) L = 6 • *(t )(1 +0(h )) n n n through our control. For a backward error u of order p. (3.1)/(3.2) implies (3.3) u = h P ip(t ) (1 + 0(h )) n n n n = h P 6, ~ n d-1 n L decreases like h and acceptance occurs essen- * n tially for (6.3) h a IIAf II si 6 n Similarly, the contribution of this step to the global error consists of two parts: The effect of the "smooth" backward error (normally very small due to the smallness of h ) and an immediate incre- n ment in e(t) due to a wrong consideration of Af. This second error is of the form h BAf, its effect on the global error further on is thus proportional llAfll. Due to the irregular dependence of o n Consider the scalar case, with ip(t) » a(t-t) in the vicinity of t. Then (5.2) and (3.1) imply n-1 t -t n n-1 so that the stepsize will grow locally until t >t; but the details of that growth will depend on the relative location of t in the grid and thus on 6 in an irregular manner. Thus there will be large u in the vicinity of t ( cf . (3.3) and note that ip(t) * 0), but their size will generally not be proportional to 6 . This may introduce a considerable perturbation of strict proportionality. 7. Starting Effects The choice-of -starting-step dilemma is well- known: On the one hand, it is unwise - and contrary to the intentions of numerical software design - to delegate this choice to the user; on the other hand, the values of y and f(t ,y ) are simply not suf- o o ■'o ficient to let the code choose a reasonable first step, cf . e.g. [6 ] . From the point of view of achieving proportio- nality, it is extremely undesirable to have an in- itial interval in which the stepsize (and hence the local error) is considerably smaller that the speci- fied tolerances would require; cf. section 4. The following suggestion seems to be practical: Design a procedure which computes an h. from y , f(t ,y ), and 6 which would be the proper step- size if some heuristic assumptions about the problem were true. Execute the first step and test for ac- ceptance. Repeat the step, with a new stepsize com- puted from (5.2), not only if IlL.II > 6 but also if IlL.II < c 6, with a reasonable value c < 1. De- tails may be found in [6]. Such a procedure will prevent ill-effects from unduly small starting stepsizes. The same procedure should also be used for restarting after a singular condition, like a jump in f, has been recognized by the code: -4 Assume that h has to be reduced to 10 in or- der that a discontinuity in f may be passed but that the smoothness of f after the jump permits steps of order 10"'' under the specified tolerance. With or- 10-3 itepsize control strategies, the computation many steps to return to the appropriate stepsize. Therefore a restart, even if tt involves a repetition of the first step will be much more jmic besides enhancing proportionality. which is indicated by the bi in fig. . without the (1 + 0(h))-fa< jrs gem. r under (8.1) and (8.2) .own. It may be ai that the differences are 0(h + it this i leading power on the level of the newly genera- ted error per step. -, 8. Output Effects So far we have assumed that steps are chosen exclusively to comply with local error require- ments (except in "emergency situations"). In many codes, however, these requirements are overridden by the requirement to reach output points precisely. If there is a majority of steps determined in this fashion, proportionality must obviously become com- pletely obsolete. But even if we consider output at the end T of an integration interval only, the situation is not trivial. A common and natural strategy is the fol- lowing: Determine h by the local error vs. toler- ance requirements; use the stepsize (8.1) h - r T-t if T € (t ,t +h _ ], n n n n+1 (T-t )/2 if T€(t +h _,t +2h A , n n n+1 n n+1 Also in codes where the output values are ob- tained by interpolation, care has to be taKen if strict proportionality is not to be abandoned in the output process. Fig. 2 shows the source of the difficulty. Some related analysis may be found in [ ]. e (t . 5)' fig. 2 n+1 if T?t +2h „, n n+1 Thus normally the last two steps before T are mo- dified. A simple analysis shows that - from the point of view of proportionality - it is better not to looK ahead beyond t +h „ : n n+1 F -t if T€(t ,t +h J, n n n n+1 (8.2) h n+1 = fig. 1 n+1 if T £ t +h „ n n+1 newly generated error at T A rvM n+2 This is best explained by fig.1. Assume . (hA). Here the local error is too small so that h is increased again. The same effects occur in one- step methods if we consider the one (principal) root for a dominant and a (mildly stiff) subdominant component of the solution. We may say that in such a situation the com- putation is controlled by stability rather than by accuracy. Since the stability requirement has no- thing to do with the specified tolerance 6, a com- putation which is predominantly controlled by sta- 10-4 bility cannot be expected to yield a global error proportional to 6. We may even define a given problem (1.1) to be stiff for a given non-stiff code at a specified tolerance , if the computation is no longer control- led by local error requirements throughout. In a high-order Adams PECE-code stiffness in this sen- se should be very common as a quantitative analy- sis shows: Define (hX) by hA (9.1) IS (hA) I 2 .9 e for hA < ( n X; paras. and let us characterize the location of (hX) by the corresponding number of "locally correct di- gits" In a "variable procedure" code we can, of course, not assume that the same procedure is used in a certain section of the integration interval for different <5 . This affects our considerations at their roots: In section 3, we were able to permit the fic- titious local error estimate (3.1), with ipfcp, be- cause we could assume (p(t) to be essentially inde- pendent of 6 ; thus a function y was generated in (2.2) which was also essentially independent of 6. If there may be different functions (p for different <5 , the argument in section 3 breaks down. Furthermore, even if we have an asymptotically correct error estimator and achieve II u II w 6 for all n, this does not imply (2.2) in the variable pro- cedure case. Different vectors Y(t) of identical norm may lead to quite different solutions of (2.1). Thus there is not much to expect with regard to proportionality from a variable procedure code, no matter whether it uses local extrapolation or not. On the other hand, this does not mean that propor- tionality, at least of the vague type, cannot occur with such codes: If there is sufficient regularity in the relation between the various high order deri- vatives, the resulting function value Y(t) in (2.2) and (2.1) may not depend strongly on the particular procedure which is used in the neighborhood of t. Actually, it has been experimentally estab- lished that, e.g., the variable order, variable step multistep code STEP which also uses local ex- trapolation exhibits a considerable degree of vague proportionality for a wide variety of problems. rTTTi (hX) ? • (hA J - e 10. princ. - log — ^— (hA) 11 . Conclusions Then a situation is stiff for a K-step Adams PECE-code (p=K+1) if the tolerance requests fewer than K 4 locally correct digits 1 .5 2.5 3.75 5 B.5 10 9.5 10. Change of order (method) So far we have made the silent assumption that the tolerance parameter 6 influences the computation through the selection of the stepsize only. However, some of the most efficient codes are based on varying order or varying basic inte- gration procedures. This is well-known for Adams multistep codes; the recent construction of families of embedded Runge-Kutta procedures ([9]) suggests a similar approach for one-step methods. This approach is further supported by results which demonstrate that lower order methods are more efficient at cru- de tolerances ( [l] j . The preceding considerations should have made it clear that the output of a robust and efficient code for initial value problems (1.1) cannot be ex- pected to exhibit more than the vague proportiona- lity between the global error and the specified tolerance 6 which has been described in section 1 . An exponent e » 1 in (1.3) should appear only if the least squares fit is based on a wide range of to- lerances . If 6 is varied only over a narrow range, with many samples taken, say 6 = 10~ 31 - ■ ' ' , and if the values of the individual error components are re- garded in place of the norm of the error vector, we should not be surprised if few traces of propor- tionality are observable before smooth ing. This means that we cannot reliably judge the size of the global error from running the problem twice with different tolerances. But this was the main advantage that we hoped to obtain. If we can- not have it we may as well forget about proportio- nality, at. least form a utilitarian paint of view. Even if our codes were more perfectly propor- tional, the ordinary user would certainly not bother to run them twice at different tolerance levels if they would give him an approximate value of the global error directly . I do not believe that he would be very fastidious w.r.t. the accuracy of this 10-5 iate value; the order of magnitude might be t in many cases. He might even accept an -;ional blunder without much fuss. (Precisely rehavior is shown by the users of quadrature lea where the reliability of the error assessment -. tter. ) If we admit that an assessment of the global ir should be incorporated into ODE-codes, the tion is whether it can be done more efficiently :-ing the codes more proportional or by a di- rect approach. From my experiences. I would favor the second choice; see [10]. Nevertheless, the analysis of the proportiona- lity properties of ODE-codes yields important in- sights into the mechanism of these codes. References ['] W.H. Enright: Using a testing pacKage for the automatic assessment of numerical methods for OOEs, in: Performance Evaluation of Numerical Software, Proceedings, North Holland, to appear . .H. Enright, T.E. Hull: A collection of pro- grams for assessing initial value methods, Dept. Comp. Sci. Tech. Rep., Univ. of Toronto, 1979. [3] H.J. Stetter: Considerations concerning a theory for ODE-solvers, in: Numerical Treat- ment of Differential Equations, Springer Lec- ture Notes Vol.631 (1978) 186-200. [4] K.R. Jackson, W.H. Enright, T.E. Hull: A theo- retical criterion for comparing Runge-Kutta formulas, SINUM r5 (1976) 618-641. [5] H.J. Stetter: Interpolation and error estimation in Adams PC -codes, SINUn 16_ (1979). [6] L.F. Shampine, H.J. Stetter: Starting an ODE- solver, to appear. [7] B. Lindberg: Characterization of optimal step- size for methods for stiff differential equa- tions, SINUM 14_ (1977) 859-887. [B] H.J. Stetter: On the maximal order in PC-codes, Computing 20 (1978) 273 - 278. [9] D.G. Bettis: Efficient embedded Runge-Kutta- methods, in: Numerical Treatment of Differen- tial Equations, Springer Lecture Notes Vol.631 (1978) 9 - 16. .J. Stetter: Global error estimation in ODE- solvers, in: Numerical Analysis, Springer Lec- ture Notes Vol.630 (1978) 179-189. 10-6 Numerical solution of large systems of ODE's with sparse Jacobians by Per Grove Thorns en Department of Numerical Analysis Technical University of Denmark. Introduction . In many technical applications, model problems are defined by systems of ordinary differential e quations f(t, y ; t e [o,t] y e R (D where n is a large number, and where the Jacobi- an matrix -fil i,j = 1,2,. ,,n (2) is a sparse matrix. In this paper we shall discuss the problem of selecting a solution method and show how sparse matrix techniques can be used to solve the systems of linear equations that arise, and how it may be implemented efficiently. Solution methods . The systems ( 1 ) that concern engineers will usu- ally be characterised by having a wide spectrum of eigenvalues for the Jacobian J, and these will of- ten vary along the solution. This means that only methods for stiff systems can be used in general purpose software. On the other hand the engineer will often have very modest requirements of accu- racy. In such cases methods of constant order will be reasonably efficient, for the present discussion we have selected order two. Often it will be advantageous for the stepsize strategy to take large steps in spite of the stiff character of the system, therefore it is important that the region of absolute stability of the method contains points at infinity in the left half plane, and that damping of errors in the stiff components will take place. In some applications one would also be interested in damping of highly oscillatory components . These objectives can be achieved if, in the case of one-step methods we have L(a )-stability , while in the case of linear multistep methods we could use the following L, (a )-stability definition: (see ref.fi]). Definition: A linear k-step method k k I a.y . = h I B.f .. j=o J n+J j=o J n+J is L, (a)-stable if (i) it is A(a)-stable, and if k applied to the testequation y' = Xy , Re(X) < , it yields |y n+k l <*(h*) ma *|y 11+i l < i < k-1 , such that |R(hX)| -» when |x| -» » , |arg(-X)| < a . Of the methods that satisfies the above require- ments we have considered two types, the Backward Differentiation methods and Semi-Implicit Runge- Kutta methods, for order two we will use the fol- lowing formulae. n+1 U 1 2 3 n 3 n-1 3 n+1 (3) for error estimation, an approximation to y'" may be found by using backward differencing. The SIRK- methods are of N0rsett-type (see ref.[2]). k, = f(t +yh,y +Thk. ) (U) 1 n n 1 \ = f(t n+ 6 2 h,y n+621 k lh+ yhk 2 ) k 3 = f(t n+ 6 3 h,y n+ 6 31 k 1 h + 63 2 k 2 h +Y hk 3 ) V1 = y n + ^VWWV (5) LE = h(c 1 -k 1 +c 2 -k 2 +c 3 -k 3 ) (6) 11-1 on . h variable stepsize implementation the step- . are in general controlled by the estimated local errors. In the stepsize selection strategy the error constants will influence the stepsizes chosen in such a way that the smaller the error jtant is, the larger the stepsize will be chosen. In our comparison the SIRK will have a smaller r constant than the BDF method. If the same error tolerance is prescribed we may expect the SIRK to do better on average than the BDF be- cause the number of steps will be smaller. On the other hand, Ngrsett has shown (ref.[2]), that it is possible to select the parameters in C*,5»6) in such a way that the k found in step n can be used as k in step n+1 , so that only two k.'s have to be found in each step, except when starting and when changing stepsize. Each of the formulae (3) and (h) is a nonlinear system of equations that has to be solved in each step, this solution is in general solved by means of quasi-newton iterations, that take the following forms , in the BDF case: in the SIRK case: Vi * 2 -■— ► t (I-YhJ)e [ * ] =k U] -f(y +vhkj s] ) 11 n i k U+l] =k [s] -e[ s] i li (9) (10) The matrices A will be changed when either the stepsize is changed or when the convergence is slow and the Jacobian updated. If our system is large the factorizations of the matrices A will r be very expensive and the main part of the total work involved is spent on that. If sparse tech- niques is used the work can be reduced but it is essential that as much information is carried a- long about the earlier factorizations. Application of sparse techniques . The subroutine SSLEST has been constructed so as to facilitate the solutions in such cases. By using ordered lists for storing only nonzero ele- ments and by selecting pivots so that the number of fill-in's is small, each factorization is made very efficient. Further if as in the case above the se- quence of matrices A all have the same structure, r then the information in the pivotal strategy may be used in subsequent factorizations, thereby saving a great amount of work, in general this amount will be about 60% of the total factorization. This sparse matrix package has been used in connection with the SIRK method in a subroutine SPARKS that will be presented, and compared to the use of BDF methods . Conclusions. In each case J is an approximation to the Jaco- bian, this is only updated in case the convergence is too slow. It is convenient to consider the pro- cess written in the form A r 6 n-1 " b n+1 (11) where A = (I - YhJ ) and b , is one of the r n+1 right hand sides (7) or (9) both formulae have this general form. If we consider the time history of the solution we will find a sequence of matrices A , r ■ 0,1, that can be displayed in the form The use of sparse matrix techniques in connec- tion with methods for stiff systems for ODE's with sparse Jacobians, can be very effective. Two types of methods have been compared, the BDF methods versus the SIRK methods , both are suf- ficiently stable to allow the use of large step- sizes. Efficient implementations with a sparse package is presented. 11-2 References . [1]: Z. Zlatev, P.G. Thomsen : An analysis of the solution of linear systems of ordinary dif- ferential equations. Inst, for Num. Anal., Technical University of Denmark, June 1977, NI-77-10. [2]: S.P. N0rsett : Semi Explicit Runge Kutta methods, Mathematics and Computation No. 6/Th. Dept. of Mathematics, University of Trondheim, ISBN 82-7151-009-6. [3]: Z. Zlatev, V.A. Barker, P.G. Thomsen : SSLEST A Fortran subroutine for solving sparse sys- tems of linear equations, Users guide. Inst. for Num. Anal. , Technical University of Den- mark, NI-78-01. 11-3 AUTOMATIC PARTITIONING OF STIFF SYSTEMS AND EXPLOITING THE RESULTING STRUCTURE W.H. Enright and M.S. Kamel Department of Computer Science University of Toronto Abstract Most methods for solving stiff systems are based on implicit formulas and require the use of Newton-like iterations. The cost of the matrix operations in the iteration scheme of these methods can be quite high. A new iteration scheme is developed which exploits the structure of the system and also allows fast updating of the iteration matrix after a stepsize or order change. The technique is particularly useful when the stiffness is due to only a few components of a large system, and is applicable to most methods for stiff systems. It is based on an automatic partitioning of the system into subsystems corresponding to stiff and non-stiff parts. This partitioning is then ex- ploited in the solution of the system of equations that must be solved on each time step. 1. Introduction In most methods for solving stiff systems, the solution of a system of equations on each time step by a modified Newton iteration usually consumes a considerable amount of the compu- tational effort of the method. One way to achieve efficiency is to take advantage of the structure of the problems and develop techniques to solve problems of the same structure. For example, problems with banded or sparse Jacobians can be solved efficiently by methods with iteration schemes based on banded or sparse (linear system solvers), rather than using general methods (see Hindmarsh [1976]). If the system can be partitioned into stiff and non-stiff parts (transient and smooth com- ponents), more effective techniques may be developed to reduce the cost of the iteration scheme. The idea of partitioning the system into transient and smooth parts is an old idea and is related to the singular perturbation problem, where the system y(t) - f(y(t)), y(0) = y Q (1.1) is transformed into Du f 1 (u,v)) f 2 (u,v) i (1.2) This research was supported by the National Research Council of Canada where D is a matrix which may be singular u corresponds to transient components v corresponds to smooth components. As early as 1952, Dahlquist used a method called Pseudo Stationary Approximation (PSA) in which the transient components are eliminated by putting D=0, thus reducing the system to a dif- ferential algebraic system. The integration can then be done by any conventional ODE solver. This method was revised by Dahlquist [1969] and analysed and tested by Oden [l97l]. The new version is called SAPS (for smooth approximate particular solution) . Partitioning of linear systems of the form, y = Ay + b(t), where A is a matrix with constant elements and with the large eigenvalues well separated from the small eigenvalues, has been considered by Lee [1967] and by Branin [1967] in connection with circuit theory. Lee used matrix filtering to divide the system into loosely coupled subsystems, one for the smooth components and one for the transients. The smooth and transient components are obtained by using the transformation y 1 = LP y, y 2 = HP y; where LP and HP are matrices which act like filters. LP passes the eigenvectors of A that correspond to small eigenvalues and annihilates eigenvectors that correspond to large eigenvalues. HP has the opposite effect. The integration procedure he used is similar to the PSA method where the transient components are eliminated (by putting y„ = 0) and the smooth components are integrated by a conventional method. Branin suggests that the matrix A be de- composed into a sum of diagonal and off-diagonal parts. If the diagonal part has large eigen- values and the off-diagonal part has small eigen- values, the restriction on the stepsize for the formula he used (which is based on an analytical solution) is relaxed. He noted that, for most practical problems, the matrix A has to be trans- formed first so it can be decomposed in this way. Branin suggests using some cycles of the QR- algorithm to do the transformation. He also gives guidelines to handle nonlinear systems, where the Jacobian matrix plays the role of the matrix A, but is changing during the course of the integration. In this case the transformation has to be frequently updated. 12-1 Methods based on implicit formulas have been widely used for solving stiff systems. However, using these methods may not be efficient for large systems in which the stiffness is caused by only a few components of the system. The cost of dealing with the Jacoblan of the full system can be quite high. Dahlquist has pointed out that more work should be done on designing software which can utilize the favorable situation which occurs when some subsystems are not stiff (Dahlquist [1974]). Alfeld and Lambert [197 7] report on a method where the correction is only done in the 'domi- nant' space (the space spanned by the transient components). They consider systems of the form, y = A(t)y + g(t), where the eigenvalues of A(t) may be separated into two parts, one of which dominates the other. The integration consists of taking one forward step by a conventional multistep method and then making a correction in the subspace spanned by the eigenvectors corres- ponding to the dominant eigenvalues. The technique requires the explicit computation of the dominant eigensystem. Application to non- linear systems is also discussed where the eigensystem of the Jacobian will be needed at each step. Recent work on partitioning stiff system has also been reported by Hoffer [1976]. Hoffer's method is designed to solve large systems in which the transient components originate from only a few, easily identified equations. His integration formula is a combination of the modi- fied midpoint rule applied to the subsystem of the smooth components, and the implicit trape- zoidal rule applied to the transient components. In the partitioning techniques described above the difficulty of solving linear systems of equations on each time step is reduced somewhat by applying an analytic expansion or implicit formula (with the modified Newton iteration) only to the stiff subsystem and applying a standard non-stiff formula to the remainder of the system. The accuracy and stability of such techniques is very sensitive to the correct identification of the corresponding subsystem as well as to the amount of coupling that is present between these subsystems. It should also be noted that, even for linear problems, the appropriate partitioning will in general be time and tolerance dependent. The technique we propose is based on applying the same implicit formula to all components of the system and exploiting partitioning only in the iterative solution of the resulting systems of equations. As a result only the cost of the Iteration will be sensitive to the correct iden- tification of the corresponding subsystem and the accuracy and stability of the integration will not be sensitive to the partitioning. In a related study, Robertson [1976] has con- sidered partitioning of the Iteration matrix in the iteration scheme. By fixing the part of the iteration matrix which corresponds to the smooth components and only updating the part which cor- responds to the transient components, the iter- ation scheme will result in using a modlfied- Newton technique only on transients and repeated substitution with correction terms for smooth components. By allowing periodic updating of the part corresponding to the smooth components further efficiency can be achieved. Robertson does not discuss how one can obtain the appropri- ate partitioning. Techniques that achieve efficiency in solving stiff systems without pre-information about the transient and the smooth components are more general and can easily be incorporated into existing packages for stiff systems. The fast update technique developed by Enright [1978] Is of this kind. The iteration matrix is reduced to Hessenberg form and subsequent updates of the representation of the iteration matrix, after stepsize or order changes, require fewer oper- ations than working with the inverse of the iteration matrix, or its LU factors. The technique to be developed and analyzed in this study is related to Enright 's approach. As the iteration matrix is reduced using a sequence of similarity transformations, a test is performed to identify the stiff components. The resulting partitioning of the system leads to an iteration scheme which adapts to the behaviour of the problem and can be interpreted as applying a modified Newton iteration to the transient components and repeated substitution to the smooth components. (The idea is actually similar in spirit to the identification and solution of rank deficient least square problems by the method described in Lawson and Hanson [l97A,pp. 77]) . Substantial saving of computa- tional operations can be expected, especially if the stiffness is due to only a few components of a large system. 2. The Proposed Technique The application of an implicit multistep formula to the system y = f(y), y(t Q ) - y Q , can be written as = G(y 1+1 ) = y i+1 -h6 f(y, +1 ) (2.1) i+r ■E Yi+i-j " h £ 3-i 3=1 3 y i+l-3* (2.2) An iteration scheme for solving the nonlinear system G (y i+1 ) ■ 0, can be viewed as: Predict y , set 6 i+1 'i+1 y i+l for l = 1 ' 2 and determine <5. by solving the linear system i+1 W6 fc-1. i+1 G(y~ + p for I - 1,2,.... (2.3) Different iteration schemes will result from different choices of W. For example, if W is the identity matrix then the iteration corresponds to repeated substitution or, in our context, a predictor corrector iteration; while if W = [l-h6 A] where A ■ 3f/3y evaluated at y. (3 * 1+1). the iteration corresponds to a modified Newton scheme. If (2.1) is stiff then there is an obvious trade-off between these two possible choices of W. The former choice is cheap (0(n) operation 12-2 per iteration) but in order to guarantee con- vergence the stepsize must be severely restricted (| |hB Q A| |<1). On the other hand, the latter choice is expensive (0(n~) operation per iteration) and 0(n ) operations when W has to be updated), but much larger stepsizes are possible. Basically the scheme we describe is a compromise between these two choices of W. Our idea is to approxi- mate the matrix A, of the Newton iteration by a matrix A that will permit efficient iterations and updates, and at the same time, include enough information for the iteration to converge. Before we discuss the technique we should review the technique suggested by Enright [1978] for implementing the modified Newton iteration. He observed that if the system is of dimension n and W ■ [i-hB-AJ is retained in its LU factored form then 0(n ) operations will be required wher- ever the scalar hB n changes. On the other hand -1 if W is retained in its LHL factorization (H is an upper Hessenberg matrix), then fast updates are possible when h6_ changes and only when A needs to be updated will W require a full update. In order to motivate how we choose X, let us rewrite the modified Newton iteration as: accomplish this by making ||R, 9 || small (A - (l/hB )I)<5* +1 = (l/hB )G(y^ + J) (2.4) Consider the factorization of A by the similarity transformation: A = S R S where R is of the form (2.5) H R 12 \° R 22/ and, H is an k*k upper Hessenberg, R.. - is an kx(n-k) matrix, R-- is an (n-k)*(n-k+l) matrix. The matrix R is of the same structure as the reduced matrix, A, , at the k-th step in the re- duction of A to upper Hessenberg form. That is, if the reduction were done using Householder transformations then A k" U k-l U k-2--- U l AU l"- U k-2 U k-l where U , I = l,...,k-l are elementary reflectors. We modify this approach by introducing pivoting (row and column interchange on each step). R is then the reduced matrix at the k-th step: R = Vi p k-i--- u i p i AP i u r-- p k-i\-i (2 - 6) if R by che matrix | is small enough, then we can replace M2 (2.7) and A can then be replaced by X = SRS -1 . Note that ||A-A|| F < ||R 22 M F . Returning to the system (2.4) and substituting A for A we get (A-(l/hB )I)6* +1 = S(R-(l/h8 )l)s" 1 6j +1 = (l/hB )G(yJ;J) (2.8) A decomposition algorithm has been developed which applies the transformation as in (2.6). The interchange of rows and columns is designed so that the columns of R „ with large norms will be the pivot columns. The value of k is determined as the smallest integer such that ||R -||_ is less 22 r than a specified threshold value, T. The cost of this decomposition is approximately 10/3 n k 2 3 -7/3 nk +4/3 k adds and multiplies. The matrix H is updated by adding (-1/hB-Jl and the resulting matrix is factored into its LU 2 decomposition at a cost of k /2 adds and multiplies. Then, to solve (2.8) we must apply the following four steps: "c 1. form the vector c 2. set z 2 = -(hB Q )c 2 L°2j B^tt/h )G(yJ ) 3. solve Ly 4 . solve U C l ' R 12 Z 2 5. set & i+1 The cost of solving one system is then approxi- mately 5 n k - 5/2 k adds and multiplies. In the remainder of this paper we will discuss the selection of the threshold value, T and the application of this technique to existing methods. We will also present some numerical results illustrating the performance of this technique. where P is a permutation matrix. Let P U = S . Then the matrix S in (2.5) will be S - S . . .S 1 k-1 -1 T and S - S . One would prefer, in our application, that the k*k matrix, H, correspond to the transient components. The interchanges are designed to 12-3 References P. Alfeld and J.D. Lambert [1977]. Correction in the dominant space: A numerical technique for a certain class of stiff initial value problems, Math. Comp., vol. 11, No. 140, pp. 922-938. F.H. Branin [1967]. Computer methods of network analysis, Proc. IEEE 55, No. 11, pp. 1787- 1801. G. Dahlquist [1969]. A numerical method for some ordinary differential equations with large Lipschitz constants, in Information Processing 68, ed. A. Morell, North-Holland Pub. Co., Amsterdam, pp. 183-186. G. Dahlquist [1974]. Problems related to the numerical treatment of stiff differential equations, Proc. International Computing Symposium 1973, eds. A. Gunther et al., North-Holland Pub. Co., Amsterdam, pp. 307-314. W.H. Enright [1978]. Improving the efficiency of matrix operations in the numerical solution of stiff ODE's, ACM TOMS, 4, No. 2 pp. 127-136. A.C. Hindmarsh [1976]. The LLL family of ordinary differential equations solvers. Lawrence Livermore Laboratory, Report UCRL-78129. E. Hoffer [1976]. Partially implicit method for large stiff system of ODE's with only few equations introducing small time constants, SIAM J. Num. Anal. 13, No. 5, pp. 645-664. C.L. Lawson and R.J. Hanson [1974]. Solving least squares problems, Prentice Hall, Englewood Cliffs, NJ H.B. Lee [1967]. Matrix filtering as an aid to numerical integration, Proc. IEEE, vol. 55, No. 11, pp. 1826-1831. L. Oden [l97l]. An experimental and theoretical analysis of the SAPS method for stiff ordinary differential equations, Dept. of Information Processing, Rep. NA 71.28, Royal Inst, of Tech., Stockholm. H.H. Robertson [1976]. Numerical integration of systems of stiff ordinary differential equations with special structure, J. Inst. Math. Appl. 18, No. 2, pp. 249-253. 12-/1 IMPLEMENTATION OF IMPLICIT FORMULAS FOR THE SOLUTION OF ODES L. F. Shampine Applied Mathematics Research Department Sandia Laboratories Albuquerque, New Mexico The material presented in the talk appears in written form as the report SAND79-069A of the same title. Copies may be obtained from the author on request. Implicit formulas are quite popular for the solution of ODEs. They seem to be necessary for the solution of stiff problems. Every code based on an implicit formula must deal with certain tasks studied in this paper: (i) A choice of basic variable has to be made. The literature is extremely confusing as to what the possibilities are and the consequences of the choice. These matters are clarified. (ii) A test for convergence of the method for solving the implicit equations must be made. Ways of improving the reliability of this test are studied. (iii) Deciding when to form a new approximate Jacobian and/or iteration matrix is a crucial issue for the efficiency of a code. New insight which suggests rather specific actions will be developed. The report closes with an interesting numerical example. It was created to illustrate in a simple way a difficulty observed in practice. It also shows that in very reasonable circumstances, (i) the implicit formula can have more than one solution, (ii) the predictor can be closer to the "wrong" solution, and (iii) any convergence test based solely on the difference of successive iterates can accept an approximate solution from a divergent iteration. 13-1 TAYLOR SERIES METHODS WITH VARIABLE STEPSIZE RATIO AND VARIABLE ORDER George Corliss Dept. of Math, and Stat, Marquette University Milwaukee, WI 53233 Y. F. Chang Dept. of Computer Science University of Nebraska Lincoln, NE 68588 Partial support for this work was provided by the National Science Foundation under Grant No. MCS78-03022. Taylor series methods compute a solution to an initial value problem in ordinary differential equations by expanding the solution in a long series. This solution is extended across the desired interval by analytic continuation. We present techniques for choosing the largest ratio of the integration stepsize to the series radius of convergence (Re) which may be taken subject to error control constraints. This optimal ratio depends on the length of the series, the orders and locations of the primary singularities, the desired error tolerance (e) and the type of error measure the user wishes to control. Then a series length can be chosen which minimizes the computational cost of solving the problem. Most state-of-the-art codes for solving ODE's vary both the stepsize and the order of the formula used in order to handle efficiently a variety of problems and a range or error toler- ances. Several authors have written translator programs which accept the statement of the system of ODE's to be solved and produce object code which can run later to solve the system using series methods. The lengths of the series used varied from 11 [1] up to 32 [5], but none of these codes attempted to determine the most efficient series length at run time. We have successfully solved examples using series as long as 1000 terms. Chang [2] gave an example for which 20 terms was found empirically to be a good length. Each author varied the stepsize to keep it substantially less than some estimate for Re to con- trol the error. By using a relatively long (20 to 30 terms) series, Chang was able to give an accur- ate estimate of Re. Then a step of approximately 1/CH-l) Re- e was used. The choices of an optimal stepsize and series length require the accurate estimation of Re from tests Involving several consecutive terms of the series outlined in [3]. Singularities in the solutions to real ODE's can occur with any order, but they must occur either on the real axis or in conjugate pairs. Without loss of generality, the solution y(x) may be expanded at x - with series OO CO Z y (J) (0) trVj! - E Y . We assume for the j-0 j-0 J+1 moment that there is only one singularity on the circle of convergence. Solutions with more than one primary singularity are considered later in this paper. If the distance to the closest secon- dary singularity is R2 , then the contribution from the secondary singularities to the series is 0(|Rc/R2| N ) as N +«. Hence the effects are neg- ligible if a series is long enough. Consequently we consider the model problem v(z;a,s) = (a - z) -s , for s t 0, -1, ... . N Let v(z;a,s) Z V, be expanded using a j=l j step which satisfies < r = h/a < 1. Let N be i u ~j c |N-l+s.hi . large enough to satisfy — <1, N a v(z) N Z V j* d = N-l+s N ' r V (1 r), and B rdV N /(l rd). Then [A] showed that < A i T N 1 B > for 1 £ s; and B < T s. A and B depend on the h which was used to expand the series, but the series can easily be shifted to a new h* by multiplying the j-th term by (h*/h) J ~ 1 . We wish to choose an optimal stepsize ratio U = h*/a as large as possible while maintaining one of the error measures: Type i) absolute truncation error per step, |T N N or Type ii) absolute truncation error per unit step, |T N /h|; less than a prescribed tolerance, e. Relative or mixed error measures are handled in a similar manner. For Type i) error if s < 1, then U is the single real root on [0, 1] of V. P(x) N N-l Let 6 « /v N . N j. x + ex Then^U 1.1/N N ) is an approximate root of p(x), so we use it as an estimate for U. Estimates for other error types can be computed similarly, and some are given below. Figure 1 shows an example of how these estimates depend on the tolerance, c, and on the series length, N. 14-1 1 Error Estimated Optimal Measure Stepsize Ratios. Type i) s < 1 6 1/N (1-|6 1/N ) s > 1 (f) 1/N U-f 6 1/N ) Type ii) s < 1 1 (a6) N-1 (l - NTi(a6) 1 ■-1, s > 1 1 .a6.N-l.. d . x . <*-> (1 " FT U6) 1 Table 1. Optimal Stepsize Ratios We have shown how to compute an optimal stepsize ratio. Next we consider the choice of the series length. The codes produced by the various translator programs use recurrence rela- tions to generate the series. Fast series generation techniques have not yet been found suited to this application, so each integration step requires aN2 + 3N + y floating point opera- tions. The approximate cost of solving the entire problem is proportional to qN 2 + gN + y U As an example, the code produced by Chang's ATS 2 translator for y" = y - (y') requires approxi- 2 mately 1.5N + 15N + 100 floating point operations per step. The translator programs can be modified to supply their object codes with the values for a, 6, and y so that C(N) can be evaluated at run time to choose an optimal C(N) series length, N opt As experience with other methods has indicated, problems with a more stringent error tolerance are solved more efficiently by higher order methods (longer series). As indicated by Figure 2, C(N) typically has a broad minimum, so the exact choice of N is not critical. In practice, a series length longer than N should be used because the increased accuracy of the estimates for Re from a longer series balances the in- creased cost of its generation. The preceeding discussion assumed that the solution to the ODE had only one primary sing- uoarity. If the solution has a conjugate pair (a and ~a) of primary singularities of order s, then its series is asymptotic to EF . , the series for f (x) = (a - z) 8 + (a - z) S . j' If EV is the series for v(z) ■ (|a| - z) 2V cos(s+J)9, where 6 = truncation error for f it then F, j+l arg a. Hence the Ls at most twice the truncation error for v and the estimates for U given in Table 1 for a single primary singularity remain valid for a conjugate pair of primary singularity remain valid for a conjugate pair of primary singularities, provided that £ is re- placed by e/2. 50 100 SERIES LENGTH 150 200 FIGURE 1. OPTIMAL STEP FOR ABSOLUTE TRUNCATION ERROR PER UNIT STEP. d cm 10" 12 o CO -— • on i <=> LU <\j 10"6- X o 1 — • CO °° o <-> o ■ n 2k+2 := n 2k+l + hf (t 2k + 2' n 2k + 2 } ' This scheme may be regarded as an implicit trapezoidal rule on a 2h-grid: n 2k + 2 " n 2k = h [ f(t 2k' n 2k } + f(t 2k + 2' n 2k + 2 } ] started by two implicit Euler steps on an h-grid. The intermediate quantities r\ are only used to con- struct useful starting points for a Newton-like iteration to solve the implicit step (details omitted here). Moreover, in the implicit steps, an additional stepsize restriction is derived on the basis of the associ- 2 ated affine invariant convergence theorem. Once more, an h -expansion exists. The three above sketched extrapolation methods have been tested extensively, among other examples also by means of the well-known sample set STIFF DETEST. The results seem to indicate that the semi-inplicit mid- point rule in combination with the implicit Euler starting step is the preferable choice. All of the methods compare quite well with the code GEAR (due to Hindmarsh) . Detailed comparisons will be given. 15-2 WHAT MIGHT ODE SOLVERS BE TESTED FOR? T.E. Hull Department of Computer Science University of Toronto Toronto, Ontario, Canada Outline I have found it convenient to divide the material of this talk into three parts. The first part is intended to provide a brief summary of the testing that took place during the "classical" period of the 1950's and 1960's. The second is intended to provide a similarly brief summary of the "modern" period, the 1970's. The third and more extensive part consists of some speculation about what we might look forward to in a "future" period of testing of ODE solvers. The emphasis is primarily on testing in the experimental sense, that is, on testing programs for solving ODEs on a computer. However, testing in the theoretical sense is included as well, as for example in testing mathematically that a cer- tain class of formulas is asymptotically unstable, since it seems to me that the two aspects of test- ing complement one another and are best considered together. In a relatively brief survey of such a large field, it is impossible to be complete, and in particular to attempt anything like a complete bibliography. However, several recent papers and books do contain substantial lists of references. I therefore decided to omit references altogether from this summary, and only mention in the text itself a number of authors' names and dates. Those that are mentioned are intended to be reasonably representative, but by no means complete. The next section of this summary is devoted to the first, or classical stage of the 1950-60's. This stage is characterized primarily by a study and comparison of mathematical formulas, as opposed to anything like what we would now call mathematical software. The mathematical results obtained during this period were quite substantial, however, and had very important practical impli- cations, particularly with regard to the stability properties of formulas (first for non-stiff, and later for stiff problems). Experimental testing provided some important results as well, for example in pointing to the possible importance of variable order Adams formulas, but the experiments were relatively modest comparisons of various formu- las using fixed step sizes. Then a further section of this summary is ' <-d to the second, or modern stage of the 1970'b. This stage is characterized by a shift in emphasis to relatively large scale testing of "real" computer programs over fairly substantial classes of problems. There were also a number of important theoretical contributions related to the testing of various approaches during this period, for example with regard to the stability of methods for stiff systems, but the emphasis in testing had certainly shifted towards large scale experiments. Based on the results obtained, a number of conclu- sions have been reasonably clearly established about the reliability and efficiency of various classes of computer programs, particularly for non- stiff problems, but also to a limited extent for stiff problems as well. However, the developments and conclusions attributable to this stage brought with them many new and interesting questions about what we should be trying to accomplish when we test ODE solvers. The final section of this summary is intended to help clarify a few of these questions, and perhaps provide some suggestions about what steps we might take to resolve some of them. It is clear that our objectives must be more clearly defined, that we must broaden the classes of problems over which the tests are performed, and above all that our testing must reflect as closely as possible the real needs of the ultimate user, for example in terms of such factors as convenience, robustness, flexibility, and so on. Such changes will of course reflect back in turn on design objectives for the development of good numerical software for solving ODEs. Classical stage As already indicated, the important develop- ments of this stage appear to be more theoretical than experimental. The beautiful results establish- ed by Dahlquist in the mid 1950's can be considered, from the point of view of testing, as providing a definitive test for determining which multistep formulas are asymptotically stable and which are not . Then in 1963 Dahlquist "did it again", when he published a paper on A-stability which provided the impetus for much of the subsequent growth of interest in methods for solving stiff problems. The ideas and the criteria were essentially theore- tical. Other theoretical developments were also taking place in the 1960's. With the help of hindsight, we see that among the more notable from the point of view of testing mathematical formulas 16-1 for solving ODEs, were the extrapolation method 1 by Cragg (and extended by Bulirsch and Stoer) and the Runge-Kutta formula pairs developed by Fehlberg and others. The experimental testing of this period was relatively limited, and was restricted mostly to comparing formulas of fixed order, using fixed stepsizes, on relatively small numbers of problems. For example, some attempts were made, independently by R.ilston and Johnston, to determine which were the better Runge-Kutta formulas. Of the testing I associated with during this period, I believe that the results obtained with Creemer and publish- ed in 1963 were the most significant. After several years of trying to determine the best combination of accuracy and stability in multi-step formulas, we finally obtained experimental evidence that con- vinced us we should use Adams formulas, and worry about order, and not about more accurate formulas of a particular order. We also concluded that one corrector per step was optimal, but that this corrector needed to be followed by an evaluation (PECE) for stability. We didn't foresee the development of variable order, variable step Adams codes which are so common nowadays, but I wish we had! Modern stage During the 1970's there has been a rapid development of techniques for handling stiff systems, including theoretical criteria for testing the stability of various formulas. Theoretical criteria were also developed for non-stiff methods. For example, criteria for variable order, variable step methods were developed at the University of Illinois. More recently a theoretical criterion for testing Runge-Kutta formula pairs was developed at the University of Toronto. Led by the publication of Gear's DIFSUB, sub- stantial codes began to be widely available during this period. Some of the better known ones in North America are associated, for example with the names of Bulirsch and Stoer, Hindmarsh and Byrne, and Shampine and Watts. But others have been developed both on this continent and in Europe. (For the purposes of this talk, Britain is consi- dered to be a part of Europe.) Experimental testing of such programs also began to be carried out on a relatively large scale. An early series of tests was carried out at the Bell Telephone Laboratories. Even more extensive testing has been done at the Sandia Corporation, at the University of Manchester, and at the University of Toronto. The emphasis in all of this testing has been primarily concerned with reliability and efficiency, in one form or another, of the various methods. The Toronto testing program is called DETEST. It has changed considerably over the years, but basically it collects statistics on reliability and efficiency over several classes of problems. Care was taken, in all the ways we thought were possible, to produce results that were not affected appre- ciably by relatively separate issues, for example by issues such as determining the starting stepsize. The original design was aimed specifically at comparing methods, as opposed to simply making In- dependent assessments of them. This caused us to Insist that each method try to accomplish the same task, namely to keep the error per unit step below the prescribed tolerance, and this led us to modify the programs accordingly. This caused a great deal of controversy. Moreover, an "error per step" rule is more efficient and, in addition, makes it possible to step over the most common kinds of discontinuities. The present version of DETEST therefore does not require any modification of the program being tested. It is aimed more at the assessment of methods, rather than their comparison. However, a new notion of "normalized statistics" has been introduced to provide an indirect method of comparison. The present version is also much more flexible. For example, it is easy to select what problems are to be used, or to add new ones, and global statistics can be obtained, as well as local ones. NEW DETEST also contains a class of problems with discontinuities. Out of all the testing that has been done on non-stiff problems, a number of general conclusions have become established. For example, one happy conclusion is that there do not seem to be any serious "class distinctions", by which I mean that methods that are relatively good for one class of problems seem to be relatively good for the other classes as well. A second conclusion, which is to be expected, is that variable order methods are needed to handle a wide range of tolerances. A third and most interesting conclusion is that func- tion evaluations need to be relatively expensive to justify the extra overhead that multi-step methods require in order to keep down the number of func- tion evaluations. This last point was the key one in comparing extrapolation methods with Adams methods according to our earlier tests, while the Runge-Kutta methods were occasionally unreliable. However, the newer Runge-Kutta methods based on formula pairs appear to have displaced the extra- polation methods as the best ones to use when func- tion evaluations are relatively inexpensive. The corresponding situation with the results of STIFF DETEST is more complicated, and the con- clusions are not nearly so definite. For one thing, the Jacobian plays an important role, and an impor- tant factor in assessing different methods is the dimension of the problem, since two methods may differ in the relative number of n J and n opera- tions they perform. Also it may matter whether the Jacobian is determined numerically, whether it is sparse or not, and so on. Other important factors appear to be the nearness of the eigenvalues to the imaginary axis, the degree of non-linearity, the kind of coupling between the smooth and non-smooth components of the solution, and the tolerance. Also, methods are more likely to fail on stiff problems, and this makes the collecting of meaningful statis- tics more difficult. Nevertheless, it does appear that programs based on backward differentiation formulas stand up very well, their main weakness being when the eigenvalues are near the imaginary axis. But there are other good methods as well, and much more needs to be learned about them. For example, there are good methods based on second derivative formulas, implicit Runge-Kutta formulas, extrapolation, and block implicit techniques. 16-2 The future Of course it is easy to see that in the imme- diate future we can expect further testing in terms of reliability and efficiency of more programs over broader classes of problems. But I would like to make some brief comments on a number of specific areas in which I consider it would be especially worthwhile trying to improve on what we have been doing. One common feature in each area is, I believe, an attempt to focus attention on those aspects of testing that reflect the real needs of the ultimate user of ODE software, at least insofar as it is realistic to do so. As a first point, consider the interpretation of the tolerance . As a consequence of earlier arguments about error-per-step and error-per-unit- step, a number of us have become convinced that a very good way to interpret the parameter TOL is in terms of proportionality. The design goal for the ODE program is then to produce answers whose global error has a norm that is proportional to TOL. (This is intended to be quite independent of whether the error is relative or absolute, or of whether com- ponents of the error vector are weighted in some special way.) The factor of proportionality is of course problem dependent, but this interpretation is easily understood by a user, and also easily exploited by him. The current version of DETEST computes a measure of how smoothly the error depends on TOL, but improved measures need to be considered. It should be noted that it is often particularly difficult to meet this design goal for large TOL. It should also be noted that the frequency with which output is requested could have a considerable effect on how well some programs meet this goal. A second point that I would like to mention concerns the scale of a problem. I have often thought that it was not unreasonable to expect the user to provide some measure of the scale of his problem, as part of the problem specification. Requiring HMAX is a possibility, but HMAX has the disadvantage that it depends on the method as well as the problem, and the user should only be required to provide information about his problem. Many users would object to any such requirement (even if we had a good idea of what we meant by it!), but I would at least like to have the option of providing a value of SCALE, if only as a "knob" that can be turned one way to make the method more reliable (and the other way to make it more efficient of course) . In our earlier tests of Runge-Kutta methods we found some of them to be quite unreliable, until we imposed a maximum on the stepsize, in which case the improvement was dramatic. The test results can therefore depend in a crucial way on how this question is handled. Of course we can simply use default values, or build heuristics into our programs to estimate the scale of a problem. In any event, I believe that more attention should be given to this question, fit. should of course be acknowledged that this question is made more complicated by the fact that the scale may change during the course of the integration. ) What other measures of reliability or robust- ness should we test for? Can we test a program for Its ability to detect stiffness , or the presence of i i ■■' or.r Iri'il tJMB, DX trouble with roundoff? Can wo test separately its ability to get started ? Or to cope with discontinuities ? What about improper use ? (A common example of the latter is undefined input variables, but our programming languages seem to make it virtually impossible to deal with this problem in a decent way.) With regard to efficiency , it seems to me that more thought should be given to standardizing the way in which we present the results. For example, with non-stiff methods for a particular class of problems, I wonder if we could present the results in terms of a small number of parameters. One possibility might be the parameters h, f, t and k, where the computing cost has been determined to be h + f*F + t*TOL Thus h is a kind of overhead cost, f is the coeffi- cient related to the cost F of evaluating functions, while t and k appear only in the contribution due to TOL. It would be particularly helpful if some such result could be normalized in such a way that it was relatively independent of the problem class. Of course it would also be helpful if the normali- zation could factor out the dependence on a parti- cular machine, and, as well, the dependence on implementation details such as how the compiler handles subscripts. Other components of cost, such as storage, and other properties of the method, such as how well it gets started, or how well it copes with discontinuities, would also have to be considered, but it would be very helpful if these could also be presented in some sort of normalized way. The same sort of objective wuld of course be desirable in the case of stiff systems as well, but the situation would be much more complicated here. Nevertheless, we are already measuring the number of times an n J process is invoked, the number of n processes, and so on, so we have some of the ingre- dients . I continue to believe that theorems about pro- grams are important. I have in mind, for example, proving that a particular method solves all pro- blems in a particular class, in the sense that it produces z(x) such that ||z' (x)-f (x, z) || < TOL. I realize that such theorems can be established only over very limited classes, and are really not very practical. However, I continue to believe that establishing of such theorems is a good necessary condition to be satisfied and, moreover, that the weakness of the conditions under which such a theorem can be established is a measure of the strength of a method. With regard to broadening the class of problems over which methods are tested, I have a number of points in mind. One has been impressed on me because none of our problems ever uncovered the inability of Fehlberg pairs of order more than 4 to solve simple quadrature problems of the form y 1 = f(x). A second point is to have more problems in which stiffness changes. A third is to develop problems for testing the sensitivity of methods to changing parameters. (I would like to ask what we should do about the situation that often arises in chemical kinetics when, due to a minor numerical error, it is possible for one of the concentrations in Ih-i t ,ini ■ i k i ■ . 1 1 i v i • ami throw the calculation ofi 16-3 ibout this dil much worse than , purely because minor event Lly have .1 negligible effect on the final n mother factoi that Id be taken Into account with new cla prob ' Be- Mm to add new classes, we Ld ones. Surely, by should b. ■ Identify some problems undant and can be removed to make room I would like to conclude with a reminder that good documentation is needed and, if we are testing -, we should include some sort of assessment t documentation. Of course I don't mean to suggest that this can be quantified, but some kind of assessment would be helpful. Perhaps this would also involve a statement about how convenient the program is to use, about what options are available (that have not already been considered) . One particularly difficult option to deal with is the incorporation of stopping criteria other than the standard one of stopping at a particular value of the independent variable. Another property that must not be forgotten is the extent to which any special knowledge that the user might have can be exploited, regarding the stiffness of the problem for example, or some measure of the scale, or whether the problem is linear or not, or whether the user can provide the Jacobian, and so on. Altogether, it appears that our ODE solvers might be tested for many different things, and some of those properties that are already being tested might be done in ways that are at least somewhat more carefully defined and useful. 16-1 The NAG Library Runge-Kutta Routines by Ian Gladwell University of Manchester England A bstract The latest version of the NAG library includes a suite of subroutines based on the Runge-Kutta- Merson formula. This suite includes a main routine with a variety of features and a number of driver routines designed, in the main, for "black-box" use. These routines and a number of auxiliary ones are described in detail. 1 . Introduction The NAG library ordinary differential equation chapter contains a set of subroutines . based on the Runge-Kutta-Merson formula for integrating the initial value problem f(x, y), y(X) given (1.1) over a range [X, XEND,. At the lowest level there is a subroutine D02YAF which is designed to take one step of predetermined length using the Merson formula. This will be described in more detail in the next section. At the next level is a general subroutine D02PAF which is designed to integrate across a range by calling D02YAF at each step. D02PAF also has a number of other features such as interrupt and interpolation facilities and these will be described in detail in section 3. At a higher level still there are several driver routines for D02PAF. In the main, these routines have been designed so that the inexperienced user can solve his problems without recourse to D02PAF. However, one driver routine, D02BDF, is designed to calculate global error estimates and to perform stiffness checks, and this is intended for general use. All these routines are discussed in section A. Specifications of all the subroutines discussed here can be found in the NAG Manual, Mark 7 (1979). calculation of the error estimate was significant, componentwise. In designing this significance test we observe that the error estimate is calculated by componentwise subtraction of two approximate solutions at X + H (as is almost always the case in Runge-Kutta routines). We flag insignificance if, in this subtraction, all the significant digits of the two solutions cancel to rounding error level. Since the step H is determined (in routine D02PAF) from the error estimate we feel that it is important to be aware if this determining factor is completely dependent on the wordlength characteristics of the computer in question. In fact, the workspace array has been so arranged that all the workspace used on the intermediate steps in the Runge-Kutta step is used again in returning the information about the solution and its derivatives, and the error estimate and its significance. There are two other modes of ral lino rmtrinp D02YAF. Both are designed for use on calls from routine D02BDF. This routine computes global error estimates for solutions of routine D02PAF and attempts to determine the degree of "stiffness" of the system (1.1). In neither of these calls is there any need to retain information for inter- polation purposes and immediate storage savings are made. Also, when calculating global error estimates the stensize is controlled by independent calls to D02PAF so that the local error estimate is not needed, leading to savings in computer time and storage. When calculating the degree of stiffness, one of the algorithms used is essentially that of Shampine (1977) which requires a local error estimate for Euler's method with step H for use in a comparison with the Runge-Kutta-Merson error estimate. This is returned in a column of the workspace used on a normal call for returning interpolation information. 2. Subroutine D02PAF In its normal mode of operation routine D02YAF takes a step of length H from a current point X at which point the solution is stored in the array Y. It produces the solution at X + H in Y and it reorganises columns of the workspace array to retain values of y' at X + H and y and y' at X for use by the interpolation routines which are to be discussed in section 3 below. The routine D02YAF also returns a local the workspace and in •■in .in indicat Lor ' her the 3. Subroutine D02PAF and auxiliary interpolation routines Routine D02PAF is designed to integrate the system (1.1) from an initial point X to a final point XEND using a sequence of calls to routine D02YAF. The stepsizes in the calls to D02YAF are designed to control the local error, with a special short last step to reach XEND exactly. In many respects D02PAF is similar to well-known routines RKFA r ) (Shampine and Watts, 1976a) and DVERK (Hull et al., 1976). We will only report on the major points of di i i el ence between D02PAF and these rout Lnes. 17-1 Subroutine D02PAF has the calling sequence SUBROUTINE D02?M\\, SEND, N, Y, CIN, . FCN, I COMM, CONST. QOUT, W, IW, 1W1 , [FAIL) GES N, IV, IW1, IFAIL REA . YtH), CIN(b), TOL, COMM(5), CONST(3), i cout(U), w(iw, ivn RSAL FCN The parameters X and XEND are discussed above, Y is used to pass the computed solution to and from the routine, TOL is used for local error vol, FCN is the usual routine used to define the differential system, the array W is used mainly for workspace, and IFAIL is an error indicator to be discussed below. There are four arrays CIN, COMM, CONST, COUT used for communication with routine D02PAF, and in tain circumstances columns of the workspace array W are also used. In general terms, the array CIN is used to communicate information to control the mode of entry to the routine, to determine the type of error control, and to supply initial, maximum, minimum and continuation step- sizes when needed. The array COMM is used to handle a variety of facilities designed to inter- rupt the integration in different ways. If the integration is controlled largely through the use of interrupts with no intention to integrate to XEND, then the value XEND is essentially used to define a "scale" for the integration. Both CIN and COMM are described in more detail below. The array CONST can be used to pass control constants to the routine to assist in the internal stepsize strategy. Two of the elements of CONST can be used to vary the maximum and minimum ratios in attempted stepsizes on successive steps. (Of course the stepsize can be reduced as much as necessary at any step by trying a sequence of stepsizes of successively smaller magnitude.) The third element of CONST can be used to indicate if the differential system is linear and then a stepsize strategy is adopted to exploit the special properties of the Runge-Kutta-Merson formula for linear problems. The array COUT is used to return information to the user. Some of this information is statistical, for example maximum and minimum stepsizes used so far in the integration, but other elements of COUT are used to assist in interpolating the solution within integration steps, and yet others carry infomation between interrupts. This eliminates the need for use of COMMON but places the user on trust not to change elements of COUT at interrupts. The elimination of COMMON is important in a library environment as it permits integrating separate systems of differential equations simultaneously using the interrupt facility. The use of the four communication arrays (five including the workspace W) rather than just one, as in DVERK, is designed to permit future development in a library environment without major reprogramming effort. We now discuss the use of the control array CIN. The element CIN(l) is used to control the mode of integration. On initial entry, CIN(l) must be set to zero or one. If it is set to zero, default values of the other elements of CIN, and the elements of COMM and CONST are used. The routine D02PAF is designed to treat this choice as a special case and a fast integration from X to resul t s . ii ciN (l ) is set I all the other elements of CIN and the eli COMM and CONST must be set (the choice zei aull values in all cases). On return, by reaching XKND or by interrupt, a value of CIN(l) between two and six is obtained (depending on the mode of exit ), and the routine can be re-entered with this value. The choice CIN(2) = 0.0 (the default) gives mi itive/absolute error control on the maximum norm of the local error. Similarly CIN(2) = 1.0 and 2.0 give absolute and relative (with floor) local error tests respectively, and the choices CIN(2) = 3.0 and 4.0 give componentwise mixed and absolute local error tests respectively (a relative test is a special case of the mixed test in this case). In these latter cases, columns of the workspace array are used to pass componentwise error cont information. The element CIN(3) can be used to specify a minimum stepsize. However, a minimum stepsize designed to ensure progress across the range of integration is calculated and used if it is larger than the user's choice. (The minimum step- size is actually used only to activate error exits.) The element CIN(4) can be used to pass a maximum stepsize, the default being 0.5(XEND - X). The element CIN(5) can pass an initial stepsize estimate (and returns the initial stepsize used). If the estimate is smaller than the minimum or larger than the maximum then the corresponding minimum or maximum is used. If a default value is requested, then an estimate approximately proportional to (HyI|/||? II)*(XEND " X) * (MACHEPS * T0L) p ' 2 is used, where p is the order of the method. Safeguards are taken against zero values arising in the norms. This estimate usually gives a step rather smaller than the largest stepsize allowed by the local error control. This is deliberate as an iteration is employed to compute a realistic initial step and care must be taken to avoid nvprflnw in this iteration. The element CIN(6) is used only for output of the next stepsize in any exit from D02PAF. This stepsize is used on re-entry. The array COMM is used to enable interrupt control of D02PAF. If C0MM(1) > 0.0, it is used as an upper bound on the maximum norm of the solution as the integration proceeds. The integration is interrupted if this bound is exceeded (this facility is designed to enable users to avoid over- flows). If C0MM(2) > 0.0 it is used as an upper bound on the number of permitted calls to routine FCN, in the usual way. If C0MM(3) > 0.0, an inter- rupt is made on a change of sign of Y(I) - Z(I), I = 1, 2, ..., N where the user's chosen values Z(I) are passed in a column of workspace (this facility is useful for root-finding). If COMM (A) < 0.0, an interrupt is made when X - C0MM(5) changes sign (this simplifies output at user- selected points). If COMM(A) > 0.0, an interrupt is made after every step in the manner of many other solvers. All the other interrupts can clearly be accommodated by using COMM(A) at the expense of repeated calls to D02PAF, but we feel this is unjustified. Because of the interpolation process to be described below two steps are taken before permitting an interrupt, on the assumption that the user will wish to interpolate. Hence if an inter- rupt of the C0MM(1) > 0.0 or C0MM(3) > 0.0 type is required on the first accepted step we adjust the stepsize so that the interrupt comes on the second step instead. 17-2 As has been previously stated, D02PAF has been designed to permit interpolation to the approximate solution within integration steps. This interpol- ation process is designed to use the minimum amount of retained information whilst being of at least the same order of accuracy as the integration. To minimise the amount of retained information we use Hermite interpolation on three points giving an accuracy of 0(h) as against an accuracy of 0(h) in the integration. This choice clearly entails using solution and derivative information at X + H, X and X - H' where H' is the stepsize previous to the current one. Clearly we could save one column of workspace storage by not retaining, say, the derivative at X - H'. However, this would mean losing symmetry in the interpoaltion process and, in any case, the author feels that it is preferable that the integration error dominate the interpolation error so that by varying the local error tolerance TOL one varies the accuracy at all points approximately equally. We could instead have retained solution values at several earlier points but this would complicate the code far more than our choice. Note that the retained information is used passively, that is for inter- polation purposes only, and hence should have no effect on the stepsize strategy. This ensures that the algorithm remains one-step in character. However, it is insufficient just to consider the order of the interpolation error - its actual size is important. The interpolation error at a point t is approximately proportional to p(t) = (t - X - H) 2 (t-X) 2 (t-X+H') 2 . (3.1) If we consider max |p(t) te[X-H',X+H] function of r = H/H' then its minimum is attained when r=l and it increases monotonically and quickly as r increases taking the value 5 near r = 2 and values of over 100 for r < 5. The maximum value of ||p|| is attained on the longer of the subintervals [X-H*,x] and [X, X+H]. We assume that interpolation is generally only required on the step [X, X + H] and therefore the large values of ||p|| imply that we should restrict stepsize increases to bound the increase in inter- polation error from one step to the next. We choose an upper bound of r = 2 on this evidence. Two subroutines D02XAF and D02XBF are provided to perform the interpolation using information retained in Y, W and COUT on exit from D02PAF. Routine D02XAF computes an approximate solution at any point in [X-H', X + H], whilst D02XBF computes a selected component of the solution in the same range. It was felt that to supply routines for these two extreme cases was justi- fiable but that any more complicated selection of components might be more expensive to calculate than repeated calls to D02XBF. It is widely accepted that the error in integrating the system (1.1) should be approx- imately proportional to TOL. To ensure this is the case, routine D02PAF, in common with most other modern subroutines, is designed to vary the stepsize to make the local error vary proportionally and smoothly with TOL. One asepct of this strategy is to avoid choosing stepsizes on the evidence of insignificant error estimates (see section 2). These insignificant estimates can arise in two ways. During an integration, the error estimate is unlikely to become insignificant unless TOL is so small that the stepsize becomes very small. If this occurs an error exit is taken (other subroutines have error exits for the related case that too much accuracy is requested). Alternatively an insignificant error estimate can be returned by D02YAF during the initial stepsize iteration. This is just taken as a sign that the current stepsize is too small and the stepsize is increased and the iteration proceeds. An indication of the dependence of the local error on TOL is given through two elements of the array COUT. A count is kept of the total number of successful steps and the number using the maximum stepsize. If the second count is large relative to the first then one might infer that this choice of TOL has not controlled the local error. 4. Driver Subroutines There are five driver subroutines for D02PAF. The first and simplest, D02BAF, integrates the system (1.1) across the range [X, XEND], It is intended for very inexperienced users and hence has a particularly simple calling sequence designed to define the problem, input a local error tolerance TOL, designate some workspace and permit error exits. The second routine, D02BBF, is similar to D02BAF but implements slightly more general error control (similar to the use of CIN(2) = 0.0, 1.0 or 2.0 in D02PAF), and permits output at user-selected points. This output is implemented through a user- supplied subroutine OUT. This routine is called at the initial point X by D02BBF and a parameter of OUT must be changed to the next output point and so on. Hence, when OUT is called at the current output point, it must redefine the output point in addition to performing its printing task. Internally D02BBF uses the interrupt C0MM(4) < 0.0 to pass output points to D02PAF. There are also two driver routines designed to solve equations posed in terms of the solution) as the integration proceeds. The simpler of these, D02BGF, solves an equation yf(x) = v for x and exploits the interrupt C0MM(3) > 0.0 to locate the position of x approximately. The second, D02BHF, solves a general equation g(x, y) = by locating a change in sign of g using the interrupt C0MM(4) > 0.0 (interrupt at every step). In both cases the equations are solved to machine precision using the interpolants defined by D02XBF and D02XAF respectively and using a secant/bisection rootfinder C05AZF based on an algorithm by Bus and Dekker (1975) The equations are solved to machine accuracy to put all the emphasis on the accuracy of the integration, so that by varying the local error tolerance TOL we can simply control the accuracy of the solution of the equations. The last driver routine, D02BDF, is designed to solve (1.1) whilst producing a global error estimate and a variety of related statistics, and to give an indication whether the system (1.1) is stiff. The global error estimate is calculated using the two- and-one technique discussed in Shampine and Watts (1976b). The actual solution output should be close to that which would be calculated by a call to routine D02PAF with the same value of TOL (though it is calculated by calls to D02YAF only). The calculation proceeds by interrupting D02PAF at each 17-3 and taking two hall stops with DOJYAF to calc- ulate .in auxiliary solution, the difference between the two solutions giving the estimate. The other statistics produced include the maximum error estimate in magnitude in the integration for each t .inJ information about the monotonic growth of the error estimate .iiui any sign changes in the estimate. The user is permitted to specify a maximum size for the solution and this can be used to prevent overflow during calculation of the various solutions. Overflow is more likely to occur during calculation of the solution using DC2YAF as there is no error control on this solut i on. SHAMPINE, L.F. & WATTS, H.A. (1976a) Practical Solution of Ordinary Differential Equations by Runge-Kutta Methods, Rept . SAND 76-0585, Sandia Labs., Albuquerque, New Mexico. SHAMPINE, L.F. & WATTS, H.A. (1976b) Global Error Estimation for Ordinary Differential Equations, TOMS, 2, pp. 172-186. The stiffness check is based on two separate calculations. The first uses the ideas of Shampine (1977) where an algorithm is described for a Runge- Kutta-Fehlberg method. The stiffness check is based on a comparison of stepsizes predicted from the Fehlberg algorithm and Euler's method. The comparison is dependent on certain properties of the absolute stability regions of the two methods. It is fortuitous that the stability region of the Merson algorithm used in D02PAF are as well- suited to this application as the Fehlberg stability region. The second stability check is based on the simple observation that if a stiff problem is integrated using a non-stiff solver such as D02PAF, the stepsizes used are, in the main, controlled by the requirement of stability rather than by the local error tolerance TOL . Hence, if separate calls to D02PAF are made with local error tolerances TOL and 2P * TOL, twice the number of steps should be used in the first case if the step- size is being controlled by the local error tolerance. We compute a stiffness check by measuring the deviation from this expected behaviour. In our tests we have found these two checks approximately equally reliable (but the second is more expensive to compute and hence the user may opt to avoid it). In fact we have also tested other possible stiffness checks. For example, we have counted the number of unsuccessful steps relative to the number of successful ones for various values of TOL. One can formulate theories about the behaviour of this ratio for stiff systems but unfortunately this behaviour is not restricted to such systems and the ratio does not provide a reliable measure. REFERENCES BUS, J.C.P. & DEKKER, T. (1975) Two Efficient Algorithms with Guaranteed Convergence for Finding a Zero of a Function, TOMS ,1, pp. 330-346. HULL, T.E., ENRIGHT, W.H. & JACKSON, K.R. (1976) Users Guide to DVERK - A Subroutine for Solving Non-Stiff ODE's. Department of Computer Science, University of Toronto, Tech. Rept. 100. NAG Manual, Mark 7(1979) NAG Central Office, 7 Banbury Road, Oxford. SHAMPINE, L.F. (1977) Stiffness and Non-Stiff Differential Equation Solvers II: Detecting Stiffness with Runge-Kutta Methods, TOMS , 3, pp.44-53. 17-4 TOWARDS ROBUST GENERAL PURPOSE SIMULATION SOFTWARE Francois E. Cellier Peter J. Moebius Institute for Automatic Control The Swiss Federal Institute of Technology Zurich ETH - Zentrun CH-8092 Zurich Switzerland The aim of this paper is to present methods to improve the robustness of a general purpose run-time system for digital simulation of continuous systems. This aspect of simulation software robustness is just one of several possible aspects. One could as well discuss the robustness of a language to formulate simulation problems (simula- tion language) or the robustness of a compiler for such a language (simulation compiler). Such aspects will, however, not be considered in this article. The run-time system robustness is a very important aspect which has not sufficiently been taken into account in the development of existing simulation software. The aim of this paper is to specify demands rather than to present optimal solutions. Solutions are presented from an engineering point of view. They improve the run-time system robustness to some extent, but better solutions can certainly be found. It is hoped that some of the Numerical Mathe- maticians attending this Conference may find' these problems of interest and may develop better solutions to them in the future. This is one of the main reasons why the authors de- cided to present this material at a Conference of Numerical Mathematics rather than at a Conference of Simulation Techniques. I) INTRODUCTION: A simulation run-time system can be robust in two senses : produced results. For this reason, it is vital that each algorithm in the system has its own "bell" which rings as soon as it is unable to properly proceed. Under no circumstances are incorrect results allowed to be displayed to the user. a) The user should never be required to provide any kind of information which he does not have at his disposal. He should be able to concentrate on those factors which have to <^o with the statement of his problem, and should be relieved, as much as pos- sible, of all aspects which have to do with the way his problem is executed on the machine. He should be able to describe his system as easily as pos- sible in terms which are closely related to his common language, but must not be required to pro- vide a step-size for the numerical integration or to specify the integration algorithm to be used. b) The run-time software itself must be able to check whether the produced time responses are "correct" (within a prescribed tolerance range). The user, normally, has a more or less precise (although often not mathematically formulated) knowledge of the system he is investigating. He has, however, hardly any "insight information" into the tool he is using for that task. He is, usually, very credulous (the obtained results must be cor- rect because the computer displays 14 digits!), and he has no means to judge the correctness oi thi II) AUTOMATED SELECTION OF INTEGRATION ALGORITHMS: A huge step towards robust simulation software has been taken in the development of step-size controlled integration algorithms. Before these algorithms existed, the user of digital simulation software was required to supply information concer- ning the step-size to be used — an information item which he clearly did not have at his disposal. Now, the user can simply provide a tolerance range for the accuracy of the results. This is identical to requesting the user to identify the smallest number in his problem which can be distinguished from zero. This question can certainly be answered by any user, independently of whether he is an expert in Numerical Mathematics or not, since it is closely related to the physics of the problem, and not to the numerical behaviour of the algorithm. Available simulation software, up to now, usually offers a comprehensive selection of different inte- 18-1 tion algorithm*. It does, however, not toll the : which would be the most appropriate one for his particular application. In this way, the user igaifl confronted with making a decision on some- thing he does not really understand. Experience has ~-n that the majority of the average users always operate with the default integration method imple- mented in the package which, in most cases, is a Runge-Kutta algorithm of 4th order. Since he does not know what to specify, he simply ignores that question, and after some time of using the software he has even forgotten that the language provides him with the facility to select among different integration algorithms. So far, no integration algorithm could be found which would be able to handle all kinds of problems equally well, and it is more than doubtful whether such an algorithm could be found at all. The user, who does not make use of the facility to select among different integration rules, will, consequently, often waste a lot of computing power. Although much research has been devoted to the development of different integration methods for the different classes of application problems [3,8], the user has, however, no means to easily determine, from the state space description, the problem class to which his par- ticular application belongs. For this reason, the selection of the appropriate integration algorithm should also be automated. For this purpose, we try to extract features from an application problem during its execution which are supposed to characterize the numerical be- haviour of that particular problem as completely as possible. These features are then combined in a feature space in which we can identify specific clusters for which a particular integration method is optimally suited. The proposed methodology for the solution to this problem originates from pat- tern recognition. What features may be used for this purpose? A first feature can be associated with the accuracy re- quirements for the problem. It can be found that the CPU-cost to execute a particular problem depends on the required relative accuracy. Low order algorithms are appropriate for the treatment of systems from the "gray-" and "black box" area where the available data and models are so vague that a precise numerical integration does not make much sense, whereas higher order algorithms are ap- propriate for the handling of systems form the "white box" area, e.g. from celestrial mechanics. Since the user is requested to specify the wanted relative accuracy, this feature can be extracted from the input data. Multi-step methods require more CPU-time during their initial phase, but are more economic than one-step methods if integration goes on over a longer interval of simulated time. This can be ex- plained by the fact that one-step methods are self-starting whereas multi-step methods need to be initialized. This leads to a second feature. Since the integration has to be restarted after "event times" (instants at which some state trajectories are discontinuous), multi-step integration is in favour for purely continuous problems or for problems with few event times, where. is one-step integration is appropriate Tor combined (con- tinuous/discrete) problems where discrete events occur with a high density. Each integration algorithm has associated with it a domain of numerical stability. The stiffer a parti- cular problem is (the more the different eigen- values (A.) of the Jacobian are separated), the smaller the step-size (h) must be in order to keep all (X.*h) within the stability region of the algo- rithm. Fortunately, special integration algorithms could be found for which this restriction no longer holds (A-stability, stiff stability). If such an algorithm is used, the step-size need not be re- duced due to restrictions imposed by stability de- mands, but is determined exclusively by the re- quirements of accuracy. For this reason the eigen- value distribution of the Jacobian determines a third feature which must influence our decision as to which integration algorithm to use for the execution of a particular problem. / These three features can now be combined to a fea- ture space as depicted in fig.l. i t Accuracy Stiffness Smoothness Fig.l: Feature space for selection of integration algorithm 18 clusters have been distinguished in fig.l, and fig. 2a and 2b show integration rules which can be associated with them. 1 I Accuracy i Aco AdaiK (high order Adams (medium order) Adams {high order) Adam [ton order) (In- order) Runge- Kutta- Fehlberg (8th ord.) Rungc- Kutta- FeMberg (Sth ord.) FeMberg (8th ord.) V a *e (i/n) where a is the largest singular value of the matrix, £ is the machine resolution (e.g. MO ' "* on a CDC 6000 series installation), and n is equal to the order of the model. For higher order models e" n approaches 1.0 and hardly any eigenvalues will then be properly computable. Smaller eigenvalues can take any value and small modifications of the elements of the matrix can place them almost any- where within the band of incertainty. As can be seen, the system becomes periodically un- stable. During such time intervals, errors will ac- cumulate, and, consequently, we must be careful in the interpretation of the obtained results. This fact should be reported by the software to the user. For this purpose, we define another feature (stability). A variable STAB is set equal to zero when all eigenvalues lie in the left half-plane and is equal to one as soon as at least one of the eigenvalues moves into the right half-plane. STAB - { 0.0 : Re{X.} < 0.0 ; i=l, 1.0 : otherwise We now collect statistics on STAB as for time-persistent variables. In this way we obtain the integral of STAB over time divided by the run length: 1.0 FF5 "fin beg 1 J fin (STAB)dt beg The fifth feature (FF5) is a real number between 0.0 and 1.0. If it is close to 0.0, the results ob- tained by simulation have a good chance to be re- liable. If it is close to 1.0, the obtained results are most probably nonsense, and they must be cau- tiously verified. If we now assume that the matrix under investiga- tion is a Jacobian of a state space description for a real physical process, then the elements of the Jacobian are extracted from measurements, and cannot be computed more accurately than e, which is a relative accuracy of measurement. We, therefore, must assume that, within that relative accuracy e, the elements of the matrix can take any values. In this case, we must also assume that eigenvalues of the Jacobian which do not fulfil the more stringent inequality: > a *e 0/n) can take any value within that broader band, although they can be much more accurately computed as soon as any particular values have been assigned to all elements of the Jacobian. This means that as soon as there exist eigenvalues for which the second (more stringent) inequality does not hold, small variations in the systems parameters which lie within the inaccuracy of the measurement can make the model non-stiff or stiff or even unstable. Physically seen, these eigenvalues correspond to merely constant modes which could as well be eliminated from the equation set resulting in a model reduction. Numerically seen, these eigen- values can lead to accumulation of errors so that these modes can drift away over a longer span of simulated time, again resulting in incorrect simulation trajectories. Together with the eigenvalues, we compute the fol- lowing quantity: V) VALIDATION OF THE MODEL WITH RESPECT TO THE SYSTEM UNDER INVESTIGATION: Another non-trivial question is whether a model, BORD and the number k indicating those eigenvalues whose absolute value is smaller than BORD: 18-4 n : | A . | >_ BORD k = 2C(j.) 5 j. = { i=l 1: I A. I < BORD k represents an integer between and n. We now collect statistics on the quantity (k/n) as for time-persistent variables, and obtain a sixth feature: FF6 1.0 "fin "beg "fin * f(k/n)dt t,-. - t, -* beg Also the sixth feature (FF6) is a real number bet- ween 0.0 and 1.0. If it is close to 0.0, the model has some chance to be valid. If it is close to 1.0, the model is most probably invalid, and it should be further investigated. Evaluation of features FF5 and FF6 requires com- putation of the eigenvalues of the Jacobian once per integration step. Since this can be expensive, it should not be done automatically, but the user must have a switch at his disposal to turn computa- tion on and off. In this way he can use these fea- tures during the development of a new model, whereas he can turn computation off during produc- tion runs. VI) DETERMINATION OF CRITICAL STATES: In section V we have discussed the case where single eigenvalues were situated close to the ima- ginary axis, and we have seen that in such a case it may be possible to reduce the order of the model . It is, however, as interesting to discuss the oposite case where single eigenvalues are located in the A-plane far to the left. We call these modes the "critical states" of the system. Very often one is not really interested in these fast transients. In such a case one could eliminate these modes from the equation set. If the fast transients are impor- tant one could at least try to utilize special integration techniques (like using singular pertur- bations) to expedite integration. One can, of course, again compute the eigenvalue distribution for the solution of this problem. However, it is not always easy to see which state equations are responsible for such an eigenvalue. For this reason we recommend the following pro- cedure. We reserve an integer array i of length n which is initialized to zero. Each time an integration step has to be rejected due to accuracy requirements not met, we add 1 to each element of the i is bounded for all n . We prove that if a k-step method is exponentially fitted to +iw (w real) of order p >^ and is A-stable, then the method is implicit and p <_ 2 . This extends the well-known result of Dahlquist (1963) which states that if a standard k-step method of order p is A-stable, then it is implicit and p _< 2 . We show that if an exponentially-fitted method is exact on (and only on) the functions 1, e V X t , P (the X. need not be distinct) then its local truncation error x satisfies hi = c h P+1 Try + 0(h P+2 ) n Try(t) = D(D - X x ) (D - X )y(t) , We then prove that the usual Milne device can be used to estimate the local truncation error of an exponentially-fitted predictor-corrector pair whenever both the predictor and the corrector are exact on the same set of functions. A k-step method is called A-stable if, when- ever it is applied to an ODE of the form y = Xy with ReX < , it produces a numerical solution u (h) which for n 19-1 The Linear Algebra of Discretisation Methods S. McKee Computing Laboratory and Hertford College, Oxford where f (t , y (t ) ) i function on 0Nt|)>i(,K+0 For example this arises when product integration techniques are applied to solve the Abel's equation (Holyhead and McKee) . Thus this paper includes a study of A in this more general formulation . Dahlquist stability and absolute stability, as defined by (4) and (5), can be characterised by requiring that the eigenvalues of a product of Frobenius matrices all lie inside or on the unit circle, those on having linear elementary divisors. They can be characterised by the approach of Donelson and Hansen (1971), and by the generating function approach of Holyhead and McKee (1976). This paper demonstrates the precise equivalence of these and other different characterisations. These results have evolved from and in parallel with earlier works by the author and his co-workers. A list of the more pertinent references can be found in the bibliography. Blbl iograph Allen, K and McKee, S. Discretisation methods for delay differential equations (in preparation). Andrade, Celia and McKee, S. On optimal high accuracy linear multistep methods for first kind Volterra integral equations, BIT (to appear). Donelson, J and Hansen, E. (1971), Cyclic composite multistep predictor- corrector methods, SINUM, 8_, No. 1, 137-157. Holyhead, P.A.W., McKee, S. and Taylor, P.J. (1975). Multistep methods for solving linear Volterra integral equations of the first kind, SINUM 12, No. 5, 698-711. Holyhead, P.A.W. and McKee, S. (1976). Stability and convergence of multistep methods for linear Volterra integral equations of the first kind, SINUM 13, No. 2, 269-292. Holyhead, P.A.W. and McKee, S. Linear multistep methods for the generalised Abel's equation (submitted to SINUM). McKee, S. Cyclic multistep methods for solving Volterra integr o-di f f er ential equations, SINUM (to appear). McKee, S. Best convergence rates for linear multistep methods for Volterra first kind equations, Computing (to appear ) . McKee, S. The analysis of a variable step, variable coefficient linear multistep method for solving singular integro-dif f erential equations arising from the diffusion of discrete particles in a turbulent fluid. JIMA (to appear ) . McKee, S. The linear algebra of discretization methods (in preparation). Acknowledgement s The author is the C.E.G.B. Research Fellow in numerical analysis at Oxford and is grateful for their financial support . 20-2 Some Experiments with STRIDE Fred Chipman Acadia University 0. Introduction STRIDE is an integrator for systems ot ordinary differential equa- tions, designed by Kevin Burrage, J.C. Butcher and Fred Chipman. In section 1, a brief description of the method will be presented, followed by a few details of the algorithm. Section 2 presents prelim- inary results of experiments with STRIDE on two problems. These results seem sufficiently promising to justify further work on this type of implicit Runge-Kutta method. 1. Description of the method The algor- ithm implemented in STRIDE is based on the singly- imp lici t methods introduced in Burrage [3] and Butcher [5]. These meth- ods are of the implicit Runge-Kutta type but with coefficient matrix so chosen as to have only a single distinct eigenvalue. When the method of implementation des- cribed in Butcher [4] is used, methods with coefficient matrices of this type require a similar number of arithmetic operations for the performance of its LU factorization as do linear multistep methods . In the usual notation for Runge- Kutta methods, the stages are given by (1.1) Yi=y n . 1+ h I a i .f(x n _ 1+ hc j ,Y j ) j = 1 (i=l,2,...,s) and the result at the end of a step by s Ci. 2) y n =y n -i +h .I b j f ( x n-i +hc j> Y j) • Let £i,...£ s represent distinct zeros of L s (z), the degree s Laguerre polynomial". Then, with A any positive number, the coefficients for the particu- lar methods we are considering are gi^cn Ci=A£i, i=l...s, ind b T =f ] ,...]/, )\ L , whci i .J/ j J . latri) i ned , has characteristic polynomial (z-A) s , and, as shown in Butcher [5 ] , can be written AT = A 1 . 1 - •-1 • 1 where, for T=(t i j) and T" 1 =(tT 1 ), t^-LjUi) and ^-- Itll^llp If 'i= Y i-y n -i. h =f(x n . 1 +hc i ,Y i ), (i=l,2,...s) and Z-, F. denote the corresponding transformed quantities given by CI. 3) Z i= I tTjz. , Fi- I t^F,, j-l ij J j = l ij J and if H=hA, then (1.1) is equivalent to the system (1.4) HF\ Z i = HF i -HF i _ 1 (i=2,...s) and (1.2) is equivalent to s 1 (1.5) y n =y n - r e I ±LJ(e)HF., where e-i/x. i = l In the special case that 6=Cj for some i, and hence c^=l, the computed result y n is equal to Yj_, and the method is of order s. This choice is the one made in STRIDE and, in consequence, the stability proper- ties of the methods are acceptable for many stiff problems. Since the expansion of the solution given by (1.5) is equally valid for other choices of 6, this formula serves as a basis for interpolation of output values and for the provision of starting values for the iterative evalua- tion of Y.,Y-,...Y . For stiff problems, (1.4) is solved using a modified Newton iteration. This requires the LU decom- position of the matrix (1-11. J), where J I notes the .lacobian matrix for the 21-1 .rential equation. As in most algor- ithms, J is not computed in Bach iteration, but only when: convergence rate is unacceptable, (b) H has changed by a factor greater than 2 or less than .S since the last evaluation, 20 steps have ellapsed since the last evaluation. This iterative scheme is modified in STRIDE bv the use of a relaxation factor I+H IT), where H' is the value of H at which (I-HJ) was last evaluated and :ored. The rationale behind this mod- ification is that, if u is an eigenvalue of J then the corresponding eigenvalue of an iterative scheme based on R(I-HJ)" 1 is 1- R( 1-Hu)/ (1-H 1 u) , and the supremum of this for all u in the negative half-plane is minimized when R is given the value we use. This minimum value is in fact equal to | (H-H ' ) / (H+rT ) | < 1 so that, for slowly varying Jacobians, eventual con- vergence would be achieved no matter how much H varies. For non-stiff problems (1.4) is solved using a fixed point iteration, applied either a fixed number of times or iterated to convergence. Error estimates are provided by em- bedding the s-stage (order s) method in an s+1 stage (order s + 1) method. As a bonus, error estimates are also provided for methods with orders from 1 up to s (and, if H has remained unchanged for two steps, up to s+1) . Thus we are able to implement the methods in a variable order, variable step size manner. A complete description of the gorithm STRIDE is given in [1], wh the theoretical aspects are presen [2]. Only a few basic parameter d tions will be given here. The pro STRIDE has only the four arguments FN, JAC and VALUES. N is the dime of the problem, FN and JAC are pro for the evaluation of f(x,y) and J and VALUES is a user defined proce by way of which all information is changed between the method and the Amongst the parameters of VALUES i the integer variable METHOD. The digits nin 2 n 3 of METHOD indicate t STRIDE the following: n,=l = - error per step is to be con error per unit step is to b trolled Newton iteration with an in approximated Jacobian Newton iteration with a use supplied Jacobian Fixed point iteration, iter convergence fixed point iteration appli fixed number of times, n, indicates, for n 2 *0, the maximum n 2 = 3 = 2 = 1 = al- ile ted in escrip- cedure N, nsion cedures (x,y) , dure ex- user. s three o trolled e con- ternally r ated to ed a number of iterations permitted whereas, for n it indicates the number of iterations to take place. Finally, in incorporating the Lncre ments to both X and Y, the effect of round off error is reduced by computing resi- duals equal to the intended increim minus the actual changes in the quantities being updated. These residuals are added to the increment for incorporation .into the following step. Hence, if < i the increment in Y, we perform the following steps in updating Y: RES:=RES+AY; TEMP:*Y;Y:=Y+RES; RES: = RES- (Y-TEMP) ; 2. The Experiments In this section we describe several experiments performed with STRDIE using the following problems. Problem 1: The equations of motion for a satellite moving under the influence of the earth and the moon in a coordinate system rotating so as to keep the posi- tions of the earth and the moon fixed [6]. y" = 2y 2 +yi-u' (yl+u)/r?-u(yi-u')/rl y" = -2y'i+y2-u'y 2 /r^-uy2/rl yi (0) = 1.2, y!(0)=y 2 (0)=0, y 2 (0)=-l. 0493575098. . . where r. = [(yi + u) 2 + yi]\ r 2 =[( yi -u<) 2 +y 2 2} h u = 1/82.45 and u' = 1-u. These initial conditions produce a closed orbit with period T=6 . 19 2169 3313 .. . The solution is required over one period. Problem 2: A Chemical kinetics problem [7, pp 268-269] yi = - .04y 1 + 10'*y2y3 y 2 = .04y 1 -10'*y2y3-3-10 7 yi y 3 = 3-10 7 y! with y,(0)=l, y 2 (0) = y 3 (0) = 0. Results are required at t=4 # 10 , k=0,l,...,9. The first problem is non-stiff, and is a good test of the stepsize control. The second problem is stiff, and is fre- quently used as an illustrative example. Experiment 1: To investigate the relative efficiencies of the non-stiff options, problem 1 was solved using METHOD = 101, 102, 112, 113 with EPS=10" 5 , 10 -8 . The number of function evaluations required were 21-2 IODj LOJ 112. 113 D\ DG ; 1070 L882 1882 1882 846 10" 7 1687 3165 3548 3193 1288 These results indicate that for this problem there is little use in doing more than one fixed point iteration. This has proved to he the case for all other non- stiff problems that have been tested, and thus METHOD=101 would seem to be a natural choice for such problems. It should be noted however that this method required considerably more evaluations than Krough's DVDQ. Experiment 2: To test the effectiveness of the accumulation of residuals to mini- mize round off error, problem 2 was solved with METHOD 101, and a realtive error tollerance of 10" 8 . The results, at the end of one period were: v i Y2 >'3 yif R 1.200000" 3.08(-6) -2.58(-6) -1.0493583 R 1.2000012 4.48(-6) -2.60(-6) -1.0493587 The results show a slight improvement with the acculumlation of residuals (R) . Experiment 5: To compare the effects of restricting the Newton iteration to a max- imum of 2, 3, or 4 iterations, problem 2 was solved with METHOD=122, 123, 124 and EPS=10~6. The following are the number of function and Jacobian evaluations at t=4xl0 9 : METH0D= 122 #FN UAC 12 59 50 125 1161 38 124 1199 EPISO DE 620 105 a maximum of Included here For this particular problem, three iterations seems best. are results using Hindmarsh's EPISODE, which required far fewer function evalua tions, but considerably more Jacobian evaluat i on . Experiment 4: The J whenever any one of described in section was found by monitor solving problem 2 th convergence rate ace of the evaluations, large change in H ac remainder. Letting value of H for which tluated, (I-HJ) is H ' /H no longe r ] i In thi n t , ■ HOD=133 and i gned the 11/2,2) and [1/3,3] , ■Jacobian was used function < of evaluations of bo t j -4 . 10 i [3/4,4/3] [2/3,3 acobian is evaluated the three conditions lis satisfies . I t ing the program while at an unsatisfactory ounted for about 1/3 and a sufficiently counted for the H' represent that (I-H.JJ was last re-computed when i n the interval [ h , 2] problem 2 was solved EPS=10" 6 . I was [3/4,4/3], [2/3,3/2], and a numer i en 1 that the number of I the number th f and J. Results ', i-0,|,...i. /2| [J/2,2] |1/3,3| 1 3 36 2 49 3 672 4 820 5 963 6 1059 7 1135 8 1196 9 1229 298 525 701 837 968 1063 1170 1225 1287 375 550 738 934 1056 1159 1273 1324 136 7 3 39 526 968 1168 1310 1391 1456 1548 For this problem, the two smaller inter- vals seem more appropriate. Experiment 5: To determine the effective ness of the relaxation factor R in the Newton iteration, problem 2 was solved using METHOD*] 23, EPS=10" 5 , 10" 8 and R=2/(l+H/H') and 1 (no relaxation factor). The results are given at t=4xl0 5 . The number of function evaluations (#FN), Jacobian evaluations (#JAC) , divergent iterations (#DIVGD) and non- convergent iterations (# CONVGD) were: EPS=10 J R=2/(l+H/H') R=l EPS=10 R=2/(l+H/H') R=l #FN 675 832 1975 2243 #JAC 27 33 45 55 #DIVGD 5 3 7 12 # CONVGD 3 13 12 20 For both values of EPS substantially more work was required in the case where no relaxation factor was used. [1] Burrag F., inte tion Repo appe [2] "An ment tics (to [3] Burrag R ung s'tif 18, Referenc es e , K. , Butcher, 'STRIDE: A stab grator for diffe s", Computationa rt ; Universi ty o ar) . efficient Runge ation" , Computat Report , Uni vers appear) . e , K. , "A specia e-Kutta methods f differential e 22-41. J . C . , Chipman , le Runge-Kutta rential equa- 1 Mathematics f Auckland (to -Kutta imple- ional Mathema- ity of Auckland 1 family of for solving quations" , BIT [4] Butcher, J.C., "On the implementation of implicit Runge-Kutta methods", BIT 16, 237-240. [5] "A transformed implicit Runge- Kutta method", computational mathe- matics Report No. 13, University of Auck land . [6] Krogh, F.T., "On testing a subrou- tine for the numerical integration of ordinary differential equations" JPL Technical Memorandum No. 217. [7] Lapidus, L., Seinfeld, J. II., "Numeri- cal solution of ordinary differen- tial equations", Academic Press, N.Y. , 1977. 189 202 21-3 EXPONENTIAL INTERPOLATIONS AND THE NUMERICAL SOLUTION OF O.D.E (SUMMARY) ARIEH ISERLES King's College, Cambridge. Exponential approximations, mainly the Pade approximations, N-approximations and exponentially fitted approximations, have attracted great interest during recent years as a means for the numerical solution of ordinary differential systems. The common feature of these approximations is that they fit to the exponential and to its derivatives at the origin. Hence, they are, in fact, Hermitian interpolations and not approximations. The high degree of interpolation at the origin is of great benefit if (in the case of either linear or mildly nonlinear autonomious differential systems) the behaviour of the solution is determined by a single principal eigenvalue of the Jacobian matrix. From the computational point of view, this is not the case if parabolic or hyperbolic partial differential systems are solved by the method of lines. The solution of the derived linear ordinary differential systems is determined (unless the solution is in an asymptotic stage) by a great number of eigenvalues, in the neighbourhood of the origin. Even if the exact spectrum of the matrix is unknown, it is a useful practical assumption that the eigenvalues are more or less spread uniformly in an interval [-T, 0] . Out of these eigenvalues, only a smaller subset, contained in [-T , 0] , < T « T, influence strongly the numerical solution of the system. Hence, it is of practical interest to investigate the use of exponential approximations which fit to the exponential along an equispaced mesh in [-T ,0] , instead of high-degree interpolation at the origin. Such exponential interpolations are derived in this paper. The first result extends the C-polynomial theory of N^rsett to the interpolatory case: A 1 k Theorem 1: If p(x,c) = £ — , p (c) x , where k=o c > 0, p (c) - 1 and p (c) , p ,(c) are ' r n r o r n-l arbitrary functions of c and (1) R(x,c,p):= Y, (-l) k D n " k p(0,c)(l-e" C ) k (- x/c ). k=o x ; £ (-l) k D n " k p(l,c)(e C -l) k (- x /c). k=0 x k where (-y). is the factorial function, (-y) = 1, rC o (-y) k = (-y)(-y+l) (y+k-1) for k ^ 1 and D p(a,c) is the (n-k)-th derivative of p(x,c) at (a,c), with respect to x, then R(qc,c,p) = e ■qc $ q N < n . Definition : The function p(x,c), as defined in Theorem 1, is called the C-function of the exponential interpolation (1). If p(x,c) is continuous for c > o (or an interval < c < C) and lim p(x,c) exists, the C-function is said to be c-»o+' smooth. Theorem 2 : If p is a smooth C-function, let p(x) = lim p(x,c) and let R(x,p) be the c-» o+ exponential approximation which is generated by the C-polynomial p(l-x). Then, if for all suffici- ently small c, and for a positive integer s, — Q C R(qc,cp) = e , O^q^n+s, then ~ , ,». -x f] , n+s+1) R(x,p) - e = U (x Theorem 3 : Every rational function R (x) = P (x)/Q (x) , that satisfies the conditions deg P c ( n, deg Q c i. q 4 n, can be expressed in the form (1), and n, and R (qc) = e -c l c , 0, P (x,c,a):= Zh ( £ (-l) P ( k )(l+acp)V pC )(--), n , k. *-" p c k k=o p=o and R (x,c,a):= P (x,c,a) / (l+ax) n n n then R (qc.c.a) = e , Ojcqxn n Further, if a satisfies o »-, , , ,k,n+L ., , .n -kc J] (-1) ( . )(l+a ck) e = , k=o then R ((n+1) c ,c ,a) = e ("N- interpolations" ) . n * A-acceptability : As yet, no general conditions have been found that imply A-acceptability. Various computer results are available, concerning the range of c which yields A-acceptability and A -acceptability for various Pade and N- interpolations. Arieh Iserles King's College Cambridge CB2 1ST England. 22-2 On Che Existence of Maximum Polynomial Degree Nordsleck-Gear (k,p) Methods for All (k,p) Roy Danchick and Mario L. Juncosa The Rand Corporation Santa Monica, California 90406 Introduction This paper sketches a fundamental theorem on the existence ot Nordsieck-Gear (k,p) methods of maximal polynomial degree. The notation (k,p) refers to the number of scaled derivatives re- tained by the method (k) and the order of the differential equation (p) . The method's poly- nomial degree is the highest degree polynomial that can be numerically integrated exactly; given exact starting values. In [1] it was shown, for fixed stepsize and for sufficiently accurate starting values, that if a certain sequence of corrector-predictor coeffi- cients was well-defined, then the Nordsieck-Gear (k,p) methods could obtain a maximum polynomial degree of k+1. The sequence was shown to be well- defined in the (k,l) case for k=2,3,A, and 5. These results led to the conjecture that such a well-defined sequence exists for all (k,p). The results of this paper prove that conjecture. Basic Definitions I is the (k+l)x(k+l) identity matrix ,(n+l) r ,(n+l) * ,'- v ~ 2 , ^"■i ! * ^«»«>^i>-40« if < i < j < k, if < j < i < k, (t 1 ; 1 )]^^), ( k ; l ),...,( k ; l )] T . Preliminary Discussion The Nordsieck-Gear (k,l) method for the numerical solution of the scalar O.D.E., y'(x) = f(x,y(x)), y(0) = y Q , for fixed stepsize, h, can be written as „(0) n+1 AY (1) I is the kxk identity matrix, I = t^,^,...,^] , e 1 - [1,0, ...0] A, is lower kxk diagonal submatrix of A, , {( i f)}= [( k t 1 ),( k ^ 1 ),...,( k : 1 )] T . K is the kxk Jordan canonical form companion matrix to the polynomial X . P is the kxk permutation matrix of columns of the identity matrix written in reverse order V, , is the Vandermonde matrix with element k-1 v., = i j (i=0,l,2 k-l,j=0,l,2,...,k-l), _ij j e = [1,1 1] • 23-1 Theorem Maximum Polynomial Degree Nordsieck-Gear (k,p) Methods exist for all (k,p) . Sketch of Proof rP +1 , The idea behind the proof is to reduce the lower kxk submatrix of A, (I-£ e ) to Jordan T canonical form leaving e AS invariant and to show that this term is non-positive. Since A, is upper triangular the problem can be reduced by one in dimension by defining the sequence. = 0, S . = (l"-ee^)[AS -{(^ 1 )}j. n+1 1 n 1 (A) Let W = L, Q, A, S . Here Q is a unique upper n k k k n n — T triangular matrix with first row equal to e.. . It — — — T takes A^(I-2.e ) into the form _ TT _1 V-L1 _T \ = \\ (I - £e i )Q k =I+K - {( i )}e i • (5) \ = (WVp^^V^ and that for the (k,p) case {( k+i )} = ( k+i )diag [( p )>( p+i )) __ >)C) . j k -1 -1 P k-p This reduces the general (k,p; case to the (k-p,l) case. References (1) Danchick, R. , "On the Non-Equivalence of Maximum Polynomial Nordsieck-Gear and Classical Method," Procedings of the Confer- ence on the Numerical Solution of Ordinary Differential Equations 17,20 October 1972, University of Texas at Austin, Lecture Notes in Mathematics, Springer, Vol. 362, No. VIII, pp. 92-106. B is similar under P to the companion matrix of k the polynomial X . L is a unique lower -T triangular matrix, whose first row is e. , that takes B into K according to K ■ Wk 1 (6) Under L, Q A, the S sequence goes into Vi = K[ v L kV (k t 1)} J • (7) Moreover e^W = ejx, Q, A, S = e^Q, A, S = e.A, S =e^A, S .(8) In 1 k x k k n l v k kn 1T( n lien The theorem is proved if it can be shown A sufficient that efw = < for n = 0,1,2, In - condition for this sequence of inequalities to hold, given (4), is that L kV (1 ? )} i° • < 9 > This is accomplished by showing that L kV (k l 1)} = (k+1)L kVk lg (10) and, by induction on k, that w; 1 - ° ■ (11) passage to the k > p > 1 case is made by showing that the general form of the Q, associated with a (k,p) method is 23-2 TWO DEVICES FOR IMPROVING THE EFFICIENCY OF STIFF ODE SOLVERS Peter Alfeld Department of Mathematics University of Utah Salt Lake City, Utah ABSTRACT The nonlinear systems of equations that a- rise when a stiff system of ordinary differential equations is tackled by an implicit linear multi- step method is normally solved using a modified version of Newton's method in which the Jacobian is kept constant for as many steps as possible. An initial approximation is obtained by applying a predictor, typically an explicit linear multi- step method. In this talk, two approaches to improving the efficiency of the above procedure are described. Firstly, it is suggested to use quasi-Newton methods that update the Jacobian approximation at each iteration-step for the solution of the nonlinear system. It is argued that the most popular such method, Broyden's first method, is unsuitable for the special nonlinear systems under consideration, and that instead Broyden's second method should be used. That argument depends very critically on the different roles that are played by the eigenvalues (of the Jacobian) with large and small, respectively, absolute values. Computational aspects are con- sidered, and numerical examples are given. Secondly, it is suggested to improve the initial approximation by interpolating correc- tions that were applied to previously obtained initial approximations. This device is motivated by considering the test equation y' = X . Nu- merical examples suggest that it works very well. 1. INTRODUCTION When an implicit linear multistep method is applied to a stiff initial value problem y' = f(y) ; y(a) = n e R m ; x 6 [ a ,b] (2) F n ( W = y n + k' hB n f( Vk'Vk ) k-1 + I (a.y -hB.f ..) = jtg J n+J J n+j' (3) has to be solved at each step of the integration (1). 2. THE USE OF QUASI-NEWTON METHODS In this talk, it is proposed to solve the nonlinear system (3) by a rank-one quasi-Newton method defined by: 1) y n+k given by some prediction operator. Let r = . 2) Solve B^slfj = - F (yf?J) where B [r J n+k n+k n y n+k n+k frl is an approximation to F'(y ,,) n n+k 3) Let y [ lt 1] = y lr l + s [ *} 'n+k 'n+k n+k u [r] = ( Ir+llj.p ( [r] n+k n^ y n+k ; n^ y n+k ; (4) [r+1] m [r] ^n+k n+k S n+k MV n+k ; n+k n+k " [r] T [r] (S n+k j V n+k [r] (where v ,, is specified below) n+k 4) Check for convergence. If not satisfied let r = r+1 and go back to step 2). then a nonlinear system of equations 24-1 A detailed discussion of quasi-Newton methods may be found in [1], [3], [A], [6]. Among others, the following properties are important in the present context: 1) Quasi-Newton methods converge superlinearly. 2) The update of the Jacobian can be replaced by an update of the inverse Jacobian. Normally, a version of Newton's method is used for the numerical solution of (3). The Jacobian is kept constant for as many steps as possible, but when it has to be reevaluated, or reestimated, this requires additional function evaluations. In a quasi-Newton method, this expensive procedure is replaced by a cheap up- date of the Jacobian at each iteration for the solution of the non-linear system (3)o This is more efficient as it requires no additional function evaluations. Broyden [3] suggests the following choices of the update (4) : v [r] = g [r] ( Broyden ' s first method) (5) n+k n+k v I r l = (b , ) T u ,, (Broyden's second method) n+k n+k n+k (6) Broyden's first method is generally consid- ered the better general purpose method (see [3], [5]). However, in [1] it is argued that for the special nonlinear systems under consideration Broyden's second method is more suitable. For lack of space the argument can only be outlined here. (A full description may be found in [1].) It is based on considering the different roles that are played by the eigen- values of the Jacobian F' with large and small, n respectively, absolute values. The argument can be summarized in three points: 1) A prediction y ,, can be obtained that is r ; n+k particularly accurate in components corre- sponding to eigenvalues of large absolute value. 2) Any quasi-Newton method reduces in a cer- tain sense to the secant method applied to each component (corresponding to one of the eigenvalues) of the solution of (3). For the secant method, the error at the r-th stage is roughly proportional to the previous two errors and inversely propor- tional to the first derivative of the function under investigation. In the pre- sent context, this means that the error in the component corresponding to an eigenvalue , say, is roughly proportional to the two previous errors in that component and Inversely proportional to A, . The above points show that the eigenvalues of F' with large absolute value have little n importance for the convergence of the iteration for the solution of (3). (They do matter for stability). The key point is the following: 3) Broyden's second method approximates partic- ularly well the eigenvalues of small absolute value. More precisely: Broyden's first method minimizes (over all rank-one updates) a sharp bound on the error of the eigenvalues, whereas Broyden's second method minimizes a similar bound on the reciprocals of the eigenvalues. Since these bounds are independent of the eigenvalues, Broyden's first method yields small relative errors of the large eigenvalues, and Broyden's second method yields small relative errors of the small eigenvalues (large and small in absolute value) . In [1], numerical examples (taken from [7]) are given that confirm that quasi-Newton methods are strong competitors of Newton's method, and that indeed the update (6) performs better than the update (5) (in the present special context). Numerical aspects considered in [1] include the following: [r] 1) Rather than updating B in (4), we can n+k fr] -1 update (Bl, £) > thus saving an LU- f actorization at each step. 2) When the step-size or order of (1) is changed, then it is a viable strategy not to change the Jacobian approximation at all. The basis for this device is the fact that the smaller the absolute value of an eigen- value, the less it is affected by such a change. 3. THE USE OF INTERPOLATION In [2], the following modification of the normal prediction process is described: y given by an explicit linear (7(i)) multistep method 'n+k y n+k + n+k (7(ii)) y ,, obtained by applying a variant (7(in)) n +k c ., , , . . of Newton s method or a quasi- Newton method to (3). d ,, is obtained by interpolating n+k d n+J " y n+j " y n+j (J < k) • This approach depends on d ,. being a smooth function of n+J x . Using an asymptotic expansion of the n+j 24-2 .runout ion error tor investigating this Ther. ie test equation y' ■ Vy (8) nsidered. It turns out that the corrections d , show the same qualitative behaviour as n+k the numerical solution of (8) by (1). It follows that, for example, Gear's back- ward differentiation formulas (see [8]) are suitable for the approach (7). On the other hand, the Trapezoidal Rule leads to oscillating global errors, and thus to oscillating corrections d . . . However, these oscillations occur only for n+j eigenvalues of large absolute value. Since we have seen in section 2. that only eigenvalues of small absolute value matter for convergence, it follows that the device (7) can also be applied to the trapezoidal rule. Moreover, particular predictors can be designed that dampen the oscil- lations in the corrections. In [2] numerical examples are given that con- firm the above statements, and the viability of the approach (7) . REFERENCES [1] ALFELD, P., Quasi-Newton Methods for Stiff Ordinary Initial Value Problems, submitted for publication. [2] ALFELD, P., A Device for Improving the 'iciency of Stiff O.D.E. Solvers, in preparation. [3] BROYDEN, C. G., A class of methods for solving nonlinear simultaneous equations, Math. Comp. 19 (1965), pp. 577-593. (4) BROYDEN, C. G. , Quasi-Newt on Methods, in W. Murray (ed.), "Numerical Methods for Unconstrained Optimization", Academic Press, London and New York, 1972, pp. 87- 106. [5] DENNIS, J. E., A brief introduction to qua8i-Newton methods, in Proceedings of Symposia in Applied Mathematics, V. 22, American Mathematical Society, Providence, Rhode Island, 1978, pp. 19-52. (6] DENNIS, J. E., and J. J. MORE, Quasi- Nevton Methods, Motivation and Theory, SIAM Review 19 (1977), pp. 46-89. [7] ENRIGHT, W. H. , T. E. HULL, and B. LIND- BERG, Comparing numerical methods for stiff systems of ODEs, BIT 19 (1976), pp. 10-48. [8] GEAR, G. W. , Numerical Initial Value Problems in Ordinary Differential Equa- tions, Prentice-Hall, Englewood Cliffs, NJ, 1971. 24-3 STRUCTURAL STABILITY AND THE NUMERICAL METHODS FOR STIFF DIFFERENTIAL SYSTEMS* Shohei FUJITA Department of Computer Science, Faculty of Engineering Tokyo Institute of Technology O-okayama 2-12-1, Meguro-ku, Tokyo 152, JAPAN Abstract — Several computational methods based on the backward differentiation formulas have been proposed for the numerical solution of stiff systems of ordinary differential equations. These methods, however, suffer from serious difficulty in that they are unable to cope with the stiff systems in the presence of the eigenvalues of the Jacobian which are close to the imaginary axis. This paper expli- cates this phenomenon from a point of view of struc- tural stability. A necessary and sufficient condi- tion is presented for structural stability for n- dimensional system of linear ordinary differential equations which is used as a test problem. This proves that the region of stiff stability contains structurally unstable systems. The poor performance of the algorithms when some of the eigenvalues of the Jacobian are close to the imaginary axis is caused by the presence of structurally unstable systems. Furthermore, a reliability of testing for stiff systems is discussed from a point of view of topological equivalence. A practical problem arising from dynamic security assessment for the large-scale interconnected power system is also discussed. 1. Introduction The purpose of this paper is to demonstrate importance and usefulness of the qualitative theory of differential equations, i.e., the theory of structural stability, in the design of numerical algorithms for stiff differential systems. 2. Structural Stability for n-dim. LODEs lity for n-dim. LODEs and present a clear proof of this result based on the topological equivalence of linear flows [7]. We shall begin by explaining a general concept of "structural stability". A flow $ on a topological space M is a continuous map: (J): R x M such that (i) (0, x) (ii) <|>(t M "!' (t 2 , x)) = >(t + t 2 , x) for x £ M and t , t„ £ R, where R is the real line. A trajectory of 4> through x £ M is a canonically ordered set {(J)(t, x) | t £ R} . Definition 1: Two flows 4> and i)j on M are said to be topologically equivalent if there exists a homeomorphism of M onto itself carrying each trajec- tory of (f) to a trajectory of \p together with preser- ving orientation. For linear flows Definition 1 is shown to be equi- valent to the following definition: Definition 1' : Two flows and i/j on M are said to be topologically c-equivalent for c > if there exists a homeomorphism: M ■+ M such that h • , where | | • | | is some norm. Definition 3: A vector field f £ V(M) is structurally stable if there exists a nbhd N(f) of f in V(M) such that for every g £ N(f), the flows of f and g are topologically equivalent. Now we can show a new result pertaining to n-dim. LODEs. Let us consider a family of n-dim. LODEs: x = Fx, x £ R°, F c L(R n ) (1) which is used as a test problem, where L(R ) denotes 25-1 n mat t i ha follow- ing theorem: Theorem 1: The n-dim. LODEs (P is structural- table if and only if Re \ ± (7) t 0, I n, Re ^.(") denotes the real part of the eigen- values of the matrix. Remark 1: It seems that there is no descrip- tion of a necessary and sufficient condition for the n-dim. LODEs (1) to be structurally stable. The structural stability for n-dim. ODEs was inves- tigated in [11] which however depends on the "c - homeoraorphism" . E xample 1: Consider the 4-dim. LODEs: x = AAA _1 x, x E R 4 , t ( [0, 20/r,] (2) A - diagli, -i, -ct+ Lui, -a - Im], I > -1 1 i ■1 1 -i i -i -i -i 1 -1 x(0) = The eigenvalues X l,2 " t 1 ' 3,4 -a + im Hence, by Theorem 1, the LODEs (2) is structurally unstable. We can find, Jn fact, tests examples where the BDF of order A to 6 for the 4-dim. LODEs (2) come to unstable [14], Proof of Theorem 1: First, we note that the vector field of the n-dim. LODEs (l) p |s f(x) = Fx, the flow is specified by 41 ( C , x) = e x, and V(M) is the set of all n x n matrices, which we denote by L(R n ). The nbhd of F e L(R n ) is any subset N(F^ in L(R n ) which contains a set of the form:{G c L(R ) 5) for some 5 > and | | • | j uniform Lemma 1 [7]: Let ■i and fy be flows on R gene- rated by x = Fx, x e R , and y = Gy, y e R , respec- tively. Then, the two flows 41 and ty are topologica- ls equivalent if and only if In(F) = In(G) and G* » cF*, where In(-) is the inertia of the matrix, c is a positive constant, and F* denotes the real Jordan form of the matrix F corresponding to the eigenvalue with zero real part. We are now in a position to prove Theorem 1: (Necessity) Assume that the n-dim. LODEs (1) is structurally stable but Re > (F) = for some i. Chat By Definition 3 it follows that there exists a nbhd N(F) of Ft L(R n ) such that for every G E N(F) the flows of F and G are topologically equivalent. By choosing any nbhd N(F) of F, however, there is an G e N(F) which does not satisfy the condition in Lemma 1. This contradicts to Lemma 1. Thus we conclude that Re A (F) t 0, 1 - 1, ... ,n, for the structurally stable n-dim. LODEs (1). (Sufficiency) Suppose Re A . (F) 4 0, i = 1, . . . , n. Then there exists a nbhd N(F) in L(R ) such that for any G E N(F), Re A (G) i 0, i « 1, . . . , n, and In(F) = In(G). Hence, by Lemma 1 it follows that the flows of F and G are topologically equi- valent. 3. Stiff Stability Contains Structurally Unstable Systems The most popular methods for solving stiff diffe- rential systems are based on the backward differen- tiation formulas. The crucial drawback of these methods is the poor stability property when the Jacobian matrix has eigenvalues close to the imagi- nary axis (e.g., [1], [2], [9], [13]). A reason of this difficulty can be interpreted as follows. First, in view of Theorem 1 and the definition of stiff stability ([5], [8]), we have the following theorem: Theorem 2: The region of stiff stability contains structurally unstable systems. To summarize: the implicit BDF of order 1 to 6 are stiffly stable, so that by Theorem 2, structu- ral instability is possible for the BDF of order 1 to 6. 4. Testing — Topological Equivalence An ultimate aim in discussing stability of numeri- cal methods is to determine the stability of the methods for nonlinear systems: x = f(x) , x c R , x(0) (3) The property we would like to demand from a method is that, for a given step-length and order, the global error remains as small as possible. In order to forecast the behavior of a method on a non- linear system, it is customary to consider the test equation: .1 x = Ax, x e R (4) for all A in some set including the set of eigen- values of the Jacobian matrix J = 3f/3x. This test problem is, however, too simple to give a reliable indication of the performance for stiff systems unless some additional conditions are satisfied. Little has been rigorously proved concerning error propagation in the numerical solution of nonlinear stiff systems [1]. There does not exist a general theory of stability of numerical methods such that the stability behavior of the method when applied to nonlinear stiff systems can be determined in a systematic manner. What are conditions under which the (local) behavior of the nonlinear stiff problem (3) is determined by the solution of the linearized system? This question is partially answered by the topological equivalence Lemma: Lemma 2 [6]: The local behavior of the NODE (3) near the singular points is topologically equi- valent to the behavior of the linearization of (3) if and only if ReA (J) + 0, i - 1, ... , n. Mo the reover, let tt , tt , and tt denotes the number of e eigenvalues of the matrix J with positive, zero, and negative real parts, respectively. The topolo- gical type of singular points is determined by the intergers tt , tt , and by the behavior of the phase trajectories on some invariant submanifold whose dimension is equal to tt . When the linearized system has an eigenvalue with real part zero, know- ing about the behavior of the linearization of (3) at singular points gives no information about the behavior of (3) itself near the singular points. 25-2 Example 2: Let us consider a practical problem arising from dynamic security assessment for the large-scale interconnected power system. A mathema- tical model for the power system is described by 6. 1 - ETC . - D.6. E A. .cos(6 . J = 1, ij (5) One of the problems of dynamic security assessment is to estimate critical clearing times. Hence, we must know the behavior of the solution of the second-order NODEs (5) which is possibly unstable due to contingency. If the effects of damping torque is approximated by "uniform damping", i.e., D./M. 1 x C, 1*1, then the NODEs (5) result in n, f(x) - Ci, x e R n-1 (6) The local behavior of the NODEs (6) near singular points is topologically equivalent to the behavior of the linearization of (6) if and only if .e A. l (K ) * 0, ' I " = • J o " CI - 1, 2n-2, On the other hand, when all eigenvalues of J_ are real, to each negative eigenvalue of J. there corresponds to either two negative real eigenvalues of K_ or to a pair of complex conjugate eigenvalue with real part equal to -C/2. Thus, in the case of without damping, i.e. Re > i (K Q ) = 0, it is possible that 0. Hence, we need a full nonlinear system to obtain information about the behavior of (6) itself near the singular points. 5. Concluding Remarks The qualitative theory of differential equations, i.e., the theory of structural stability is an important factor which must not be over looked in designing algorithms for stiff differential systems. It may be hoped that the theory of structural stability makes it possible to give a theoretical foundation to the idea of P. E. Zadu- naisky in estimating global errors; and may be also useful to study the dangerous property ([1], [10]) of numerical methods for stiff differential systems with eigenvalues which change from large negative values to large positive values during computation. [2] Enright, W. H., Hull, T. E. and Lindberg, B. Comparing numerical methods for stiff systems of O.D.E.s. BIT, vol.15, pp.10 - 48, 1975. [3] Fujita, S. Structural stability and genericity for n-dim. dynamical systems, with applica- tions to system theory, in Research Reports for Coordinate Research Program (A) (suppor- ted by the Ministry of Education, Japan, under Grant 239002), pp.117 - 144, 1978. [4] Fujita, S., Tsuji, K. and Fukao, T. Topolo- gical equivalence of state spaces in the interconnected power system in emergency control. IEEE Trans . Automat . Contr . , vol. AC-24, no. 3, 1979 (to appear). [5] Gear, C. W. Numerical Initial Value Problems in Ordinary Differential Equations . Prentice- Hall, 1971. [6] [7] [8] [9] [10] [11] [12] [13] [14] Hartman, P. Ordinary Differntial Equations . John Wiley & Sons, 1964 & 1973. Ladis, N. N. The topological equivalence of linear flows. Diff . Urav . , vol.9, pp.1222 - 1235, 1973. Lambert, J. D. Computational Methods in Ordi - nary Differential Equations . John Wiley & Sons, 1973. . The initial value problem for ordinary differential equations, in The State of the Art in Numerical Analysis (Jacobs, D. A. H. ed.), Academic Press, 1977. Lindberg, B. On a dangerous property of methods for stiff differential equations. BIT , vol.14, pp.430 - 436, 1974. Markus, L. Structurally stable differential systems. Ann . Math . , vol.73, pp.1 - 19, 1961. Peixoto, M. M. (ed.). Academic Press, 1973. Dynamical Systems , Skeel, R. D. and Kong, A. K. Blended linear multistep methods. ACM Trans . Math. Soft . vol.3, no. 4, pp.326 - 345, 1977. Skelboe, S. The control of order and step- length for backward differentiation methods. BIT , vol.17, pp.91 - 107, 1977. References* [1] Dahlquist, G. Probelms related to the numerical treatment of stiff differential equations. International Computing Symposium 1973, (Gun- ther, A. et al. eds.), North-Holland Publ. Co., * An extensive bibliography up to 1973 can be found In Stiff Differential Systems (Willoughby, R. A. ed.), Plenum Press, 1974. 25-3 On the Theory and the Construction of Cycl Le Mult is i op Peter Albrecht Nuclear Research Center Julich neral theory of high order cyclic methods is given in /l/ and /2/; in this contribution we want to present some practical consequences. The theory is based upon the result that the functional ("A- '') A-. - °snp • si J K ■ 1-1 \i o. £ R s 1 J (1) ty which depends upon a real (s.s)-matrix S with |A J |i?:const . , jflH, is (the most adaquate) stabili functional for a very wide class of discretization methods, namely the A-dependent methods ("A-methods") defined by (h); o o z.=Az. +h?(x. ,,z. .,z.;h), j=l(l)p (2a) where ; (h),z.eR S , h£(o,h ] ( MJ-2M+1 MJ-2M+2 +h Mj-M / K oMj-M 4 B 4 f • M Mj ? 2 If . +...+S f . U Mj-M P l Mj ^ p Mj-M 3 M f M Mj L j-i with j=l(l)p, z =(o, ,o,n ) dR . +h <2>(x. „ ,z. „ ,z. ;h) M o ' ' ' o These one-step methods have been previously studied by Watts and Shampine /9/ and by Bickart, Burgess and Sloate IV, Ibl. Let us assume, that all stages have order q except the last one which is supposed to have order (q+1); then t=(c • ... ;c .,;o) (c . : error const. q+1 q+1 q+1 of the r stage) and the method converges with order (q+1) because (6a) is satisfied. Another application of the order condition (6a) is the following: If we construct a k-cyclic k-step method with (unstable) stages of order (2k-l), determining the k free parameters such that the method is stable and satisfies (6a), then we obtain Donelson and Hansen's /4/ k-cyclic k-step methods of maximum order 2k. Thus, our theorem provides a rather simple and straightforward approach to them, discarding the concept of "auxiliary systems" introduced in I hi. As the stability regions of these methods are fairly small for k>3, we may ask, whether methods with slightly reduced order possibly have better stability properties. This is the case and, in fact, one can even achieve stif f _stability, maintaining unusual high orders. In particular, we have found stiffly stable 2-cyclic 2-step methods of order 3 2-cyclic 3-step methods of order 4 . . 3-cyclic 3-step methods of order 5 4-cyclic 4-step methods of order 6. The i h methods La rcm.-irkable thei ible, Id near pnumber -wn. We now consider M-cyclic k-st:p methods of the iorm z =C(h); lz.=Uz. +h(B I F.+B I ,F. J , j=l(l)p o L j U j-i k-V r s with F : = (o,...,o,f , . . . , f , , ) dR ; s : =max(M,k; : o o •'-i j Mj+k-s Mj+k-1 and (s ,s)-matrizes L,U,B ,B , where B and B have r a similar structure as L and U and contain the B, . Tendler, Bickart and Picel /8/ have constructed stiffly stable methods of this type with_order k for the special case B J . = C / . In the general case one obtains the maximum order methods mentioned in (7) they were constructed in the following way: Let be p.(i=l(l)s) the zeros of P (u) : =det (pL-U) and (j. (i=l(l)s) those of P(y) : =lim H~ M Q(y;H) , 1 H-«° where H = Ah and Q is the characteristic polynomial Q(p;H):= det (u (L-HB L )- (U+HBy) ) . The order q of the stages is chosen as high as possible and such that a sufficient number of free parameters is available to satisfy the following conditions : 1. |u.|< 1, i=2(l)s ("stability at H=o") (8a) 2. |y.|< 1, i = l(l)s ("stability at H=°°") (8b) and, if possible, 3. rank (L-U,c) = s— 1 (order condition) (8c) ,1 M T r c:=(c .,...,c ,) ; c , :error const, of q+1' q+1 q+1 th r stage (8c) is obtained from (6a) by multiplication with L using Lt=c. Each method that satisfies (8a-c) is then testet on stiff stability performing the mapping of the unit circle into the H-plane by Q(y;H)=o. Up to now, a general existence theory of stiffly stable, cyclic methods with maximum order doesn't exist. Also the question of practical performance and implementation have not yet been sufficiently investigated. A series of tests, run at the Nuclear Research Center in Julich, Germany, furnished promising results. From consistency we have always \i„ *1 . 26-2 References Ubrecht, P.: Explicit, optimal stability funccionals and their application Co cyclic discretization methods. Computing 19 (1978). '. brecht, P.: Die numerische Behandlung gewohn- licher Dif f erent ialgleichungon , Hanser, Miinchen 1979. /3/ Bickart, T.A.; Burgess, D.A. ; Sloate, H.M.: High order A-stable composite multistep methods for numerical integration of stiff differential equations. Proc. 9 Annual Allerton Conf. on Circuit and System Theory, Univ. of Illinois, 465-473 (1971). Ikl Donelson, J.; Hansen, E. : Cyclic composite multistep predictor-corrector-methods. SIAM J. Numer. Anal. 8, 137-157 (1971). /5/ Skeel, R.: Analysis of f ixed-stepsize methods. SIAM J. Numer. Anal. 13, 271-3oo (1976). Ibl Sloate, H.M. ; Bickart, T.A.: A-stable composite multistep methods. J. Assoc. Comp. Mach . 2o, 7-26 (1973). Ill Stetter, H.J.: Analysis of discretization methods for ordinary differential equations. Berlin, Springer 1973. /8/ Tendler, J.M. ; Bickart, T.A. ; Picel.Z.: A stiffly stable integration process using cyclic composite methods. Assoc. Comp. Mach. TOMS Vol. 4 No. 4, 339-368 (1978). /9/ Watts, H.A.; Shampine, L.F.: A-stable block- implicite one-step methods. BIT 12, 252-266, (1972). 2&-3 An automatic differential equations solver based on linear multistep methods , by ZAHARI ZLATEV s Royal Veterinary and Agricultural University, Institute of Mathematics and Statistics Thorvaldsensvej UO , DK - 1871 Copenhagen V, Denmark and PER GROVE THOMSEN inical University of Denmark, Institute for Numerical Analysis, DK - 2800 Lyngby, Denmark :der the non-stiff system of ordinary dif- ferential equations (1) y'-f(x.y); y(a)=n; x£[a,b]j f£R s ; s>l. . r.e the grid a=x r| and the computational work per step are smaller than those of the other linear multistep methods. However, if the absolute stability requi- reme Lag the numerical integration of the particular problem are dominant over the accuracy see Shampine [5] ) then the possi- bility of using large stepsizes will be restricted, because these methods have small absolute stabili- ty regions. An attempt to improve their performan- ce in such situations has been done in [4] and [6] where P.EC, ..E schemes (which have better stability properties, see the plots in [6] ) are adopted instead of P EC E schemes. Some linear multistep formulae with large absolute stability regions have been derived by Thomsen and Zlatev [7]. A code based on these formulae is discussed in this paper. * Consider the following scheme : 'VV - p k EC k + i E k-1 (2) CAVH +(1 AW h .VikVi ! i=0 (3) Ck^wC^ k-1 (U)y n + k = Vn + k-1 + (1 -\ )y n + k-2 +h . J: n e ik f n + i 1=0 +h6 kkW (5) W f 1 !6 - a) g n+k,0 = W Ck,o"Ck J {6 ' b) g n + k,i =(g n + k,i-rVk-1,i-1 )/U n + k"Vk-i )- ' C?,a) r 0,T X n+k Vk-1' r o,1 = Vk x n+k-2 ; (7 - b) r o,r r o,i r o,o-1 ( ^ 1)/ J' 0=2(1)^*1^ (7 - c) fo.j-'o.l'oj-i"- 1 ^' Ml*** (7 - d) r iJ- 2 +(1 - a k )f i,2 i *"«>**■ Then the ( a k .\) _ VSVFM ie based on the following formulae: (9) y r**, = i * y n + k-1 + (1 - a k ) Vk-2 + . Z n Y i,k g n + k,i i 1=0 10) y =i / „ + M-ci )v ♦ I y . & ; 'n+k k*n+k-1 V y n+k-2 . i Y i,k fe n+k,i * * E n+k,0' Fig. 1 (the differences g ,, . in the sum in (10) can be found by the use of (6.b) with g replaced * ty g )• In our implementation many ideas from the most efficient codes ( [1], [3], [4], [6] ) were used. The code can be used with different va- * * lues of a and a , . By the choice a =a =1 , k=2(1)l2, an Adams VSVFM may be obtained. This allows us to make comparisons in which only the use of different formulae in the VSVFM, but not the stepsize and order selection strategy, will have influence on the results. Many problems from [3], [4] and [6] were run. As an example consider problem 12 in [4] with B^Uv, B 2 =3v, B,=-v, B^O.OOlv, v=1(1)12, e-10 , [a,b]=[0,50]. The numbers of function evaluations (ND) are given in Table 2. The results of the comparison between the * Adams VSVFM 's and the (a, ,a, ) VSVFM 's can be k* k' summarized as follows: (i) The (a ,a ) VSVFM 's require more storage and computational work per step. Moreover the error constants of (2) and (k) are larger than those of the corresponding Adams methods. Therefore it is not efficient to use these methods if the accuracy requirements are dominant. But note that very often the number of function eva- luations was not larger ( a careful comparison 27-2 V Adams VSVFM (a k ,a*) VSVFM 1 289 227 2 1*59 3U9 3 629 un I* 801 597 5 871 719 6 1139 839 7 1309 963 8 mi 1085 9 1577 1210 10 1815 1332 11 198U 1U52 12 2155 1 Table 2 showed that this is so because higher orders were selected by the (a^,o ) methods ). (ii) If the absolute stability requirements are dominant then the (a ,a ) methods will nor- mally use a smaller number of steps in the inte- gration. (iii) The zero-stability properties of the a, ,ol ) methods are the same as those of the k k Adams methods, see Zlatev [8]. (iv) A one step predictor-corrector scheme ( the Euler formula as a predictor and the tra- pezoidal rule as a corrector ) is used in the starting stage and in the situations where the order should be dropped to one (as an example consider the case where discontinuities in the right -hand side of (1) appear for some values of x in [a,b]), see more details in [6]. Our main conclusion is that ( a v» a v) methods :ar. s^::essfully be used instead of the Adams ods in the case where only the absolute sta- - S requirements restrict the selection of large stepsizes. REFERENCES 1. Gear, C. W. : Algorithm 1*07, DIFSUB for so- lution of ordinary differential equations. Comm. ACM '\U , 185-190(1970) 2. Gear, C. W. , Watanabe, W. S. : Stability and convergence of variable order multistep methods. SIAM J. Numer. Anal. 11, 101*1*-1058(197^) 3. Hull, T. E., Enright, W. H., Fellen, B. M., Sedgwick, A. E. : Comparing methods for or- dinary differential equations. SIAM J. Numer. Anal 9, 603-637(1972) 1*. Krogh , F . T . : On testing a subroutine for numerical integration of ordinary differential equations. J. Assoc. Comput . Mach. 20, 5U5-562(1973) 5. Shampine, L. F. : Stiffness and non-stiff differential equations solvers. In: Nume- rische Behandlung von Deferentialgleichungen (L. Collatz, ed.). Int. Series Numer. Math., Vol. 27, pp. 287-301. Basel, Switzerland: Birk- hauser 1975 6. Shampine, L. F., Gordon, M. K. : Computer so- lution of ordinary differential equations: The initial value problem. San Francisco, Calif. : Freeman 1975 7. Thomsen, P. G., Zlatev, Z. : Two-parameter families of predictor-corrector methods for the solution of ordinary differential equa- tions. Institute for Numerical Analysis, Technical University of Denmark, Lyngby, Denmark. Report No. NI-77-08, 1977 8. Zlatev, Z. : Stability properties of variable stepsize variable formula methods. Numer. Math. 31, 175-182(1978) 27-3 SOFTWARE DEVELOP. IIT FOR STABILITY ANALYSIS OF NONLINEAR DIFFERENTIAL SYSTEMS R. Leonard Brown Dept. of Applied Mathematics and Computer Science University of Virginia Charlottesville, VA 22903 1. PROBLEM STATEMENT - Since the exact solution of a system of m first order ordinary differential equations y* = f(y,t), y(t Q ) = (1) can be numerically approximated at only discrete points, the comparison of analytic and numerical solutions to determine accuracy and stability can only be made at discrete values t 0, ,n. Since most numerical methods attempt to keep the stepsize h. = t. „-t. constant for as many steps as l i+1 l practical subject to accuracy and stability constraints, the study of local accuracy and stability can be carried out for constant h. In particular, define the analytic k-step solution sequence Y( t - wv k+1 .>, y< t - 1 f 2 >'--'* ( V ) and the numerical k-step solution sequence to be Y n = (y n-k +1 ' Vk + 2' — ' V Where y i i3 , the numerical approximation of the solution y(t ) of (1) and each of y(t ), y is an m vector. This paper will describe a software package of interconnected programs and subroutines which can be used to compare the stability properties of Y generated by a particular multistep method with the properties of Y(t ) for a particular nonlinear differential function f(y,t). Fig. 1 shows a flow diagram of this package. The usual analysis of the stability of ( 1 ) involves two concepts: the equilibrium solution and the domain of attraction. Define the equilibrium solution z(t) to be the analytic solution of ( 1 ) where z = z(t ) is chosen such that z'(t Q ) = f(z Q , t Q ) = 0. (2) For autonomous equations, the solution is z(t) = z . The domain of attraction can be found if there exists a function v(y) > for y i 0, v(0) = 0, such that dv(y)/dt i in a doaain D about the equilibrium z(t). This region Is the domain of attraction and, as t increases, the difference between any two solutions of (1) generated by different initial values w and x within D decreases with respect to v, i.e. Iv(w(t+h))-v(x(t+h))| <|v(w(t))-v(x(t))|. If the function v(y) = y Qy is used, where Q is a positive definite matrix chosen by known techniques [1], then v(y) is for the scalar product u fc Qv. In fig. 1 , box 2 represents a routine which solves for the m(m+1)/2 unique values of Q by dividing f(y,t) into its linear and nonlinear parts so that f(y,t) = By + g(y,t) Then a symmetric (or Hermitian) Q is found such that B fc Q + QB = I where I is the m by m identity matrix. While this routine has not yet been written, the linearization can be automated by finding the approximate Jacobian B at t = using column by column finite differencing, requiring m evaluations of f(y,t). Box 3 of fig. 1 represents a routine which, given an initial guess z inside D will find either z , the equilibrium at t = 0, or else z •, the equilibrium solution in a plane specified by holding m-2 values constant, if such a restriction does not produce algebraic contradictions. This routine is also being written and is discussed in section H below. In the sequel, z will represent either z or z * as appropriate. 2. COMPUTATION OF STABILITY REGIONS - A concept similar to the analytic domain of attraction is sought for the analytic and numeric solution sequences. Dahlquist [2] has devised a function V G,K,n ( V ' ill |l g ij for some positive definite matrix G = (g,^, where is the scalar product above. LinigSr and Odeh [3] have considered how to choose G to obtain the nonlinear equivalent of A-stabllity, and Dahlquist [1] mentions in a recent report that a program to compute G has been written at Stanford. This is represented by box 1 of fig. 1 . To facilitate computational and analytic manipulations with respect to multistep methods, it has been shown [2,5] that the numerical solution 28-1 usinj k-step multistcp formulas = * ay ,+hb.f(y t . ) r=0 i'n-i 1 'n-i n-i stability properties similar to numerical solutions using one leg k-step formulas k has power series x(t) = jj. a^ , y(t) b.t defined recursively by i=0 i Vn-i + hf( i=0 V.-i'i=0 Vn-i 5 if .: ', b. = 1. Software is being developed which, given a point within the resulting domain finds the boundary of the region D 1 where V Q k n (Y.,) - V k h ^V - °' which is the discre te analog of the domain of attraction -here v'Cy.) <. 0. This can be used to compute D 1 for either the analytic or the numerical solution sequence, provided the analytic sequences Y(t ) and Y(t ) are available. The k-step numerical solution sequence is also computed using y and k-1 exact (analytic) previous points to compute y using the one leg k-step method under consideration. Box 1 of fig. 1 represents a compiler module which transforms a specification of the differential function f(y,t) into two FORTRAN subroutines S and F. S(T,Y0) returns the correct solution vector at t = T given that Y0 is the initial value at t = 0, and using techniques described in section 3 below. F(T,Y) returns the derivative of Y at time t = T given that y(t) = Y. Box 5 of fig. 1 represents working software which uses brute force search and then bisection on up to 80 lines radiating from z along a 2-dimensional plane, storing the result in an array STAR. The stability region D 1 ' such that V h^ Y n^ <. V f can also be computed. V» is max V * f Y) over all Y on the boundary of D 1 . This corresponds to all initial values that generate solutions that remain inside D 1 after k steps. Box 6 of fig. 1 represents the working software that computes and stores up to 80 points of the boundary of D' ' . These points are computed using bisection between z and the boundary of D' stored in array STAR. It ha3 been shown [5] that, if a domain of attraction D exists for v(y) = y Qy that is convex, then the boundaries of D' and of D'' can be computed using rays from some point within D' ' . 3. AUTOMATED ANALYTIC SOLUTIONS - Since the analytic solution is required for k-1 steps before y Q , means of obtaining y(t ), , y(t ) are needed. Translation software, box 1 of fig. 1, is being developed to manipulate the specification of a system of m non-linear first order ordinary differential equations to produce a FORTRAN subroutine which, given y at t = and some value t, outputs y(t) using power series solution techniques provided t is in the radius of convergence of the power series. For example, x' = -x-2y+x 2 sin(t) y 1 = 5x-y+y/(t+JO ( " a i-r 2b i-1 < V b = V i-1 ln( V jio. a j b i-i-j )/i - b i '- C5a 1-1 -(t +3)/(t +«)b^ 1 )/l. The regions D, D' , and D' ' for the analytic solution with q^ = 37 A4, q = 16/44, q = q ?1 = 3/1*14 , G = I, h = .1, and k = 2 appear in fig. 2. 1. EQUILIBRIUM - Although any 'point inside D' • can be used to generate the boundary of D' and D 1 ' , if the equilibrium (2) is used then good results for the analytic case are insured. Software corresponding to box 3 of fig. 1 is being written to handle computation of z , possibly using an n dimensional chord method starting with a value that is known to be inside D' ' . This is an application of the software described in secton 2 since a message that V(Y ) - V(Y Q ) > can be generated to indicate that the point being used as the origin of the rays used to find the region boundary is itself not in the domain of attraction. A special case occurs when, for an n-dimensional system, n-2 initial values are held constant and the equilibrium initial value for the other 2 variables is found. This is both computationally simpler and more practical since only 2-dimensional regions can be displayed on most graphics equipment. Since this routine is not available in the current software package, a z is picked using knowledge of the function f(y,t) being investigated. 5. GRAPHICAL OUTPUT - All of the currently " working software has been gathered as subroutines to an interactive control MAIN program which uses a simple command language (the first n positive integers) to perform the n available functions. For those subroutines not yet available, such as the computation of the equilibrium, a default (0) is provided and execution of the command will allow the user to input a usable value. The command produces a list of available commands; each command prompts the user in its use. The graphics package, using TEKTRONIX PL0T10 software connected to a Control Data CYBER 172, prints the left, right, lower, and upper bounds of the figure to be drawn, then allows the user to override these limits. This allows all figures to be drawn to the same scale for comparison. A hardcopy device can be used to save the figures. The example in section 6 was drawn using this package. 6. EXAMPLE - Consider the analysis of the complex equation y' = Ly where L = (In y Q )/h, which has two equilibria, a stable one at z = and an 28-2 unstable one distance from r 1, for v(y) = y y, the Euclidean I. The analytic solution and G sequence is y. i y , so the analytic domain of attraction and stability region are both bounded by the unit circle, and the unstable equilibrium is on the boundary. Fig. 3 shows D' about 1 for the numerical sequences generated by the Euler predictor (dotted), the Euler predictor with one backward Euler corrector (dash), and the 2-step Adams explicit formula (solid). Only implicit formulas that are stable at ™ have domains of attraction that include z = 0. In this case, there is no stability region D' ' since the circle of radius V* about y = is not inside D' . 7. ACKNOWLEDGEMENT - This work has been supported in part by grant NSG 1335 from NASA through its Langley Research Center in Hampton, VA. 3. Liniger, W. And Odeh, F. , On Liapunov Stability of Stiff Non-Linear Multi-step Difference Equations, AFOSR-TR-76-1023, IBM Thomas J. Watson Research Center (1976). 4. Dahlquist, G. , G-stability is Equivalent to A-stability, Report NA-7805, Dept. of Information Processing, Royal Institute of Technology, Stockholm (1978). 5. Brown, R.L., Stability of Sequences Generated by Nonlinear Differential Systems, to appear, Hath. Comp. (April, 1979). 6. Brown, R.L., Stability Analysis of Nonlinear Differential Sequences Generated Numerically, to appear, Int. Jnl. for Computers and Mathematics with Applications . REFERENCES 1. Lehnigk, S.H., Stability Theorems for Linear Motions . Prentice-Hall, Englewood Cliffs, N.J (1966). 2. Dahlquist, G., On Stability and Error Analysis for Stiff Nonlinear Problems, Report NA-7508, Dept. of Information Processing, Royal Institute of Technology, Stockholm (1975). Fig. 2 - Nonlinear time dependent example 0.1 T rig. i J - Linear example; D' for three numerical methods. . u Fig. 1 - Flow graph of Software System. f(y,t) power series & function E S(T,Y0) LU F(T,Y) equilibrium _. P '■* *-. compute Q foi y'Qy compute boundary of n compute boundary of D' ' _0 t compute G foi V G,k,h (Y) - — ;UD 1 graphics display 28-3 SECOND DERIVATIVE EX- TENDED BACKWARD DIFF- ERENTIATION FORMULAE FOR THE NUMERICAL IN- TEGRATION OF STIFF SYSTEMS J . R . Cash Department of Mathema- tics, Imperial College, South Kensington, London S.W.7, England In two recent papers the present author has introduced a class of extended backward differentiation formulae suitable for the approximate numerical integration of stiff systems of ordinary differential equations of the form &= f(x.y). y(x o ) (l.D Extended backward differentiation formulae take the general form L 5, j-o y . - h(6. f . « n + j k n+k 5 k+l f n+k+1 (1.2) where the a.'s and B.'s are constants, J J y is the approximate numerical solu- n ♦ j tion obtained at x . and f ■ n + j n + j £ ( x .,y .). These formulae are simi- •*• n + j n + j lar in structure to the advanced multi- step methods of Williams and de Hoog and to a subclass of the composite multistep methods considered by Bickart and his co- workers. However, the way in which the required solution of (1.2) is generated is entirely different from that consid- ered by others. In an earlier paper it was shown that if these extended formulae are used in a certain precisely defined 'predictor-corrector' mode, with a con- ventional backward differentiation formu- la being used as a 'predictor' and (1.2) being used as the corrector, then it is possible to derive L-stable schemes of orders up to 4 and A(a)-stable schemes of orders up to 9. Numerical experi- ments with the resulting algorithm on a wide class of stiff test problems have indicated that it is often competitive with Hindmarsh's Gear Rev. 3 which is widely accepted as being one of the most efficient general purpose codes for the integration of stiff systems currently available. Numerical experience has also shown that when integrating stiff systems of equations with jacobians having large complex eigenvalues lying very close to the imaginary axis using extended backward differentiation formulae it is desirable that all of the basic integra- tion formulae should be highly stable (not just A -stable but A(a)-stable with a close to it/2). This is because an algo- rithm employing A(a)-stable formulae with a small sometimes selects a high order formula with a small value of h to pre- serve stability rather than correctly choosing a more stable lower order formula with a larger value of h. The requirement that all of our formulae should be very highly stable leads us naturally to con- sidering extended backward differentiation formulae employing second derivatives. These take the general form y n+k " y n+k-l k+1 h / B.f . + h 2 y, g , t- j n+j 'k 6 n+k j-0 (1.3) , _ df (x,y(x)) where g , = -r— 6 n+k dx n + k ,y=y n+k Note that we choose a. = 0, j < k-1 to ensure stability when h = 0. Following the general strategy adop- ted for the implementation of (1.2) we use our formulae in a predictor-corrector mode with (1.3) being the corrector and a con- ventional second derivative formula of the form k y , -y , , = h / B-f .+h 2 y,g , (1.4) 'n+k 'n+k-1 c_ M j n +j 'k 6 n+k j-0 being used as the 'predictor'. This ap- roach may be summarised in a step by step form in the following way: (1) Compute a quantity y , as the v t •'■'n + k solution of the conventional linear k- s tep f ormu la 29-1 : 'a.k" h3 k f n + k" h2Y k 8 n + k = y n+k-l k-1 + h y e . f ^ . j=0 (1.5a) (2) Compute y , , as the solution of v 7 n+k+l v -ha f -1i 2 yc ="7 'n+k+1 k n + k + 1 Y k B n+k+l y n + k It has been shown by Enright that for any integer k satisfying 1 < k < 7 it is possible to derive stiffly stable schemes of the form (1.4) having orders k+2. It is straightforward to verify that for all positive integers k it is possible to derive formulae of the form (1.3) having order k + 3. We now state the. following lemma regarding the order of accuracy of the procedure defined by steps (1) - (4) above. k-2 + he, ,f(x . ,y ( "£) + 2Z h6 - f .-., (1.5b) k-1 n+k n+k * j n+j+1 j-0 (3) Compute ? n + k+1 = f ( X R+k+1 ,7^^ ) (4) Compute y , from (1.3) written in the form y n + k" h Vn + k- h ^k g n + k = y n+k-l k-1 Z j=0 hg.f .+h6, ,f , , j n+j k+1 n+k+1 (1.5c) We note that in implementing steps (1), (2) and (4) it is necessary in general to solve systems of nonlinear algebraic equations for the required solutions. In each case these are solved using a modi- fied form of Newton iteration iterated to convergence. If the usual procedure of 8 2 f neglecting the term - — 7 is adopted, the coefficient matrices associated with steps (1), (2) and (4) respectively are he. J „. - h- k n + k Y k J n+k I - hB, J , , - h 2 y, J 2 , , k n+k+1 'k n+k+1 and I - h?, J , - h 2 y, J 2 , k n+k 'k n+k where J . is some suitabl n+k tion to the relevant Jacob x , . Exhaustive comparis n + k r earlier for the case where shown that if, with step ( coefficient matrix (A) rat then the iteration scheme variably converges if the scheme in step (1) converg experiments seem to indica also the case when y + 0. it has been found that ite (2) normally converges if matrix (A) is used. Thus only coefficient matrix wh compute and LU decompose i leads to a substantial sav .1 ■• f f r t . (A) (B) (C) e approxima- ian matrix at ons conducted Y, = have k 4 ) , we use the her than (C) , in step 4 in- i te rat ion e s . Numer ic a 1 te that this is Fur thermore ration s cheme the coefficient in practice the ich we need to s (A) and this ing in compu t a- L e mm a 1 Given that (1) Formula (1.5a) is of order k+2 (2) Formula (1.3) is of order k+3 (3) The implicit algebraic equations , r ■ ■ — (n) , — (n) , defining y , and y , , are solved usins n+k n+k+1 an iteration scheme iterated to conver - gence , then scheme (1.5c) has order k+3. This lemma is easy to verify by standard arguments. Having defined our basic formulae we are now in a position to examine their regions of absolute stability. We do this by applying them to the usual scalar test equation y =Ay, where X is a complex constant with negative real part. It was found that these algorithms are very highly stable for 1 < k < 6 and the cor- responding regions of A (a ) -s t ab i 1 i t y are given in Table 1. Note in particular that these formulae are L-stable for or- ders up to and including 6 whereas En- right's conventional second derivative formulae are L-stable only up to order 4. If we consider still higher values of k we derive A(a)-stable formulae with the angle a decreasing to zero as k is in- creased. However, we have not investiga- ted such formulae since, as stated pre- viously, we are interested in deriving a package where all of the formulae are very highly stable. Furthermore, when solving stiff systems of equations it may not be desirable to have too many differ- ent order formulae as quite often a change in order results in the coefficient matrix of the Newton scheme having to be re-evaluated and LU decomposed and this is, of course, generally very expensive. The exhaustive testing of a new algorithm is a long and complicated busi- ness and at the present time we are un- able to draw any firm conclusions con- cerning the general efficiency of our algorithm. However, preliminary results obtained using a fixed order (and some- times a fixed step) seem to indicate that the formulae discussed in the present paper are promising, especially when com- pared with conventional second derivative formulae. It is hoped that further re- search, particularly into choosing the 29-2 order anJ sccpsize ot' the particular for- mula being used, will lead to the produc- D efficient code incorporating the algorithms described in this note. TABLE 1 Order Region of A ( a ) -s t ab i 1 i t y ,o 90 90 90 89' 87" 83 v 29-3 COMPUTER STUDIES OF PLANETARY-TYPE EVOLUTION Donald Greenspan Department of Mathematics University of Texas Arlington, Texas 76019 John Collier Computer Science Department University of Wisconsin Madison, Wisconsin 5 3706 ABSTRACT In this paper a new, computer approach to the study of the interactions of particles with differing masses is applied to the study of planetary type evolution. The formulation contains an inherent self-reorganization property in which particles self-stratify in accordance with their masses. Computer examples are described and discussed. 1. INTRODUCTION Recently [4] , a new computer approach was devel- oped for modelling the interactions of particles with differing masses. In this approach, matter is approximated by finite sets of particles, each of which is treated as an aggregate of molecules, and the active forces are appro- priately rescaled. Inherent in the formulation is a natural, self -reorganization property, in which the particles self-stratify in accordance with their masses. After a precise mathe- matical and physical extension of the method, so as to include both long-range and short- range forces, we will apply it to the study of planetary-type evolution. 2. BASIC MATHEMATICAL DEFINITIONS AND FORMULAS For positive time step At, let t kAt, 0,1,2 For i = 1,2,3, . . . ,N, let particle P. have mass m. and at time t, let P. k j be located at r. , = (x. , ,y. , ) , i,k i,k i,k have velocity v. = i,k (v. , ,v. ) , and i,k,x i,k,y have acceleration a . , = (a . , ,a. , ). i,k i ,k,x i ,k,y Let position, velocity and acceleration be related by the "leap-frog" formulas ([3] , p- 107): ■* . At -* V i,l/2 = v i,o + T a i,o i,k+l/2 i,k-l/2 (At) i,k+l r i + (At) v i,k+l/2' (2.1) a. , , k = 1,2,. i,k . (2.2) , k = 0,1,2.. (2.3) If F. . is the force acting on P. at time t , where I l,k (F. . ,F. , ), then we assume that i,k,x i,k,y ■ and acceleration are related by F i,k i,k" (2.4) Once an exact structure is given to F . , , the i,k motion of each particle will be determined recursively and explicitly by (2.1)- (2.4) from prescribed initial data. In the present paper, we will want F . to incorporate both short-range and long-range components, and this will be implemented as follows the distance between P. and P At time t, , let r. . , be k 13 ,k Let G. . 3 13 (coefficient of molecular-type attraction), H. . (coefficient of molecular-type repulsion) , 13 (exponent of molecular -type attraction), a. . (exponent of molecular-type repulsion) and G* . (gravitational constant) be constants determined by P. and P., subject to the constraints (see 1 3' [6]): G. . > 0, H. . > 0, a.. > B. . > 2, G* . > 0. 1] - _13 - _ 13 - 13 - 13 - Then the force (F ; ,_ .. ,F^ u tr ) exerted on P . by i,k,x' i,k,y P . is defined by -G. F. H. . i /k,x ■ij,k 13 r. . 13 13 ,k G*. (m.m.) (x. ', - x. J-3 1- 3 J-.k 3>' ij ,k J r ij,k (2.5) i ,k,y JJ_ 13, k -G. JJ_ 13 ,k JJ_ r. . 13 13, k (m.m. ) (y. . 1 3 i,k y 3,k> (2.6) •ij,k The total force (F . , ,F. . ) on P. due to all 1 ,k,x 1 ,k,y 1 the other (N - 1) particles is given by N i ,k,x I»< 3=1 i ,k,x ' i,k,y N j=i i ' k ' y 3Vi (2.7) The formulation (2.1)-(2.7) is explicit and economical though nonconservative. Conservation of energy and momenta can be achieved [3] , but 30-1 the rule light-side particle: , -v 1.001 v. , i,k i,k (3.1) =1 Fi». l. only through an implicit, less economical approach. Throughout, the time step to be used in (2.1)-(2.3) will always be it = 10~ and a comprehensive FORTRAN program for implementation of (2.1)-(2.7) is given in the Appendix of [5]. It should be noted that the way in which (2.1)- (2.7) will be applied lends itself directly to parallel computation also [2] . 3. BASIC PHYSICAL ASSUMPTIONS AND DEFINITIONS We will consider a system of 137 particles, so that N = 137. This parameter was determined solely by economic constraints. The initial particle positions will always be those shown in Figure 1. Next, while in motion we will require a rule for determining the physical state of each particle. Such matters can be exceedingly complex (1) so that, for the present, we will use, as in [5] , the following intuitive notion. We will use "temperature" as a phenomenon of a particle's "local" velocity, that is, its velocity relative to neighboring particles. Thus, when a particle is rotating within a large system, the gross system velocities should have no effect on the particle's temperature. As a last consideration, we will allow for radiation into and out of the system. This will be accomplished as follows. Consider the vertical strip regions x_^(k + 1)y; y>0, k = 0,±1,±2,±3,... . In each strip, the particles with maximum and minimum y coordinates are called outer particles. An outer particle will be called a light-side_ particle, that is, it faces a sun, if y. . >_y.. Otherwise, it is called a dark-side particle. Light-side particles will receive radiant heat while dark-side particles will emanate radiant heat. This is implemented as follows. Every k* steps, outer particle velocities are reset by dark-side particle: EXAMPLES v. , ♦ 0.9955 v i,k i,k A variety of examples were run with various combinations of parameter choices selected from 6 = 4,5; E = 5,10; k* = 50,100; D = 2.1,2.3; Y = 0.1,1.0,2.0; and G*. = 0.01,0.001. From these we will first describe two, in both of which the parameter choices are 6=4, e = 10, D = 2.3, k* = 100, y = li and G*. = 0.001. The results of the variety of other examples run can only be summarized easily as follows. Only minor changes in input parameters often resulted in dramatic changes in dynamical behavior [4] . This is, of course, consistent with the diversity one observes among the planets and moons of the solar system. Other computations, aimed at duplicating the development of galaxy arms, all ended in failure. The systems invariably self-reorganized into one or more relatively elliptic or circular subsystems. This suggests to us the strong possibility that more than gravitation is involved in galaxy arm formation. Finally, note that initial calculations with a 239 particle configuration were begun, but these had to be discontinued due to a lack of funds. REFERENCES 1. BARKER, J. A. and HENDERSON, D.: "What is "liquid'? Understanding the state of matter", Rev. Mod. Phys., 48, 1976, pp. 587-671. 2. FITZWATER, D.R.: "The Formal Design and Analysis of Distributed Data-Processing Systems", TR322, Dept. Comp. Sci., U. Wis., Madison, 1978. 3. GREENSPAN, D.: Discrete Models, Addison- Wesley, Reading, Mass., 1973. 4. GREENSPAN, D.: "Computer studies of inter- actions of particles with differing masses", Jour. Comp. and Appl. Math., 3, 1977, pp. 145-154. 5. GREENSPAN, D. and COLLIER, J.: "Computer studies of swirling particle fluids and the evolution of planetary-type bodies", TR 299, Dept. Comp. Sci., U. Wis., Madison, 1977. 6. HIRSCHFELDER, C.F.; CURTISS, C.F. and BIRD, R.B.: Molecular theory of gases and liquids, Wiley, N.Y., 1954. 30-2 On the Reliability of Error Estimators for Explicit Runge-Kutta Procedures J.H. Verner Department of Mathematics and Statistics Queen's University, Kingston, Canada Abstract . Arbitrary parameters for (s,p-l(p)) Runge-Kutta procedures should be selected so that resulting codes yield reliable results efficiently. Shampine suggests as a compromise a process for selecting parameters which attempts to optimize efficiency subject to a constraint designed to provide reliability. This process is adapted to choose parameters for procedures with p = 5,6, and some comparisons with known procedures are given. 1 . Introduction Currently (s ,p-l (p ) )-procedures are the most efficient Runge-Kutta processes available. Pairs of formulas of successive orders of accuracy provide, for example, an approximation of order p-1 and an error estimate (EE) of order p to the local error (LE). It is known that families of these procedures may be represented in terms of arbitrary parameters, yet there is little agreement on how they should be selected. Production codes are implemented so that some function of EE is bounded by a user specified tolerance. A basic assumption in using such a code is that the behaviour of the global error is reflected by that of the local error. In -.howing that this occurs for various strategies, Shampine £3] needs to assume that the actual local error is bounded by the specified tolerance. Accordingly, parameters should be selected so that EE reflects the behaviour of LE. A second more prevalent criterion for selection is that of minimizing coefficients of the pi iiu;ipal local truncation error in order • '- increase efficiency. To enhance both tcy and efficiency Shampine [2] in'. 1 compromise , a constrained optimization of the parameters. This process is modified below in order to make the computation feasible for higher order formulas. Thus we propose to derive reliable (|LE-EE| is small), efficient (|EE| is small), accurate (|le| is small) procedures. Because the optimization can- not be viewed as minimizing these quantities, the resulting formulas cannot be considered optimal. While potentially useful formulas may be derived, the analysis is more appropriate for indicating which procedures should be avoided. 2. Selection Criteria For elementary differentials (D . } evaluated at the current approx- qi imation, error expressions in the next step are LE - hP N I S T P-1 T s h q r t^d i=l pi pi q=p+l i=l Q 1 and N N 00 c + S h q £ i=l Pi Pi q=p+l i=l ! 1 EE=h P S T P D (T p - 1 -T P )D . qi qi Q 1 32-1 where _[' and _T P are vectors of nunerica] conatanta determined i>\ the ■ethodfl of orders p-1 and p respert- i\cl\. While for h small, there will be little "contamination" of the error ea timet e, efficient solutions for moderate tolerances will require moderate to Large stepsizes. Hence, the contam- ination may be substantial, for example, mv coefficients {r } are zero, or pi i i t lit- magnitude of T p is substantially Larger than that of T^~ These comments indicate that the reliability and efficiency of a procedure may be assessed by considering the parameters a = max | T^ 1 | , i P 1 I T?~ _ . | = max | T p 7 | , q^p , q-li qi y max | T p ~ | = max T P <1 £ qx i qi It 13 : 1 -/. [ qx qi "q 4 T P"1 > qi z = number of coeffici< q jnts q>p q>p j JP -1 = qi Following the intent of Shampine's proposal, new procedures are derived to minimize a subject to making e p+l' p+l» T p+1 sma11 << T if Possible) and z = . (Even this optimization P may not yield reliable procedures, for example, if T p , = T^ . for some qi qi indicies. ) 3. Derivation and Comparison A four-parameter family of (6,4(5)) procedures is known [5]. If the nodes are confined to [0,1] and >6± *» a is minimized if c, m i±£E 5 A choice of c - so that 5 minimi zea 6,, and a further choice of Cg minimizes ^5 • The procedure V3 has nodes (0,-,^—,^-,-^-, QQ V ) . This is contrasted to other procedures : I - Fehlberg [l],Vl [4], V2 [5~|. The fact Procedure z 5 P ^6 Fl F2 V2 V2 V3 2. 08 (-3) 1.67 .20 1.9 2 1.28(-3) 2.59 .55 .46 8.57(-4) 3.77 .67 .36 5 6.25(-4) 4.44 1.00 .62 1.87(-3) 1.00 .33 1.16 Table 1: Parameters for (6,4(5)) procedures Procedure FCN STEPS MAX ERR FRAC DEC Fl 134925 21970 1.3 .002 F2 120406 19541 10.4 .022 V2 110803 17903 10.6 .039 T V2 Disaster for y = -y V3 147757 24157 1.8 .002 Table 2: DETEST for ( 6,4( 5 ) )-procedures that V3 does not perform as well as Fl on DETEST indicates the deficiencies of this analysis . To obtain an (8,5(6)) procedure with j3_ < 1 is possible only for a Si 9. 2 (-4). As such a procedure was unlikely to be as efficient as DVERK, others were considered. Without con- straining $- , a is minimized by the , , n 4 8 20 9 n 1 , ■> .. nodes (0 ,—,— ,— ,— ,\ f —,\) to the value 1.98(-4); however 3 = 6.15 . As a more reliable alternate to DVERK we suggest the procedure (V3) with nodes rn 1 10-2/10 5V10 4 , 1 -n (0 '"o"'— H ' i5~'T jl 'iT ,1) ' Observe that the results of Table 4 are supportive of the indicators of Table 3. For p = 7,8, similar com- parisons are not so consistent: while results from DETEST indicate procedures of [5] are reliable and efficient, the parameters 3,Yj5, do not adequately 32-2 distinguish them from other known procedures . Procedure z 6 a P 7 ^7 5 7 F 14 4.b2(-4) 2.07 1.00 3. 50 DVERK 5 4.62(-4) 4. 61 .91 3. 50 VI 1.69(-3) .90 .15 1.40 V2 5 4.62(-4) 1.43 .69 8.01 V3 3.67(-4) 1.54 .60 2.85 Table 3: Parameters for (8,5(6)) procedures Procedure FCN STEPS MAX ERR FRAC DEC F 75091 8913 22.3 .063 DVERK 73206 8 603 15.3 .134 VI 9 2139 10995 2.5 .00 2 V2 80282 9568 10. 7 .022 V3 76184 9012 4.1 .022 TABLE 4: DETEST for (8 , 5( 6) )-procedures In conclusion, we can only advise that procedures with z =f 0, or any of .£> Xt .§ relatively large should be avoided. Acknowledgements : This research was supported by the National Research Council of Canada under grant A8147. REFERENCES 1. W.H. Enright, R. Bedet, I. Farkas, T.E. Hull, Tech. Report 68, 1974, Univ. of Toronto 2. L.F. Shampine, Report SAND 76-0370, 197 6, Sandia Laboratories L.F. Shampine, Appl. Math, and Comp . 3, 189-210, 1977. 4. J.H. Verner, SIAM J. on N.A. 15, 772-790, 1978. 5. J.H. Verner, Preprint 1978-19, 1978, Queen's Univ. 32-3 A COMPARISON BETWEEN SEVERAL PIECEWISE CONTINUOUS ONE STEP INTEGRATION TECHNIQUES J. P. HENNART AND R. ENGLAND Universidad Nacional Aut I I MAS (Institute for Research in Appl Apartado Postal 20-726, Mexic In a series of papers by the first author [ 1 -3] a class of one-step schemes has been developed which is particularly well adapted to the numerical integration of stiff systems of ODE'S. These schemes lead to piecewise continuous approximations in time with the added possibility of exponential fitting, and provide therefore a very general frame_ work for the development of finite elements in time in connection with parabolic evolution problems [A, 5] . In this paper, details of some of these schemes are given, for cases which do not require second derivative evaluation, and they are compared with similar alternative schemes providing piece- wise polynomial solutions on the one hand [6,7], or exponential fitting on the other [8-10]. Given a system of n ODE'S in autonomous form: y = f ly(t)), (1) "0 \. ^ onoma de Mexico ied Mathematics and Systems) o 20, D. F. (Mexico) the unknown parameters y. , r=1 , s=l,...,q-1, being finally e pl (t iple collocation at t=t. and t ; 1 y (i) fori and V (S) (t. +1 )=f (s - 1) (y (t. + J), s=1, o-pq i+l % %pq i+l .,P-1, iminat :t i+1 : ..P-1, and ed by (6. a) ,q-1, (6.b) with v (t.) known either from th tions \i=0) or from integration interval [ t . , , t .] . After such a su (3) becomes: e ove 1 n 1 1 1a r the bsti tu t . i+l *i+i = ^i + I % q ( lrh^' t] )dt 1 condi- prev ious t ion.Eq. , (7) t. 1 one gets after in teg rat ion over [ t . , t . , . ] : 1 1 +1 y U i+1 )=y(t.)+ f i+l f(y(t)) dt. (2) In the integral on the righthand side, y(t) is then approximated by its two-point Taylor interpolant y (t) , lead ing to ■\-pq i+1 i+1 = V; + [ '" f ( Vnn (t))dt > 1+1 *' j % o,pq (3) i.e. a system of nonlinear equations in the unknowns y, + 1 Eq. (7) defines a class of discrete one-step ^ integration methods, depending on the quadrature formula which is used to approximate the integral, in addition to the choice of p and Assuming that y exists and is unique, each of its components^alpends on (p+q) parameters which are the corresponding components , of y> r ^, r=0,...,p-1 (if p > 1), and of y|*{ ,^=0, . . . ,q-1 ( if q > 1) . Exp I icitly W l) p-1 I r=0 .Cr) t- t-t . v r! where p, q are nonnegative integers, while each component of y belongs to the 1 inear space -X/Pq S ={t n expU t) |n=0 N . ;m=l , . . . ,M; IN =p+q}, m m-1 m m (A) where the X 's are given distinct real numbers, one of them usually being zero, v is then chosen to satisfy the p+q interpolation condi- tions y (') ( t .) = y .( r >, r=0 p ., 'pq « ;i and (sj .(s) ' (t ) = y ! . S=Q q-1 ^jjq i+l 'J+l (5. a) (5-b) q-1 s=0 ^ (s) t-t. i+1 u^ ID = , r,£=Q, 6 r£ , r,£=0,...,p-1 .q- and v^(0) = vf (1) ■s;: , s,£=0,...,p-1 , s,£=Q,.. .,q-1 (8) where h= t.,.-t. and the functions u , v e S . , i+l 1 r' s sat 1 sfy , (9- a) , (9.b) , 19. c) , (9-d) so that (5) is automatically satisfied. Then, when 33-1 some k point quadrature formula is used, Eq . (7) becomes k y. , = y. + h E w,,f (y (t. -M3»h)) with < 0. <•• •< 0, < 1 , that is k (10) i+1 y- «i,-rf l ^l r) ". ( v^!^l:! ».««"*» ft*- s=ov (ii) Restricting attention to the case p,q <2 with srn , the space of polynomials of degree less thaR 8r equal to (p+q-1), (11) leads to a class of (k+2) stage Runge-Kutta schemes which can be characterized by the following array in Butcher's notation (and with p=q=2) : W Vo (( V Yo 10 ! 1 W (12) u (e k ) w lVo (0 k ) w k v o (0 k ) Vl (0 k ) „, • • • "k The special cases 0. = and (or) 0, = 1 , and those with p and (or) q < 2, are easily derived from this general class by noting the properties (9) of the f unc t ions u , v . r s A number of particular cases are of interest, among which the following may be mentioned. p=1 , q=0, and any quadrature rule exact for a cons- tant integrand, gives the forward Euler method. p=0, q=1 , and any quadrature rule exact for a cons_ tant integrand, gives the backward Euler method. p=2, q=0, and any two-point quadrature rule exact for a linear integrand, gives a 3 stage, second order, explicit Runge-Kutta method. The family with O.=0 gives all the 2 stage methods of that type. p=q=1 , and any two-point quadrature rule exact for a linear integrand, gives a 2 stage, second order, implicit Runge-Kutta method. The family with 0.=O consists of semi -expl ic i t methods, with Butcher's array: (13) ■1 1 2 >-k 1 which include the well known cases:©-, ■ y (implicit methodj I (trapezium metnod) . p=0, q=2, and any two-point quadrature rule exact for a linear integrand, gives a 3 stage, second order, implicit Runge-Kutta method. The family with 0-= 1 consists of 2 stage methods, but these show no apparent advantage, in terms of effort or stability, over methods of third order with p=1, q=2. p=2 , q=1 , and any three-point quadrature rule exact for a quadratic integrand, gives a 4 stage, third order, implicit Runge-Kutta method. The case with 0.= 0, 0_= 2/3 is a 2 stage, semi-explicit, Radau method, of higher order than those given by (13), but of limited interest because it has a bounded stability region. p=1 , q=2, and anv three-point quadrature rule exact for a quadratic integrand, gives a 4 stage, third order, implicit Runge-Kutta method. The family with 0,= 0,0,= 1 consists of 3 stage methods, amongst which should be mentioned the Radau and Lobatto methods obtained with 0, = 1/3, 1/2, 2/3 having Butcher's arrays: 1/3 5/12 -1/12 1 3/4 1/4 3/4 1/4 (20) 1/2 1/8 1/2-1/8 2/3 2/9 2/3 -2/9 1 1/6 2/3 1/6 1 1/4 3/4 1/6 2/3 1/6 1/4 3/4 the first showing possible advantages both for storage and convergence during solution of the nonlinear equations (7) at each step. p=q=2, and any three-point quadrature rule exact for a cubic integrand, gives a 5 stage, fourth order, implicit Runge-Kutta method. The case with 0=0 ,0, =1 is a 3 stage, Lobatto method with Butcher's array: (21) 1/2 5/24 1/3 -1/24 1 1/6 2/3 1/6 1/6 2/3 1/6 an A-stable method which has been given by other authors. At this point, it is interesting to point out that the "optimal" Runge-Kutta schemes based on Gauss- Legendre quadrature [6] do not appear as special cases of the above formalism. These optimal schemes are more finite difference oriented in the sense that they achieve superconvergence at the mesh points. However, as is well known, they lose accu- racy between them. Our schemes, which were studied in detail in [ 1] for general p and q, do not pres- ent that disadvantage: order p + q is achieved not only at the mesh points but also between them, under the only assumption that the quadrature rule is exact when appl ied to the members of S. The basic idea for developing such a formalism was to provide a piecewise continuous representation in time that could match in a fully consistent way the high order accuracies offered by finite ele- ments in space for parabolic evolution problems. 33-2 It turns out that in the polynomial case, the stability properties on I > depend on p and q (if the quadrature scheme is of sufficient accuracy): namely the schemes with p=»q are A-stable, while the schemes with q=p+l . p+2 are strongly A-stable. In the general case where exponential basis functions are present (Eq. (M), it is possible to show [5] that the above schemes when applied to the test equation y= \y are equivalent to the multj_ point Pade approximant schemes interpolating the exponential N times at x=h ^ m « m= ' M > and at least once at m x=0. These schemes have been studied and compared with other exponentially fitted sche- mes [8-10] in reference II in relation with systems of nonlinear ODE'S and some of them applied in ref- erence 12 to the nuclear reactor point kinetics equations. To be applicable, they required that the equations (if nonlinear) be linearized, and that mean values of the parameters be taken: both approx imations introduce third order local errors, and the corresponding schemes were therefore limited to second order. This is not the case for the new schemes presented here which retain order p + q (or at least p+q-j . if j is the number of inter- polation points that remain far from the origin) for any problem. As far as stability is concerned, all the results previously obtained for the multi- point Pade approximants are valid here [8-12]. For a general space S of approximants as given in CO, it is again required that the quadrature rule be exact for its members. In the presence of exponential basis functions, this problem is less documented than in the polynomial case, although it is known that Gauss- type formulae exist for such Chebyshev systems [13l. Since very few realistic problems would require fitting of more than one exponential node, an attractive alternative may be described as follows. First multiply both members of (1) by z(t) to get z y = z f (y(t)J, or equivalently (22) (z y) z y + z f(y(t)). (23) If (23) instead of (1) is integrated over [t.,t. + .] and y replaced as above by y in the right nandside of tt»e resulting express ion^, i t is easy to realize that if z(t) is chosen to be proportional to exp (-At), 1 ^0, the components of" y behaving as exp (At) will be integrated exactly independently of S (and of the associated quadrature rule). Apart from (i) a family of schemes for advancing the solu- tion of Eq. (1) from t. to t. + ., a practical implementation also requires: (ii)a method for estimating the local truncation error; (iii)a strategy for selecting the step length h, and possibly also the order or particular member of the family to be used. In addition, any implicit scheme involves the solu_ tion of a system of nonlinear equations such as (7), and the implementation also requires: (iv) an a Igor i trmfor solving a system of nonlinear equations, while if output is required tit points much more closely spaced than the grid points used for the integration: (v) an interpolation process should ideally be provided. Some comparisons have been perfor_ med on complete implementations, without separating the five aspects mentioned above [ 1^1 . However, in relation with this work, a modular structure is being developed with each aspect treated in a separate routine, thus making it easy to keep the elements (ii), (iii), (iv) constant, and vary the schemes (i) used -(v) depending upon (i). It will also facilitate a later study of the effects of varying (ii), (iii) or (iv). Detailed run statistics will be reported later. REFERENCES. 1 . 2. 5. 6. 7. 9. J. P. Hennart, "One-step piecewise polynomial multiple collocation methods for initial value problems", Hath. Comp . 31, 24-36 (1977). J. P. Hennart, "Piecewise polynomial approx- imations for ngclear reactor point and space kinetics". Nucl . Sc i . and Engng .64, 875-882 (1977) . J. P. Hennart, "Some recent one-step formulae for integrating stiff systems of ODE's", Se- minaire d'Analyse Numerique No. 253, Universite Scientifique et Medicale de Grenoble (1976). J. P. Hennart, "Multiple collocation finite elements in time for parabolic evolution pro- blems", to appear in The Mathematics of Finite Elements and Applications I I I -MAFELAP 1978 , J. R. Whiteman, Ed., Academic Press, London. J. P. Hennart, "Finite elements in time and space for transient reactor neutronics analy- sis", to appear in Proceedings of the National Topical Meeting on Computational Methods in Nuclear Engineering , American Nuclear Society (1979). B. L. Hulme, "Discrete Galerkin and related one-step methods for ordinary differential equations", Math. Comp . 26, 881-891 (1972). K. Wright, "Some relationships between implicit Runge-Kutta, collocation and Lanczos t methods, and their stability properties", B IT 10, 217-227 (1970). W, Lini-ger and R, A, Willoughby, "Efficient integration methods for stiff systems of ordi- nary differential equations", S I AM J. Numer . Anal . 7, 47-66 (1970) . B. L. Ehle and Z. Picel, "Two-parameter, arbitrary order, exponential approximations for stiff equations", Math. Comp . 29, 501-511 (1975). 10. F. H. Chipman, "A note on implicit A-stable R-K methods with parameters", B I T 16, 223-227 (1976). 33-3 11. J. P. Hennart, "The multipoint Pade approx imant approach to exponentially fitted one-step forrnu_ lae", Comun icac iones Te"cnicas, Vol. 6, Serie 3, No. 119, ClflAS-UNAM, Mexico (1975). 12. J. P. Hennart and B. Torres Barrios, "Generajj zed Pade and Chebyshev approximations for point kinetics", flucl . Sci . and Engng . 59, 170-187 (1976). 13- S. Karl in and W. J. Studden, Tchebycheff Sys - tems: With Applications in Analysis and Sta - tistics , John Wiley S Sons, New York ( 1 966) . \k. W. H. Enright, T. E. Hull, and B. Lindberg, "Comparing numerical methods for stiff systems of O.D.E. 's, BIT 15, 10-48 (1975) . 33-4 ON THE A-STARILITY OF NUMERICAL HMLTISTEP METIinnS r n» SOLVING IVPs IN ODEs BY HASSAN NASR "o^artnent r> c ''athenatics, military Technical Collepp Cairo 1079 tbstrnct : is well known (Dahl^uint 1163) t for the studv o c the A-stability of linear Multislep Methods for solvinp IVPs 3Es, A-stable method for order hipher than two do not exists. This paper pives the necessary and sufficient condition for the general multistep methods to be A-stable. It includes an example f or a new A-stable -lethod o c order u. Introduction : Let us recall the definition of the "eneral Blocking Over- Impl ic it '-'ultistep "1 "ethod, by Hassan Kasr in Aolikace "athematikv (under printing). The method is used r or solvinp the ODE jr=f(x,y), xe[a,b] with the inital condition wher cont cond that the in it in t this one solu k po s po e th inuo it io the solu ial he w met step t ion int s ints a+i v(a) = e f unc us and ns w . r exist t ion o condi t hole [ hod li of th is co sucpo (1 < mh, j = t ion sat . t .y ence f eq ion a,b] es i e me mput sin? s <_ 0,1, f(x isf v in and uat i (2) . T n th thod ed s tha k). ,y) inp [a,b uni on ( be p he e e fa , th imul t it Let wher is d the ]x(. nuen 1) u uara ssen ct t e ap tone is the e m (1) (2) ef ined , Linsch it z 00 ,°°) so ess of nder the nteed ce o* hat in proximate ously in known in po ints is a positive integer am! h > is the integration step, will be called basic points and the points x.. .=x +w.h, ] k + 1 ■) k ! i=l,...,k-l where w. are any real numbers, will be called intermediate points while the approximate solution the point x., . ( i = , . . . , k-1 , ) will be denoted bv y(x. ), lk+ i 1=0,1, the s-d iment ional vector ik" * * ' X ik + s-l ) and the s-dimentional x . = ( x - 1 vector of the approximate solution at these points bv VC x . ) = (v'C x . , ),..., y(x., , )) and analoricall v the vector "I k + s-1 z . and v( z . ) by -1 -1 — j i k + s ' ' X (i+l)k+s-l ) and y(0=(y"(x ik+s ),...,r(x ( . +1)k+s _ 1 ))' and usinp the vector- valued functions f(x./v(x.)) and ^(z., v" ( z . ) ) at the —1 —1 —1 - -1 correspond inp points. If B,D be any kxs matrices and C any k*k matrix of constants, then the system: ^(z. )=B ^(x. )+hC f(z..,y'(z.))+hn f(x.,f(xJ) •1=0,1,. .. (3) represents an alporithm o f the so called HBOM method. We define the local trunction error of the method by the k-d imens iona 1 vector L(v(x),h)= y ( s+w s h ) y(x+w k-l h) y( x+(m+w )h ) y(x + (m + w , )h ) L/ s _ i _1 -hC v ' ( x+w h ) s V ' (x + w, h ) k-1 V ' ( x+ (n + w Q )h ) y * ( x+(m + w ,)h) v( x+w h ) o y(x+w _,h) ■ hn y ' ( x + w h ) o v ' ( x + w , h ) _ s- 1 (4) 34-1 If the c irst k-s components of the local truncation error is of order lv while the remaining s component are of order hP +1 we can rive a .weak definition o^ the order o c the method as the order p(p _> 1) w.r.t.s if and only if the following kp+s eauations are satisfied. s . £ h • • =1 j=l H for i=l, ,* (5) the necessary and sufficient condition f or the convergence of the solution o f the iterative system of equations (10) i.e. |w(z) | 0(z) i;L (z) , (z) si ■ n («) ,0 (z) ss < 1 s ' " i • - b ,• a s-l+i j=i il k-s n q r y o-l 'i-l =qC i=l w S-l + 1 for Re z < Example for a GBOM Method : (11) k i s ^ c . . ( m + w ,-,.•)" + , n s-k-l+i 1=1 ii -i-i i=k- s +l for i=l,...,k-s, q=l,...,p-l : d,_. w q_1 ] (6) Takinr k = 2, s=l, m=l, w =0, w =£ o 1 2 and solving equations (5), (6) and (7) as 9-equations in 8-unknown parameters i.e. the method is of order 4 and it can taVe the form ( m + w ) - - E b . . w . s-k-l+i j = i n 1-1 y = v + — ( 5f + 8f -f ) n+1 n 24 n n+1 n+2 (12) k-s a-1 =q[.^ c..w + c..(m+w . . i=x n s . 1+i i= k _ s+1 it s-k-i+- o-i y o = y + - (f + i-f +f „) ^ J n+2 J n 6 n n+1 n+2 The A-Stability of the Method : ,q-i +.Z d.. w q-x ] , for i=k-s+l, . . . ,k, i=l il i_i q = l, • . • ,p (7) A-Stability Condition for a GBOM Method : Solving the equation y ' =ciy , Rert< 0' eauation (3) can take the form (1-z C jy( z . ) = ( B + z D)y'(x.), where z = h 1 1 (I-zC) 21+1 21+2 v 2i+l y 2i+2 = (B+zD)y 12 •6z + 12 2j 2 -i z z^+6z+12 21 The last equation is (13) (14) *y(z.) = (I-z C)" 1 (B+zD)y'(x.) (8) If 0(z) = det (I-zC) The lest s equations of this vector is ^21+2 z'+6zt !2 z 2 -6z+12 y 2 .=»U)y 2 . (15) where w(z) is the complex variables f unct ion ik + k 0(z) _ y jk+k + s-l Or the form Q,,(z).. .0,_(z) 11 In + l '.'(z) Y Is . (z) . . .n (z) si ss ik ik+s-1 (9) do) Uaine the fixed-point theorem, it is clear that the necessary and sufficient condition for the A-stability of GBOM method (3) is , N z +6z+12 „ . n r(z) = Re z < z 2 -6z+12 (16) Applyinp the maximum modulous principle. Then the modulous of this function attains Its maximum value at the boundary (i.e. at z = iy , z=«>) which fives |w(z)| < 1 from (10), (15) and (17) the methods (12) which is of order 4 is A-stable. (17) .34-2 e'erences: 1) Dahlnuist, R. "A foecial Stahilitv p roMe"i f or linoar Multisten Method : T 3 (1063) 27-H3. 2) "ear, T.H. "Numerical Initial Value Problems in Ordinary "^inferential Equations "Print ice- Hall , Inc. (1071) 212-222. 3) Hassan N'asr. "Generalized Blocking Over- implic it Multistep Method" APLIKACE HATEMATIKY, Under Printing (1079) . u ) Hassan Ilasr. "Necessary Conditions for the Generalized Dlockin" Over- inplicit Multistep Method. "APLIKACE "MATIKY Under Printinr (1079). 34-3 Local error estimates for linear methods for ordinary differential equations G.J. Cooper University of Sussex 1 . Introduction Most methods used for the numerical solution of ordinary differential equations may be characterized by a pair of matrices (A,B) . Here, local error estimates for such methods are obtained as linear combinations of computed values. (This is equivalent to constructing pairs of imbedded methods.) The local error estimates depend on the assumption that exact solution values are used to start the step but it is shown that this assumption is not required for certain estimates. Some hybrid methods have such independent estimates, and also have simple step length change procedures. Consider an initial value problem x' = f(x) , x(t ) = x , f : R n -> R n , and suppose that approximations are required on [tg,t()+H] . Choose a step length h so that Mh = H , M a positive integer, and let t = tn + mh , m = 1,2,...,M . A linear m v s-stage method computes, in step m , s n-element vectors, (m) s y i a-.y^-h I b..f(y< m) ) , j=l l J J j=l lJ J ,(m) for i = 1,2,... ,8 , where y> ' approximates x(t + c.h) and c = (ci ,c-> , . . . , c ) is a m-1 i ' ' * ' s consistency vector. The method may be expressed as Y (m) = AY(m-l) + hB^On)) , (1) where Y (m) = yi (ln) e y 2 (m) ». ..•y< n) and F(Y (m >)= f(yf m )) 9 f(y£ m) ) • ...«► f(y( m >) are elements of R ns . The linear operators A and & are defined by matrices A and B , the rows of these matrices determining linear combinations of the s components of y^" 1- *-' and F(Y (m) ) . If step m is started with exact solution values _(m) for s s (m) E a..x(t ,+ c.h) + h I b..f(z) ') , j-i u m-2 = 1 !J ],2,...,s , and this may be expressed as (2) z (m) _ AX (m-l) + hBF(z (m) } _ The local error in step m is X (m) Z but the essential local error is defined to be .w /v .(m) ~(m). , . . . M A (X - Z ) where w is a given integer For a Runge-Kutta associated with the method method , . .. 1 1 ° 1 A = ... 1 and w = 1 so that the essential local error is merely the error in the final stage (repeated s times). For predictor-corrector pairs using methods of equal order, w = , but because these are multi-step methods some components of v (n) (m) X - Z are zero. The aim is to determine a matrix Q so that, (m) for the corresponding linear operator, QZ estimates A w (X (m) - Z (m) ) . For a Runge-Kutta method the rows of Q will be identical. For predictor-corrector pairs some rows of Q will (m) may be That is, some estimators be zero. For some estimates (JZ replaced by QY (m) Q are independent of the localizing assumption ™,(m) . »w.„(m) „(m). and QY estimates A (X - Z ) . 2. Conditions for error estimates T Let p = (pi, P2,...,p ) be a given vector with positive integer elements and define hi ->n . mfx M Pi ||y.|| , Y = viGyo*. . .toy 1 M l " l " ' s where is a norm on Suppose a 35-1 method is used to compute Y ,Y ,(M) (which depend on h and hence M ). If (for suitable starting values) ||X - Y m || < k , m = 1,2,...,M , for all sufficiently large M , M i M' , then l|x(t .+ ch) - yf m) || £ Kh Pi , i = 1,2 s , m-1 1 i and the method is said to be order p convergent. A method is order p consistent if for M >, M' * K , , x (m) _ z (m). M , A w (x (m) ,(m). K m = 1,2 M , ll M * M ' for some integer w i . (A more general definition is possible [3].) A stable hybrid method is order p convergent if it is order p consistent. In the present general context, a hybrid method is defined by , max Ae. i = , i = 1,2,. . . ,s , is the natural basis for of A , this where ej , e 2 , . . . ,e s . . R . Since Ae. is column l l means that stage i in step m-1 is used in step m only if this stage is of maximum order. To estimate the essential local error, a matrix Q is required so that (for M >, M') | A w (x (m)_ Z M ) _ QZ (m), M £1 M 2 1,2,. ,M Since A W (X (m) - Z (n °)- QZ (m) m .w> ,,,(m) „ (m) s ,.„(m) ... . (Q + A )(X - Z ) - QX V , sufficient conditions for Q to be an estimator may be obtained directly from the algebraic conditions for order p consistency [2]. Let T e = (1,1,..., 1) and let C be the diagonal matrix whose diagonal elements are the elements of the consistency vector. Theorem 1 Suppose a stable hybrid (A,B) method is order p consistent for a given • ^i_ ^ min order vector p with p. < 2 . p. , J J_ i = 1,2,. ..,8 . Then Q is an estimator if, 'for each i = 1,2,. ..,s , eT QC T e =0, t =0,1,..., p. +1, t w T ,r 1 T < -1 e! fO+A w ) BC " . BC [C T ° - A(C-I) T ° Bc/^e Tq+Tj ♦...+ t £ p. + 1 , where u = 0,1,2, and t r; , t , , . . . ,t take all possible positive U 3. Independent estimates An estimate is independent of the localizing assumption if k' || Q(Z (m)_ yW)! Let U (m) = Z (m) - Y (m) , m = 1,2, ,M . Then (1) and (2) give U (m) ^A^U (0) ♦ m i l L l CL<*+>- Z^) + i=l ~ h \ A i BrF(Z m i} i=0 F(Y (m_i) )] Let m > v , a positive integer, and assume QA =0 so that the expression for QU simplifies. The function differences may be replaced by a derivative expression and (1) applied again. With further assumptions it becomes possible to give conditions for independence. The following result is a particular case. Theorem 2 Suppose the conditions of theorem 1 hold. Suppose that QA = and QBA = , v a positive integer. Let p. = p i = 1,2, ... ,s , and T x -1 Tj-1 QBA BC v . . . BC r L n [C °- A(C-I) where t = 1,2,3,..., t BC ]e 0,1,2,. , and T 0» T 1» take all positive integer values. Then (for suitable starting values) IQ (Z (m) - Y (m) )|| S 4 . m * v . M 5 M' . ii M Such an estimator is called a Q(v) estimator. The condition m >, v is crucial because the estimate is independent of the localizing assumption only after v steps have been completed (and this applies to each step length change). The new starting conditions have to be chosen to maintain order p convergence. Consider the Adams-Bashf orth and Adams- Moulton methods of order k used in the PECE mode. A Q(k) estimator may be obtained and this is just the Milne estimate. An extension of theorem 1 provides estimates for Runge-Kutta methods of order four (including the classical fourth order method). These estimates are obtained by imbedding, using three extra stages 35-2 (two extra function evaluations), but may not be independent. The following method, represented as p|A|B|Q|c,isa hybrid method with a Q(l) estimator (and w » 1). The first three stages give a hybrid method requiring two evaluations of f in each step. Since the second stage is of order three, and is not 3 3 1 I 3 u 4 1 8 "7 7 5 16 21 21 1 8 T 7 used in the next step, it may be used to halve the step length. The final two stages are needed only for error estimation and one function evaluation is required in the next step; overall one extra evaluation of f is required. A similar method of order four is known. 3 4 2 7 8 7 452 567 184 567 8 81 8 ~7 16 7 3 7 27 7 -1 3 ~4 _16 21 16 21 The method is also given by the equations (m) (m-1) yi - Y3 . (m) 1 (m-1) 3 (m-1) 3 ... (m), y\ = i yi + i y\ * ■% hf(y} ') , 1 (m-1) 8 (m-1) 2 ... (m), x 8 ... (m) , y3 s -,yi + - y3 - ■? h * (yi ) + - hf ( ?2 ) , 21 1 7 7 "■ 7 J3 7 * w * ' 7 5 (m-1) 16 (m-1) 452 ... (m), 184 ... (m), J 8 ... (m), y* - tt y\ + rry3 + TA7 hf( yi )-^r hf (y2 ) + -5rhf(y 3 >, 567 81 567 1 (m-1) 8 (m-1) 8 ... (m), 16 ... (m), 3 ... (m) , 27 U£/ (m) , '4 y y§ '- - hf(yj ) - — hf(y 2 ) + j hf(y^ ') + — hf(y^ y ) Since 21 = x(t ) each component of the essential local error is proportional to x(t ) - Z3 m— 1 m and the estimator gives just one estimate „/,. s _ (m) „(m) x(t m ) z 3 = y 5 yi m) + 0(h 5 ) J.C. Butcher, On the convergence of numerical solutions to ordinary differential equations, Math. Comput., 20 (1966), 1-10. G.J. Cooper, The order of convergence of general linear methods for ordinary differential equations, SIAMJ. Numer. Anal., 15 (1978), 643-61. R.D. Skeel, Analysis of fixed step-size methods, SIAMJ. Numer. Anal. 13 (1976), 664-685. 35-3 IN SEARCH OF A ROBUST INTEGRATION ALGORITHM FOR GENERAL LIBRARY USE : SOME TESTS, RESULTS AND RECOMMENDATIONS M.B. Carver Atomic Energy of Canada Limited Chalk River Nuclear Laboratories Chalk River, Ontario KOJ 1J0 ABSTRACT The design or selection of an integration algorithm suitable for general use in a subroutine library or within a simulation package is a demanding process, particularly as it is unrealistic to support a number of routines which purport to perform the same task. In addition to satisfying the criteria of excellence of documentation, reliability and efficiency which a numerical analyst would require of an integrator to be used in a specific project, a routine destined for general use must be exceptionally robust, as the average user has neither the desire, nor the specialized knowledge to interact effectively with the algorithm. This paper examines the degree of robustness of selected integration algorithms by exposing them to a test profile of diverse problems. A number of common areas of difficulty are identified, and methods used to promote self sufficiency in these areas in the author's sparse matrix integration package GEARZ, are discussed. SUMMARY An extremely large number of algorithms have been proposed for the numerical solution of the ordinary differential equation initial value problem. The design or selection of such algorithms for use in a library of mathematical subroutines for general scientific and engineering use is difficult, as it is impractical to maintain a large number of subroutines which attempt to do the same job with varying degrees of success. Criteria for selection must include accuracy, efficiency and ease of use, thus one can consider realistically only those algorithms which have options to determine an optimal integration step size, and possibly order, by means of a built-in estimate of the associated truncation error, and are presented as quality software, complete with detailed internal and external documentation which emphasizes the weaknesses as well as the strengths of the method. These restrictions considerably ice the possibilities, but a number of candi- dates remain. Typical amongst these are algo- rithm i ,od by Gear [ 1 ] , Hindmarsh [2] , Byrne[3], ne(4) and Krogh[5J, all of which satisfy the Because these algorithms are reliable, quality software, numerical analysts have incorporated them in a multitude of applications packages, in which the user is shielded from the complexities involved in management of the integration. This type of application has in fact been their greatest success . To the uninitiated user of a subroutine library, however, their correct implementation can be seen as a prohibitive task. In order to focus further on quality for general use, the criterion of robustness must be added, a robust algorithm being defined as one which pro- duces the result to the desired accuracy or a clear indication of failure, requires a minimum of effort and does not demand clairvoyance from the user. This criterion is essential because the majority of projected users are not numerical differential equation experts, and therefore, not qualified to make decisions on fundamental issues such as the initial step size, error base, spar- sity, and error recovery. Unfortunately, as such decisions are frequently left to the user's dis- cretion in the guise of generality, the applica- tion of this robustness criterion eliminates most of the above-mentioned candidate algorithms. This paper attempts to show that such algorithms can be made considerably more robust at the ex- pense of a negligible loss of generality by fur- ther automation and simplification of the decision process during start up, integration, error pro- cessing, and discontinuity handling. Several specific areas where traditionally required user interaction can be either eliminated or reduced to optional status to improve effectiveness are discussed. They are: (a) Initial Step Size Selection: A routine with efficient step size adjustment is far better qualified than a user to provide this selec- tion. (b) Tolerable Error Imposition: An algorithm normally accepts a step if the estimated local truncation error in Y. satisfies € .v Lde - • ' |y| is ef- if this differs from the • iluation and Handling: Most irators require estimates of the in matrix from time to time, and the often left to speculate on a choice ;ow these should be obtained. Sophisti- fficiency by codinq their own .bian, but within a general package this na more remote, and could end it than numerical evaluation. ■ .1 nunerical procedure is to perturb Y to obtain the set 3F./3Y-, 7 n function evaluation calls. This • -zed; for example, a triJiaqonal matrix may be evaluated using only three calls, a:.d proportionate savings may be lized in .sparse matrices [6] . Sparse matrix Jacobian evaluation and handling optimized in this way is invaluable for eral PDE/ODE packages in which structure is arbitrary. ontimiit' -ion and Negotiation: A ntion definition occuring at an unknown tae dependent on critical rela- r.ships c .teqration is a ,-al phenomenon, but is traumatic • an integrator and frequently leads to •• limit or termination. A robust inte- jr should detect such a discontinuity and transcend it, issuing a suitable terse report. • .onal user interaction to facilitate smoother handling of discontinuities should be made available, for example see reference [7J. (e) Stiffness Detection: The most common way a standard non-stiff integration algorithm detects stiffness is to run to time limit. This is unsatisfactory as the results pro- duced will be of little value. A far less expensive approach is to report suspected stiffness ;4] and either refuse to waste e computing, or if possible, switch to a stiff algorithm. (f) Calling Sequence and Common Blocks: In a package, the user is shielded from direct communication through the integrator param- eter list, but this is mentioned here to facilitate library or package interface design. Any integrator calling sequence requires ten parameters, the derivative subroutine, dependent variable and derivative vectors and their size, independent variable and interval, tolerance, warning flag, one working storage array and current step size. these parameters are mandatory and all other communication may be through optional common blocks. The integrator itself can assign parts of the working storage array in predictable sequence, and choose a suitable method for handling the Jacobian from an assessment of the equation structure and available working storage. A number of results of both academii and ap] tions problems are given to illustrate ti modifications proposed above are not only more robust, but frequently more efficient. Finally, a brief description is given of the Gl package, which incorporates the modifications discussed into a combination of the Hindmarsh- algorithms [2] with the Curtis-Reid sparse matrix routines[6], and is currently used on the Atomic Energy of Canada Limited subroutine library and the F0RSIM(8], and MACKSIM[9] simulation pack, REFERENCES 1. C.W. Gear, The Automatic Integration of Ordinary Differential Equations, Comm. V14, 3, p. 170, March 1971 . 2. A.C. Hindmarsh, The LLL Family of Ordina Differential Equation Solvers, UCRL'-78129, April 1976. 3. G. Byrne and A.C. Hindmarsh, A Polyalgonthm for the Numerical Solution of Ordinary dif- ferential Equations, ACM TOMS, VI, p. 71, 1975. 4. L.F. Shampine, H.W. Watts and S.M. Davenpi d , Solving Non-stiff Ordinary Differential Equations - The State of the Art, Siam Re- view, V18, 3, 1976. 5. F.T. Krogh, Variable Order Integrators for the Numerical Solution of Ordinary Differ- ential Equations, Jet Propulsion Lab., Sec- tion 314 Subroutine Write-up, May 1969. 6. A.R. Curtis and J.K. Reid, FORTRAN Subroutines for the Solution of Sparse Sets of Linear Equations, AERE-R-6844 , 1971. 7. M.B. Carver and S.R. MacEwen, Analysis of a System Described by Implicitly Defined Ordi- nary Differential Equations Containing Dis- continuities, Applied Mathematical Modelling, V2, 1978, p. 280. 8. M.B. Carver, D.G. Stewart, J.M. Blair and W.N. Selander, The FORSIM VI Package for Automated Solution of Arbitrarily Defined PDE/ODE Systems, Atomic Energy of Canada Limited report AECL-5821, February 1978. 9. M.B. Carver and A.W. Boyd, A Program Package Using Stiff Sparse Techniques for the Auto- matic Solution of Mass Action Chemical Kinetics, submitted to Int. J. Chem. Kinet. 10. M.B. Carver and A. P. Baudouin, Solution of Reactor Kinetics Problems Using Sparse Matrix Techniques in an ODE Integrator, Atomic Energy of Canada Limited report AECL-5177, 1975. 36-2 ON AN 1 1T.R.N THE APPRO [ION OF SOME ORDINARY DIFFERE1 I [ONS P. E. Zadunaiskv rio Naciona] de Fiaica Cosmica, San Miguel, Argentina) G. Lafferriere tciaa Exactas, L'niversidad de Buenos Aires, Argem ABSTRACT Let us consider a system of ordinary differential equations of the form ■ , v ' v ) = where y and f are vector functions of x. By introducing an operator T such that Tu = f(x, u, u',..., u ) we have Ty = 0. Assuming that y is an approximation of the solution y(x) a generalisation of Newton's method can be applied to improve, under certain conditions, such approximation by the recursive algorithm y i+1 = y - I*(y ) "Ty 1 (i = 0, 1, -,...) and assuming that T'(y ) exists. In the present case we can use such an approach in a numerical fashion as follows. After obtaining by any method of integration numerical approximations y on a discrete set of points x (n = 1, . . . , N) we interpolate them by a convenient function P(x). By taking this interpolant as the first analytical approximation y Newton's process is applied pointwise in order to correct by iterations the discrete approximations y . This procedure may become rapidly convergent especially in some stiff problems where we have obtained so far promising results. 37-1 A Portable Automatic Taylor Series (ATS) Compiler for Solving Ordinary Differential Equations* by Y.F. Chang, Roy Morris+, Computer Science Department, University of Nebraska Lincoln, Nebraska 68588 and George Corliss, Department of Mathematics and Statistics, Marquette University Milwaukee, Wisconsin 53233 EXTENDED ABSTRACT A portable compiler computer program is de- veloped for the automatic numerical solutions of ordinary differential equations. The inputs to this compiler are FORTRAN statements of the equa- tions (differential and algebraic) and of the initial conditions. The output of this compiler is an object program written in FORTRAN. The Automatic Taylor Series (ATS) method of solving ordinary differential equations is based on the premise that differentiation of functions can be performed at specified points by simple additions and multiplications. The solution of very com- plicated ordinary differential equation can be expanded into a long Taylor series with ease. The availability of very long Taylor series has led to new understanding of the convergence pro- perty of the series in relation to the locations of singularities and to the radius of convergence. It is then possible to set an error limit before- hand and maintain local error control over the computations within that limit up to almost machine accuracy. The power of the Automatic Taylor Series (ATS) method lies in the speed and accuracy that is obtainable. The foundation of the Automatic Taylor Series (ATS) method is the replacement of diff- erentiation of functions by simple arithmetic operations. The first publication on this simplification was by Wilsonfl] in 1949. Gibbons[2] obtained the recursive relations for the fully automatic differentiation of some simple fucntions in 1960. A summary of power series method was published by Leavitt[3] in 1966. Application of power series to the auto- matic solutions of ordinary differential equa- tions was first discussed by Moore[4] in 1966. Barton, et al [ 5] in 1970 first reported imple- mentation of an automatic solution. That auto- matic solution, by Barton, et al , was obtained through a compiler program that was written in a machine-dependent language. One of the pre- sent authors published an account of a compiler program written in FORTRAN for the automatic solutions of ordinary differential equations[6] * This work is in part supported by N.S.F. Grant ' j 78-03022. ♦ Present address: Department of Computer Science, Iowa State University, Ames, Iowa 50011 in 1974. The principal differences between this and the one by Barton, et al are the calculation for the series' radius of convergence and the concomtitant error control. This compiler was written in FORTRAN with hope that it would be portable. This proved to be optimistic, not factual. There are, at present, two versions of this compiler: one for the CDC-computer, and one for the IBM-computer. These two versions are quite distinct. The present paper is the report of a por- table compiler program for the automatic solu- tions of ordinary differential equations. This program is written in FORTRAN with special treat- ment of machine-dependent parameters to obtain portability. Parameters such as bits/character, character/word, single-precision limit, and double-precision limit are stored in a single Blockdata location for easy access and quick change. Since the bits/character and word length are machine-dependent, we have also Instituted a numerical code for character strings. This allows for easy character string manipulations us inn integer arithmetic; thus, we have overcome the difficulty of character string manipulation on certain computers. As development proceeds, testing is being done on computers from several manufacturers: IBM, CDC, PDP, and Xerox Sigma. The input to this compiler program consists of four blocks of user specified information: the statement of differential equations and three blocks of instructions to be inserted into the object program. Each of the data blocks ends with a data terminator ($), and only one data terminator may appear on a line. The differential equations are entered in the first block, with (DIFF) as the differential operator. For example, DIFF(Y,X,2) denotes y"(x). The rest of the ex- pressions are entered in typical FORTRAN state- ments. For example, to solve the problem: y" = y' sin(y) exp (x); y(0) = 1, y'(0) = -1 on [0,4] the input would be 38-1 Column 7 DIFF(Y,X,2) ■ DV*SIN(Y)*EXP(X) DY ■ DIFF (Y.X.I) $ START =0.0 4.0 1.0 - 1.0 S END Y (1) Y (2) S This is an example of the simplest possible input. The use of the first and third data blocks is man- datory, while the use of the second, and fourth data blocks is optional. The sophisticated user may control the solution of very complicated problems through these optional data blocks. The second data block may be used to enter non-executable FORTRAN statements, such as COMMON, DATA, or SUBROUTINE, into the object program at the appropriate location. Using the SUBROUTINE statement allows the user to interlace many com- plicated programs together. The third data block is used to define the solution interval and initial conditions. The fourth data block may be used for such things as changing the amount of printout, starting and stopping the solution at any given point in either the dependent or inde- pendent variables, and the calling of other sub- routines. Problems with piecewise analytic solu- tions can be solved with ease. This compiler program reads the input and generates a FORTRAN object program, which sub- sequently is compiled and run to solve the speci- fied problem. The object FORTRAN program gen- erated by this compiler conforms as closely as is possible to the general structure suggested by Hull and Enright[7]. As the first data block of equations are read in, this compiler produces a list of tokens. Each token represents a con- stituent in the equations. Variables and con- stants are stored in arrays as integers with the corresponding tokens pointing to the proper element in the arrays. The token produced for a function is an integer code representing that function. The list of tokens are checked for gramma- tical correctness and put in Polish postfix form. The grammar used is an LL(1) grammar, which has the property that the top stack symbol and the input token uniquely determine which production of the grammar should be applied next. The postfix expression is then arranged in triplets of operator-operand-operand, from which the final object program is generated. This compiler is a significant improvement over the earlier compilers for the automatic solution of ordinary differential equations. Not only is it portable between different com- puters, by the simple change of an easily access- ible blockdata, there are other important features. The calculation for the series' radius of con- vergence has been improved such that the order of the singularity in the complex plane of the solu- tion function can be estimated with a certain degree of accuracy. We have also included in this compiler the ability to vary the ratio of step- size to the radius of convergence and variable order in the numerical integration. Other improvement include the choice of single- or double-precision, and comprehensive diagnostic error messages . Copies of this compiler will be made avail- able to researchers in tlie field of numerical ordinary differential equations upon reouest, when all of the major bugs have been removed and full documentation is ready. REFERENCES (1) E.M. Wilson, A Note on the Numerical Inte- gration of Differential Equations, Quar. J. Mech. & Appl. Math., Vol. 2 (1949), p. 208. (2) A. Gibbons, A Program for the Automatic Integration of Differential Equations Using the Method of Taylor Series, Comp. J., Vol. 3 (1960), p. 108. (3) J. A. Leavitt, Methods and Applications of Power Series, Math. Comp., Vol. 20 (1969), p. 46. (4) R.E. Moore, Interval Analysis, Prentice- Hall, Englewood Cliffs, New Jersey, 1966. (5) D. Barton, I.M. Willers, and R.V.M. Zahar, The Automatic Solution of Systems of Ordin- ary Differential Equations by the Method of Taylor Series, Comp. J., Vol. 14 (1970), p. 243. (6) Y.F. Chang, Automatic Solution of Differ- ential Equations, Lecture Notes in Math. #430, Constructive & Computational Methods for Differential and Integral Equations, D.L. Colton and R.P. Gilbert, ed. , Springer- Verlag, New York, 1974, p. 61. (7) T.E. Hull and W.H. Enright, A Structure for Programs that Solve Ordinary Differen- tial Equations, Dept. of Computer Science Technical Report #66, Univ. of Toronto, Toronto, Canada (1974). 38-2 ROSENBROCK-TYPE METHODS FOR THE NUMERICAL SOLUTION OF STIFF SYSTEMS Peter Kaps University of Innsbruck Abstract : Generalized Runge-Kutta methods for the numerical integration of stiff systems are investigated. A simple way of obtaining the order conditions is given. Numerical results of two embedded (3)4 methods (GRK4A, A-stable, and GRK4T, A(89°) -stable with small truncation errors) are presented. Comparing computing time both methods are competitive with the GEAR program. In this paper I want to give a survey of the results of my thesis [ 1 ] and two articles I wrote together with P. Rentrop [2] and G. Wanner [3]. For the numerical solution of the initial value problem y' (x) = f(y(x)), y(x Q ) = y Q defined in an n-dimensional real or complex space the following method is considered: E = I-Yhf ' (y Q ) Ek 1 = hf(y Q ) Ek 2 = hf(y Q +a 21 k 1 ) + hf(y )Y 21 k l Ek s = hf 'y + giV--- +a s,s-i k s-i ) + hf, ) = f c J e *p (x j (x - x o»-fyr if x j * Cj + b\(x-x Q ] if A. = 0. v J ' is the eigenvector corresponding to A., the d. are constants determined by b and A, and the c- J J are constant determined by the initial vector y . In our problem we assume that some of the initial values (components of y ) are known and others are unknown. We wish to specify these unknown values in such a way that rapid, short-term fluctuations in the solution are minimized. These short-term fluctuations are caused by those terms of the form cexp(x.(x-x )) for which A. has large negative real part. We can eliminate these terms by specifying the initial values in such a way that c- = 0, j = k+1 ..... ,n . A solution having this property is said to be in stiff equilibrium . Not surprisingly, in order to obtain a solution which is in stiff equilibrium, the number of initial values which we allow to vary must be equal to be number of stiff eigenvalues. This restricts the class of problems which can be solved, but the restriction is not as severe as it might at first seem because there is usually some latitude in the choice of stiff eigenvalues.. For example, if A has eigenvalues 1, 0, -1,-10 ,-10 , 40-1 -10^, we may choose to designate one, two, or three eigenvalues as stiff. By a straightforward procedure we can deter- mine the initial values which give a solution in stiff equilibrium. First, yV(x ) = c.-y.b., _ 1 J o J J J j=l,....,n, where u, = \~ if A, / and u. = if • • = 0. These equations can be expressed in matrix form as y*(x ) = c - Mb\ where M is the diagonal matrix whose jth main diagonal entry is p.. Thus y = V(c-mId), where V is the matrix 3 ° (i) whose jth column is v . This matrix equation can be written in block form (1) where 11 '21 V 12 '22 d-r e-s z consists of those components of y which are unknown, e consists of those components of c which are to be set to zero, and r and s are components of the known vector Mb. For any specified w we wish to find z such that this o system has a unique solution c wi th e = 0. Such a z must satisfy the equations V n (d-r) V 2] (d-r) V ]2 s V 22 s obtained from (1) by setting e = 0. These equa- tions have a unique solution if and only if V,, is nonsingular. If V,, is nonsingular we can solve the first equation uniquely for d-r. We can then substitute d-r into the second equation to calculate z . The procedure which we have just outlined in- cludes the calculation of the eigenvalues and eigenvectors of A and the solution of two linear systems. This would be costly for a large system of differential equations, but for a system of modest size (one of ten equations, say) these operations are not expensive. Now consider the nonlinear problem y' = f(y)> y( x n ) = y > where f is assumed twice continuously differentiable. We are concerned here with the behavior of the solution in the vicinity of y. Near y the solution to the non- linear problem is well approximated by the solution to the linear problem y" = f(y o } + fy (y o )(y - y o ) ' y(x o> = y o • Ay+b, where As in the This problem has the form y' A = '¥- (y ) and b = f(y ] 3f / v ■z- (y )y ;jy u o' ; o linear case, y contains some unknown components, so we cannot simply proceed as in the linear case. Instead we must use an iterative procedure. We make an initial guess of y , use this initial guess to calculate A and b,°and proceed as in the linear case to calculate a new guess of y . We repeat this procedure until (hopefully) it con- verges. In the nonlinear case the partitioning of the eigenvalues is complicated by the fact that the matrix A changes from iteration to iteration, and the eigenvalues change with it. Thus, a partition which is reasonable on one iteration may not be reasonable on the next. It is im- practical to change the number of initial values which are being adjusted from one iteration to the next, so we must decide at the beginning how many we wish to vary and stick with that number. If we decide to let i initial values vary, then at each iteration we place the i "stiffest" eigenvalues in the "stiff" set and all others in the "small" set. Once the iterations have converged to some y , we can check whether the partition is reasonable by examining the eigen- 3f values of — (y J . We can also monitor the ay partition at intermediate iterations, but this is not necessary. All that really matters is whether the final partition is reasonable. NUMERICAL RESULTS We present here some numerical results for the system of transport equations discussed in the introduction. Consider the following initial conditions: altitude = 1500 km, temperature = 2000 K, heat flow = -1.943 x 10 -5 erg(cm)-2 (sec) . The unknown initial condition is stress. The following table shows the stress values calculated by the method proposed here. Iteration 1 2 3 Stress (K) 0.0 -5.159232 -5.163496 -5.163496 We see that the procedure converges to seven decimal places in two iterations. At each iteration the Jacobian matrix was estimated by numerical differencing. The eigenvalues and eigenvectors were calculated by the IMSL sub- routine EIGRF, which uses a shifted QR algo- rithm. The initial values have been partitioned into three known values (counting altitude, the independent variable) and one unknown value. In order for this to be a reasonable partition we must be able to partition the Jacobian into a set of three "small" eigenvalues and a set of one "stiff" eigenvalue. The following table shows the eigenvalues of the Jacobian computed on the final iteration. eigenvalues 7.18 x io" 10 0.0 g -4.47 x 10 * -7.80 x io" 6 scale heights 2.24 x io cm = 2240 km 1.28 x lo 5 cm = 1.28 km If we designate the first three eigenvalues "small" and the fourth one "stiff", then the "stiff" eigenvalue is over 1000 times larger in absolute value than any of the "small" ones. Thus the partition is reasonable. The fourth eigenvalue shows that the system 40-2 is quite stiff. It looks small but it is actually large, given the units of the problem. It; the solution the term which corresponds to eigenvalue \ decays by a factor e in the distance -■' . This distance is commonly called the scale f the term. We are using cgs units, so ale height is in centimeters. If we express tnis distance in kilometers, then the scale height of the term corresponding to the fourth eigenvalue is 1.28 km. This is an extremely short ?nce, given that we are integrating out some ten thousand kilometers. By contrast, the scale tit associated with the next largest eigenvalue er 2000 • Finally, we have to ask whether the answer en by tnis procedure is meaningful. Is 63496 really the "right" initial stress value? Yes, it is. If we integrate outward using this initial value we get smooth profiles wnich do not show any short term fluctuations. The use of any other initial stress value leads to Iarqe initial fluctuations in the profiles are unrealistic from a physical standpoint. 40-3 A COMPARISON OF MKTHODS FOR COMPUTING STURM- LIOUVILLE EIGENVALUES J.W. PAINE Australian National University An important source of Initial value problems arises from the use of shooting methods to solve Sturm- Li ouvi lie eigenvalue problems of the form x e (a,b) (l.a) - (pu) ' + qu = ' ru (l.b) a u(a) + Bu(a) = y u(b) + 6u(b) = (1. c) u, pu e C[a,b] Because of the large variety of methods available for approximating the eigenvalues of this problem, and because of their different underlying assumptions on which the methods are based, it is important to know and to understand the different properties of the eigenvalue approximations which will be obtained. The most important property that any method should possess is convergence, and (if we assume that each method is parameterized over the set of positive integers) in it's simplest form this is expressed as Um X (N). (N) where X is the <'th exact eigenvalue and X < " K is the K'th eigenvalue of the N'th approximat- ing problem. But although such basic converg- ence results are necessary, they are not altogether useful since they do not really give much information about the approximations that will be obtained for a fixed finite value of N. As a result of this, the approximation results for most methods are refined, and given in the form IX - X (N) | < C N" P ; N > N (k) where C generally depends on < but is independ- ent of N. For many situations such estimates are sufficient to guide one's choice of method. However since these results are local convergence results in the sense that for a fixed value of k, we can find an N large enough such that the convergence estimate is valid for every . , there are situations when the information given by such local convergence results will be insufficient or even misleading, when comparing different methods. One such situation is when there is a large number of eigenvalues to be approximated and there is a definite practical need to keep the size of N as small as possible. An example of this type of problem arises in the use of spherically symmetric models of the earth to obtain information on the structure of the mantle from known values of the torsional free oscillation frequencies; i.e. the inverse eigenvalue problem. One approach for solving this problem is to choose a set of possible structures and by finding how well the eigen- frequencies of each of these models match the desired values, the "best possible" model (or models) from the given set can be chosen. The difficulty is that the number of eigenvalues over the whole set of possible structures becomes quite large and so the time and computing resour- ces spent in finding the model's eigenfrequencies needs to be optimised in terms of the information obtained. What is needed to improve the error bounds is to give them a more global character, and one possible means for extending the bounds to show this is to refine the bound so that the dependence of C on k becomes explicit. Thus for the purposes of comparison, we adopt (N) X' < C N F X "Pi«l N > N (K) o v ' (where C is now independent of both k and N) as the standard error bound. If we consider the bounds obtained in this form for a number of representative methods (Rayleigh-Ritz, Prufer substitution, modified Prufer substitution, finite difference and a similar operator approach) a number of features of the methods and the resulting approximations become apparent . The first point is that since the rate of growth of the error of many methods is directly related to some high order derivative of the eigenfunction (the larger the order of convergence p, the higher the order of derivative) it can be preferable to use a method with a lower order of convergence to obtain better approximations of the larger eigenvalues at the expense of losing the high accuracy approximations 41-1 of the smaller eigenvalues. The second point is that some methods (such as modified Prufer substitution) perform better (i.e. the error behaves more uniformly) if the original eigenvalue problem is restated in Liou- ville normal form (2. a) (2.b) Av t € (0,1) a*v(0) + 6v(0) = Y*v(l) + 6v(l) = (2.c) v, v e C[0,1] The final point is that for the considered approximation methods, it generally occurs that the improved uniformity of the approximation (if any is possible) is obtained at the expense of requiring stronger regularity conditions on the problem (e.g. the coefficients p and r should be sufficiently highly differentiable so that the problem can be transformed to normal form (2)). Thus it is desirable to find which methods will give the best results for a given set of regularity conditions. Although such an extens- ive analysis has not been undertaken, particular emphasis has been on the similar operator approach (based on transforming (1) to a quasi- normal form and then replacing the single coefficient function by a piece-wise constant approximation, the resulting problem is then analytically solvable) which gives uniform approximations under very weak regularity conditions. m-2 THE NUMERICAL SOLUTION OF SECOND ORDER SYSTEMS OF ODES Organizer: W. H. Enright Dept. of Computer Science University of Toronto Panelists: D. G. Bettis Dept. of Aerospace Engineering and Engineering Mechanics University of Texas at Austin P. Deuflhard Institut Fur Angewandte Mathematik Universitat Heidelberg F. T. Krogh Jet Propulsion Laboratory California Institute of Technology The numerical solution of second order systems is one of the most important application areas for initial value methods. Special numerical methods for second order equations have been in use for over fifty years and these methods continue to be improved and developed in parallel with methods for first order systems. Two areas where these problems frequently arise and where special methods are heavily used are celestial mechanics and structural dynamics. Consequently the analy- sis and development of new methods is often tailored to classes of problems from these areas. Because of the gain in efficiency that is possi- ble when the differential equation is inde- pendent of y' f there are generally two classes of methods developed for second order equations. The first class of methods are suitable for systems of the form: y"=f(x,y). y(x )=y , y'(x )=y' ; (I) and the second class is suitable for general systems: y"=f (x.v.y'), y(x n )=y., y'(x.)=y" (ID The members of this panel have each contri- buted to the improvement and development of numerical methods for second order problems. Professor Bettis has recently developed (with K. Horn) efficient embedded Runge-Kutta-Nystrom methods for second order systems (Bettis and Horn (1978)). He has also recently been involved in the development of exponentially fitted Runge- Kutta formula pairs. Professor Deuflhard has had considerable experience with extrapolation methods for seccnd order equations and has recently de- veloped a class of exponentially fitted extrapo- lation methods (Deuflhard (1978)). Dr. Krogh is well known for his work in the development of efficient multistep methods for systems of ODES. The multir::ep methods he has developed can handle second order equations directly and he has in- vestigated the advantages of doing so (Krogh (1975)). This discussion is to be an informal discussion of the state-of-the-art in this area. Each panelist will discuss aspects of this area that he has been Involved in or finds to be of particular interest. It is of course impossible to anticipate all the issues that will be addressed, but we will list a few of the issues and questions that are currently under active investigation and which are likely to arise. The first issue is the obvious one of detailing the advantages direct methods have over the alternative approach of reducing the system to a larger system of first order equations. Of course this is very much related to the general topic of comparing numeri- cal methods (to be discussed in another panel session). The difficulties involved in comparing methods are particularly relevant here since we are not only interested in comparing different direct methods, but also in quantifying the ad- vantage of direct methods. A second topic which has recently received attention from several investigators (Oear (1978), Dahlquist (1978) and Thomas and Gladwell (1978)), is the question of stability of direct methods. Are there classes of problems where a large region of absolute stability is important, and if so, what is a suitable definition of sta- bility for direct methods? This issue seems particularly crucial to problems arising in structural dynamics. Another issue that is of considerable interest is the question of what an appropriate error control strategy should be for second order equations. Although this question is reasonably well understood for first order systems, the situation for second order systems is not as clear. Is it necessary to monitor and control the local error in y' as well as y for problems '13-1 vpe I (or II). Some theoretical as well as computational Investigations would be helpful here. As a final example of issues that will likely be discussed, we will mention exponential fitting. This idea of choosing a parameter (or parameters) of the formula to match certain known character- istics of the solution is a topic that has received much attention in the literature for both first and second order equations. It is important to realize, that for this approach to be viable, one must take care in the selection and implementation of compatible error estimators and stepsize strategies. These components of a method must recognize the 'built-in' accuracy of the underlying formula and also account for the fact that the 'error coefficient' of the formula will now depend on the stepsize. It is also important to consider whether these formulas will general- ize to systems of equations, and whether this will require an approximation to the matrix expo- nential and the derivation of formulas with matrix coefficients. References: D.G. Bettis and M.K. Horn (1978) Embedded Runge- Kutta-Nystrom algorithms of order two through six, T1C0M Report 78-3, Texas Institute for Computational Mechanics, University of Texas at Austin. G. Dahlquist (1978) On accuracy and unconditional stability of linear multistep methods for second order differential equations, BIT 18, pp. 133-136. P. Deuflhard (1979) A study of extrapolation methods based on multistep schemes without parasitic solutions, ZAMF, to appear. C.W. Gear (1978) The stability of numerical methods for second order ordinary differ- ential equations, SINUM 15, pp. 188-197. F.T. Krogh (1975) Summary of test results with variants of a variable order Adams method, Computing Memorandum No. 376, Jet Propulsion Laboratory, Pasadena, California. R. Thomas and I. Gladwell (1978) Stability of some methods for second order ordinary differential equations with engineering applications, Numerical Analysis Report No. 29, Dept. of Mathematics, University of Manchester. 43-2 PANEL DISCUSSION ON NUMERICAL METHOD OF LINES SOFTWARE G. D. Byrne University of Pittsburgh A. C. Hindmarsh Lawrence Livermore Laboratory M. B. Carver Atomic Energy of Canada Limited M. Minkoff Argonne National Laboratory G. K. Leaf Argonne National Laboratory W. E. Schiesser Lehigh University R. F. Sincovec Boeing Computer Services, Inc. 1. Introduction (G. D. Byrne) The numerical method of lines (NMOL) involves the following. A system of partial differential equations (PDE's) is discretized in its space-like variables to obtain a system of ordinary differ- ential equations (ODE's) in the time-like variable, which is then solved by an ODE package. In a numerical method of lines package, the spatial discretization may be by a predetermined finite difference scheme, a user specified finite difference scheme of a certain general type, a Galerkin procedure in conjunction with B-splines, or a collocation procedure in conjunction with B-splines. In the case of Galerkin and collocation procedures, essential boundary conditions lead to algebraic equations or constraints, which are interspersed with the ODE's which must be solved. For the very simple heat equation, u = u , finite difference-NMOL t xx involves the solution of a system of ODE's of the form dy/dt = f(y,t), while collocation-B-spline or Galerkin-B-spline NMOL involves A- (dy/dt) = g(y,t) and the coefficient matrix A may be singular. Systems of these types of ODE's also arise from less simple PDE's for these discretizations. Of course, other discretizations are possible. The panelists are associated with several general purpose NMOL software packages which use each of the discretization procedures cited above. Several of the panelists have also developed special purpose NMOL software for the solution of specific problems in heat conduction, oil reservoir simulation, and combus- tion. Some panelists have also been involved in the development of ODE software and in work with various nonlinear and linear systems solvers for use within the ODE software as in an NMOL package. In the abstracts which follow, the partici- pants are concerned with the requirements imposed on the ODE solver by the types of ODE's solved in NMOL packages, desirable features of the ODE solver, the interface between the ODE solver and the rest of the NMOL code, the interface between the user and the NMOL code, and possible directions of research in packages for large, stiff systems of ODE's, which arise in the NMOL solution of multidimensional systems of PDE's. 2. The Ordinary Differential Equation Algorithm as Core Component of a Method of Lines Partial Differential Equation Software Package (M. B. Carver) The properties required of an ODE integrator within a method of lines PDE package are different from those required of an integrator offered in a general subroutine library only in one respect: that the user is one more step removed from hands- on control of the integrator. Even in a subroutine library, the average user lacks sufficiently specialized knowledge to guide integration by parameter selection. In either context, therefore, one requires sufficient robustness to automatically negotiate predictable difficulties without the guidance of a clairvoyant user, but enough scope for the informed user to provide tuning adjustments. A number of efficient, well -documented, error-controlled, variable-step-size, integration algorithms are available, for example those discussed in references [1 to 4] and it is a prerequisite to provide such a routine as the core of any integration package. However, there are several specific areas where traditionally required user interaction can be either eliminated or reduced to optional status to improve effectiveness. They are: (a) Initial Step Size Selection: A routine with efficient step size adjustment is far better qualified than a user to provide this selection. W-l (b) Tolerable Error Imposition: An algorithm normally accepts a step if the estimated local truncation error in Y. satisfies T0L*Ybase., where Ybase. is frequently defined as the maximum attained value of Y. or made available as an array for user speci- fication. Again users have no criteria on which to nominate each Ybase., and the automatic assignment Ybase- ■ max( | Y. | ,Ysig) J J is general, and requires the user only to provide Ysig, the value below which any jY| is effectively zero, if this differs from the machine unit round off. (c) Jacobian Evaluation and Handling: Most general integrators require estimates of the Jacobian matrix from time to time, and the user is often left to speculate on a choice of how these should be obtained. Sophisti- cates may gain efficiency by coding their own Jacobian, but within a method of lines package this possibility grows more remote, and could end up less efficient than numerical evaluation. Standard numerical procedure is to perturb each individual Y. to obtain the set 8F./3Y., requiring n function evaluation calls. This may be optimized; for example, a tridiagonal matrix Jacobian may be eval- uated using only three calls, and propor- tionate savings may be realized in sparse matrices [5]. Sparse matrix Jacobian evaluation and handling optimized in this way is invaluable for general PDE/ODE packages in which structure is arbitrary. The mechanics are further discussed in a paper to this conference [6]. (d) Discontinuity Detection and Negotiation: A switch in equation definition occuring at an unknown time dependent on critical relation- ships evolving during integration is a common physical phenomenon, but is traumatic for an integrator and frequently leads to time limit or termination. A robust integrator should detect such a discontinuity and transcend it issuing a suitable terse report. Optional user interaction to facilitate smoother handling of discontinuities should be made available, for example see reference [7]. (e) Calling Sequence and Common Blocks: In a package, the user is shielded from direct communication through the integrator parameter list, but this is mentioned here to facilitate library or package interface design. Any integrator calling sequence requires ten parameters, the derivative subroutine, dependent variable and derivative vectors and their size, independent variable and interval, tolerance, warning flag, one working storage array and current step size. Only these parameters are mandatory and all other communication may be through optional common blocks. The integrator itself can assign parts of the working storage array in predictable sequence, and choose a suitable method for handling the Jacobian from an assessment of the equation structure and available working storage. The GEARZ integration algorithm used in the F0RSIM[8] and MAKSIM[9] simulation packages combines the above consideration with the work described in references [2, 5 and 10]. References [1] C. W. Gear, The Automatic Integration of Ordinary Differential Equations, Comm. ACM, V14, 3, p. 176, March 1971. [2] A. C. Hindmarsh, The LLL Family of Ordinary Differential Equation Solvers, UCRL-78129, April 1976. [3] G. D. Byrne and A. C. Hindmarsh, A Poly- algorithm for the Numerical Solution of Ordinary Differential Equations, ACM TOMS, VI, p. 71, 1975. [4] L. F. Shampine, Practical Solution of Ordinary Differential Equations by Runge- Kutta Methods, SAND760585, 1976. [5] A. R. Curtis and J. K. Reid, FORTRAN Routines for the Solution of Sparse Sets of Linear Equations, AERE-R-6844, 1971. [6] M. B. Carver, In Search of a Robust Integra- tion Algorithm, paper submitted to ACM SIGNUM Conference on Numerical Ordinary Differential Equations, 1979. [7] M. B. Carver and S. R. MacEwen, Analysis of a System Described by Implicitly Defined Ordinary Differential Equations Containing Discontinuities, Applied Mathematical Modelling, V2, 1978, p. 280. [8] M. B. Carver, D. G. Stewart, J. M. Blair and W. N. Selander, The FORSIM VI Package for Automated Solution of Arbitrarily Defined PDE/ODE Systems, Atomic Energy of Canada Limited report AECL-5821 , February 1978. [9] M. B. Carver and A. W. Boyd, A Program Package Using Stiff Sparse Techniques for the Automatic Solution of Mass Action Chemical Kinetics, submitted to Int. J. Chem. Kinet. [10] M. B. Carver and A. P. Baudouin, Solution of Reactor Kinetics Problems Using Sparse Matrix Techniques in an ODE Integrator, Atomic Energy of Canada Limited report AECL-5177, 1975. 3. Method-of -Lines Considerations in ODE Software Design (A. C. Hindmarsh) Among the most frequent uses of ordinary differential equation (ODE) software are those for the solution of partial differential equations (PDE's) by the numerical method of lines (NMOL). Typical of such problems are system^ of parabolic PDE's in time t and a space vector x, given initial values at t = and boundary conditions in W-2 x. The resulting ODE problems are among the most challenging because of their potential size, complexity, and stiffness. PDE-based problems are frequently (but not always) stiff, either because of physical processes present (e.g. chemical kinetics), or because of the particualr features of the spatial discretization (e.g. a simple 1-D diffusion equation finite-differenced on a fine mesh). The combination of size and stiffness imposes severe restrictions on the choice of ODE methods and on their implementation — restric- tions that differ greatly from those for nonstiff problems. Without going into the details of NMOL discretizations, the final result can be given in one of two forms. For finite difference approaches, one generally obtains an explicit ODE system y = f(y,t), where y is a vector of length N, and an initial value vector y is known. For weighted residual methods (collocation, Galerkin, fintie element, etc.),_ one gets instead a linearly implicit ODE system Ay = g(y,t) where A = A(y,t) is an N x N matrix (usually nonsingular) wnose form depends on both the discretization method and the geometry of the problem. In the latter case, cost considerations require that the ODE algorithm treat the system in a direct manner, rather than A-V invoke an algorithm for y = f with f When the ODE is stiff, accuracy and numerical stability dictate the use of implicit methods. Nearly all implicit methods, when implemented so as to retain the desired numerical stability, solve linear algebraic system that involve the N x N matrix J = 3f/3y, or both A and J = 3g/3y in the case Ay = g. While the literature now contains numerous stiff ODE methods, most NMOL applications use methods based on the Backward Differentiation Formulas (BDF's). The reasons are partly historical, and partly based on the following facts about BDF methods: a modified Newton (chord) iteration on the corrector, but handle only y = f and assume J to be dense. They save LU factorizations of the required matrix, and recompute this only as needed. In NMOL, J and A invariably have much sparse structure, and property 4 becomes crucial. To capitalize on it, variant versions of GEAR and EPISODE have been written, dictated by structures that commonly occur in PDE-based ODE systems. The modular design of the solvers easily allows the dense linear system treatment to be replaced accordingly. The implicit ODE case requires some alterations to other parts of the algorithm as well. The following variants exist currently: 1. GEARB treats J in band form. It is used in the NMOL packages PDEPACK, M0L1D, PDETWO, DSS/2, and FORSIM. 2. EPISODEB is the band variant of EPISODE, analogous to GEARB. 3. GEARBI treats J in blocked form and uses a block-iterative linear system algorithm (block-SOR). It is used in various 1-D and 2-D applications. 4. GEARS assumes a general sparse form for J and uses the Yale Sparse Matrix Package. 5. GEARIB treats the implicit ODE Ay = g with A and J = 3g/3y in band form. It is used in the NMOL packages PDECOL and DISPL. 6. EPISODEIB is the variant of EPISODE analogous to GEARIB. In the explicit ODE case, these solvers retain a nonstiff (Adams) option for convenience, and some of the NMOL packages relay this option back to the user. 1. They form a group of stiffly stable methods of various orders, easily allowing for variability of step size and order. 2. They require the solution of only one N x N nonlinear algebraic system on each time step, and so require only one copy of J or a related matrix (involving A, if present ) at a time. 3. They do not need J (or A, if present) to be very accurate, and so do not require a reevaluation of J at every step. 4. They easily allow for sparse structure of J (and A, if present) to be utilized to gain both speed and storage economy. Few other methods have all of these properties. There are two ODE solvers which have seen particularly heavy use in the context of NMOL: The GEAR package uses (for stiff problems) fixed- step BDF' r > with interpolator step changing. The newer EPISODE package resembles GEAR externally but is based on truly variable-step BDF's. Both advantage of properties 1 to 3 above, using Some applications have involved more special variants. For example, GEAR has been combined with a tailored preconditioned conjugate gradient algorithm for use on radiation diffusion- scattering problems. - Alternatives to chord iteration exist but have not been tried extensively in stiff ODE software. Experiments with nonlinear block-SOR iteration in an NMOL context have shown considerable promise. Another promising approach, in 2-D and 3-D, is to combine a BDF (or similar) method with a multiplicative matrix splitting of the type used in ADI. 4. Some Comments on the Impact of PDE Software on ODE Solvers (G. K. Leaf and M. Minkoff) A common feature of software for solving parabolic differential equations is the use of the method of lines for the reduction of the partial differential equation to a system of ordinary differential equations. A typical PDE package will then use an ODE package to solve the resulting system of differential equations. The use of ODE solvers in conjunction with PDE software imposes certain requirements on the ODE W-3 solver. First we note that the method of lines usually leads to systems of stiff ordinary differential equations; hence implicit time differ- encing is generally required. Such implicit method use Newton-type nonlinear equation solvers which require estimates for the Jacobian. Particularly in two spatial dimensions, the storage of this Jacobian becomes the dominant consideration in the implementation of PDE software. With the use of local approximations for spatial effects, the storage problem is ameliorated somewhat by the banded structure of the Jacobian. When a Galerkin type procedure is used, the resulting system of ODE's will be implicit in time. This makes it desirable to have solvers which can treat implicit time dependence in an efficient manner. When the method of lines is used, the boundary conditions associated with the PDE's lead to system of algebraic-differential systems. Finally, with PDE's a substantial cost savings can be achieved with the use of dump/restart features. Such features are almost mandatory in ODE solvers used in conjunction with PDE software. 5. Two Proposed Standard Interfaces to Facilitate the Use of NMOL Codes (W. E. Schiesser) The remarkable flexibility and versatility of the numerical method of lines (NMOL) will most certainly facilitate the integration of systems of ordinary and partial differential equations over a spectrum of application areas. This is particularly true for scientists and engineers who do not wish to become experts in numerical analysis. If some general guidelines are developed for interfacing various components of a NMOL code, the user will be able to take better advantage of the flexibility of the method, and the exchange of applications between users will be facilitated. In this presentation, two interfaces will be considered: (1) The interface between the ODE integrator and the routines which approximate the spatial (boundary-value) derivatives in partial differential equations will be considered. In particular, the freedom to independently select the methods of integration for initial-value and boundary-value derivatives should be emphasized to users. (2) The interface between the general -purpose NMOL code and the routines which define the problem should be standardized as much as possible to facilitate the exchange of application problems. The essential elements (e.g., routines) to define a problem will be considered, and some suggestions given for a relatively standardized interface between the problem statement provided by the user and the NMOL code. 6. Comments on the Numerical Method of Lines (R. F. Sincovec) differential equations is a remarkable success. Moreover, this approach to solving partial differential equations has allowed the development of some very useful and modestly reliable general purpose software for PDEs. In large measure these successes are due to the availability of good software for solving ordinary differential equations. Unfortunately, one dimensional models are not adequate for simulating many physical phenomena and quite often one must use a model which is described or characterized by multidimensional partial differential equations. The fundamental concept of the numerical method of lines is essentially dimension independent, so in principle the technique may be directly applied to the higher dimensional problems. However, some of the practical computational aspects of the approach have prevented its widespread use of the higher dimensional problems which require the use of stiff ODE methods. For these problems, the limitations occur in the form of storage requirements and Jacobian matrix generation and solution time. The most frequently used stiff ODE techniques are those of Gear and the Jacobian matrix is required because a modified Newton method is used to solve the implicit nonlinear corrector equations. When a direct method is used to solve the resulting linear algebraic systems, the computation time per equation per time step increases drastically as the number of spatial dimensions increases. It would be extremely useful if reliable and efficient stiff ODE methods could be developed which do not require so much storage and computational cost per time step. In particular, it would be advantageous to use a nonlinear method, which requires only some portion of the associated Jacobian matrix. Alternatively, explicit methods with larger stability regions than those now in common use could be explored. It is conceivable that such ODE methods, which are basically different from the current backward differentiation formulas, could be developed - at least for some classes of problems. Another approach would be to investigate and develop other solution techniques for the nonlinear corrector equations. The goal would be to find adequate nonlinear equation solution methods whose cost in both storage and execution time is less affected by dimensionality considerations. Nonlinear SOR-Newton, Newton-SOR or ADI techniques could be potential candidates. The use of any of these methods would probably impose restrictions on the type of problems which can be solved and this needs to be investigated. Also, the relationships between the convergence criteria and the error and stepsize control in the ODE solver would have to be explored. Many problems encountered in practice can be reduced by the method of lines to a coupled differential -algebraic system. More attention needes to be paid to such systems. Finally, the ODE software itself should feature dynamic dimensioning. The use of the numerical method of lines for the solution of one dimensional partial W-M PANEL ON TESTING Moderator: T.E. Hull, Department of Computer Science, University of Toronto, Toronto,' Canada Panelists: M.B. Carver, Atomic Energy of Canada Limited, Chalk River, Ontario, Canada W.H. Enright, Department of Computer Science, University of Toronto, Toronto, Canada I. Gladwell, Department of Mathematics, University of Manchester, Manchester, England F.T. Krogh, Jet Propulsion Laboratory, Pasadena, California, U.S.A. L.F. Shampine, Sandia Laboratories, Albuquerque, New Mexico, U.S.A. Each panelist was asked to write a brief summary of his point of view and/or any comments he would like to make on this topic, and the following is what was received: M.B. Carver There is no doubt that a program of testing will produce the most useful information if it investigates the widest possible spectrum of equa- tion set type and size, accuracy level requested, range of magnitude of variables, etc., but, given a comprehensive test program and amassed results, the criteria used to assess or rate algorithms will depend heavily on the type of use for which an algorithm is to be selected. While a considerable market exists for extremely precise results, proba- bly the greater use of integration algorithms is in general engineering and scientific calculations in which three to five reliable digits suffice. Such algorithms usually reside in general sub- routine libraries or are embedded in user-oriented simulation packages, and the average user will not be equipped to select parameters to guide the integration algorithm. Thus the primary judgement criterion for such use is reliable performance in this accuracy range, and robust or self sufficient behaviour on encountering problems, such as dis- continuities, which are common in practice but traumatic in integration. W.H. Enright During the past several years a group of us at Toronto have been involved in the development, assessment and comparison of numerical methods for ODE's. These investigations have convinced us of the need for a collection of testing programs to aid in the evaluation of numerical methods, and the design and use of such a package is the subject of this paper. The collection of programs actually consists of two very similar packages; one for the assessment of non-stiff methods and another for the assessment of stiff methods. Before we discuss the test programs in more detail we must describe the type of applications we envisage for it and dis- tinguish two points of view that will be taken when interpreting the results. With the number of new methods being developed • publication, it Ls .'i formidable task for an at a referee to evaluate a now method. In order to publish an algorithm in a technical journal the author is expected to provide evidence of its performance characteristics compared to other methods. This evidence all too often takes the form of results of the new method relative to out-of-date methods on a few carefully selected problems. Refereeing such submissions is difficult and a referee is often reluctant to ask for more evidence since this often would require that an author obtain and implement a more appropriate reference method and run several new test problems. On the other hand, if a catalogue of standard bench- mark results were available for existing software, and if these results could be duplicated and results for new methods readily obtained, then we could require much stronger evidence of performance for new methods. Another use of the testing package would be that of a potential user with a particular type of problem. He may wish to determine which method is most suitable for him by benchmarking the performance of available methods on a set of his own model problems. There are two points of view that may be taken when assessing the performance of a numerical method. One may be primarily interested in the evaluation or assessment of the performance of an individual method, or alternatively he may be in- terested in comparing the relative performances of different methods. In the former case it is essen- tial that the performance criteria and performance monitors that are introduced have no effect on the method being assessed. That is, we must ensure that the method is not modified in any way, as such modi- fications might degrade the performance of the method. In order to satisfy this condition the testing package must be flexible and able to assess a particular method on how well it (the method) is accomplishing what it was designed to accomplish. On the other hand, when one is comparing methods, the various measures of reliability and efficiency can only be compared if the methods being compared are accomplishing the same task. One way to ensure that this is the case is to modify the methods so that they are all attempting to accom- plish the same task. This was the approach adopted in Hull et al (1972), Enright et al (1975) and Enright et al (1976). This approach allows one to compare the relative performances of different techniques or approaches but it must be acknowledged that modified methods are compared, and these modi- fications might have a greater effect on some methods than others. 45-1 Lng Initial value solvers ran a varietv These include r visit co the University of Toronto, stions by inexpi ■; users which I fv Le to inswer without extensive testing, □ my involvement as contributor oi the ordi- L equations chapter of the NAG iv. I have tested initial value solvers in two v IJNCl.ASSllll.l) 20. Sec urity < lass (This I \< LASS1F1ED 21. No. oi Pages 2?. ■ USCOMM-DC 40329-P7 m