OPTIMALLY STABLE FIFTH-ORDER INTEGRATION METHODS : A NUMERICAL APPROACH

In this paper, the results of the first detailed and systematic study of the family of fifth order implicit linear multistep methods requiring function evaluations at four backpoints are given. Also described is a practical and efficient nonlinear optimization procedure which made it possible to locate in a precise manner the methods of this type which possess optimal ranges of relative stability in the sense-that the corresponding relative stability regions contain maximal disks centered at the origin.


I. INTRODUCTION.
The advantages of using linear multistep methods in the numerical integration of ordinary differential equations have long been recognized. While the analytic theory of these methods has been discussed adequately by several authors (see [2] or [5], for example), the results concerning the practical problem of locating optimal methods have been less than satisfactory. In particular, while it has been recognized that a determining factor in the choice of such a method is often that of relative stability, most authors have limited their attention to a study of absolute stability. As a result very little information concerning relative stability is now available, especially for methods of order greater than four.
In [8] an algorithm was given which may be used to determine numerically the largest open line interval I on which a given method is both convergent and relatively stable. The corresponding interval was computed for several well-known methods. In [i0], from several families of Adams-style predictors and correctors, the algorithm was used to locate methods for which the size of intervals of the form (-,) contained in I is a maximum. The method was also modified and used to locate methods with increased complex stability regions. In particular, the complex stability characteristics of fourth order correctors requiring three backpoints and fourth order predictors requiring four backpoints were studied. Methods whose relative stability regions conrain maximal disks centered at the origin were located in each case. In this paper the stability characteristics of fifth order methods will be considered from this same point of view. For the convenience of the reader, several pertinent definitions and some basic terminology will first be discussed.
Many details will be omitted since they are quite lengthy and are also relatively straightforward. Any omissions are discussed at length in [i0]. It is emphasized that use of the word "optimal" in the present study applies only within the context of the methods considered. Indeed, the conclusions of this paper are not necessarily applicable to more general families of methods (see [9], however), particularly if the desired optimization is in respect to a different method characteristic (e.g., absolute stability). It is noted, however, that methods obtained as in the present study also tend to preserve other desirable method characteristics (see [9], for example).

DEFINITIONS AND PRELIMINARIES.
Consider the initial value problem where y' dy/dt and y, fer n The existence of a unique solution to this initial value problem will be assumed. An analysis of the stability of a numerical method applied to this problem may be carried out in the usual manner by studying the stability of the method applied to the collection of initial value problems (2.2) where i is an eigenvalue of (f/y) the Jacobian matrix of f with respect to y.
A numerical approximation to the solution of (2.2) may be obtained as are constructed, which, under suitable conditions follows. Sequences Yn+lj 1 _(0) is obtained using will yield the desired solution to (2.2). For each n, Yn+l a predictor method. Given Yn+l" x 0 + kh, Yk Y(Xk)' Yk f(xk,Yk) a k and b k are suitably chosen constants with b_l# 0, and h is the fixed step size. The points (x,yk) are referred to as backpoints. If equation (2.3) is satisfied exactly for poly-r+l nomials of degree less than r+l but not for the polynomial y x the method will be said to have (polynomial) order r. In this case, the method will also be said to be of type (r,q+l). When applied to the initial value (J) ( It is expected that the results of the present study will also serve as a guide in future studies of fifth order predictor-corrector combinations. At any rate, it will be assumed for the remainder of this paper that (2.3) is a genuine corrector method. (See also section 5 of [6] .) The definition of stability which will be used, that of relative stability, is taken from [5]. The  and one of the roots, which is denoted by r0( may be written as r0() 1 + e + 0(e 2) exp(e) + 0(52) for II sufficiently small, r0(e will be referred to as the principal root of (2.4). Denote the remaining roots by r I(),. .,rq(e). (2.2) is said to be relatively stable for lh in case i l,...,q subject to the restriction that when equality holds, r i() must be simple.
3) is relatively stable for e 0, it will be said to be initially stable. An analytic description of the methods under consideration and the necessary local error analyses will now be given.

ANALYTICAL DESCRIPTION OF FAMILY OF (514) MULTISTEP INTEGRATION METHODS.
Throughout this section, the concept of the relative stability disk for (2.3) will be employed. By this will be meant the largest disk centered about the origin in the complex plane and contained within the region for which ( (See [5]. For each value of (al,a2,a3) this system of equations has the unique solution where (a 0 b_l,b0,bl,b2,b3)  It is pointed out that in [4], Newbery has given a collection of (q+2)xq  (3.7). In this case, the principal coefficient of the local truncation error is given by (3.5). Furthermore, if y y(x) belongs to C6[a,b], the local truncation error is given by E(al,a2,a3)h6y (6) () (3.8) for some in (Xn_3, Xn+l) if -19 + lla I 8a 2 + 243a 3 < 0.
It is remarked that the last condition is satisfied for all points (al,a2,a3) in A with a 3 < 0.
A discussion of the last assertion of the theorem depends on the notion of the influence function w(s) for (3.1) ( [5]). w(s) is given by   Therefore, it follows that none of the planes intersect A at points with a 3 <_ 0.
[If 8 0, the plane intersects the last two lines at the point (1,-1,0)  In both cases it is possible to verify that none of these planes intersect A at points for which -19 + lla I 8a 2 + 243a 3 < 0.
In doing so, it is helpful to observe the "motion" of each family of planes as increases from 0 to I. This completes the proof of the theorem.
The purpose of the present study is to locate the method described in the theorem above for which the corresponding relative stability disk is as large as possible. It may be argued that more appropriate criteria exist by which to optimize. For example, the roots of (2.4) rather than the coefficients may be chosen as the free parameters, as in [3]. However, it is noted that the manner in which linear multistep methods are implemented in several widely available automatic computer software packages permits an easy implementation of other, methods via the redefinition of the method coefficients. The present approach may therefore offer advantages for purposes of testing and comparison.
Related to the question above, many readers will also find the following remarks to be of interest. The roots of (2.4) must satisfy several natural constraints. In particular, they must satisfy the well-known root-condition ( [2]) when h=0 (that is, they must be initially stable). This condition imposes constraints on allowable values in the ala2a3-parameter space of present interest. For example, (2.4) has a double root of r=l for h=0 if and only if 1 + a I + 2a 2 + 3a 3 0. In the present study, we found it necessary to avoid this plane. Indeed, when other families of methods were considered in related studies, it was found that the optimization paths tended to follow this plane in certain regions (see [9]).

NUMERICAL RESULTS.
In this section several numerical results obtained for the family of methods under consideration are presented. The manner in which these results were obtained will first be discussed. In order to study methods in the plane a 3 e, linear combinations aP* + (l-a)P (a(l-a3) a2,a3) where P* (l-a3, a2, a3) and P (0, a2, a3)were considered. As a first approximation to the manner in which the algorithm will ultimately be used to study families with several parameters, the algorithm was also modified and used in a modified steepest descent optimization procedure. Although the algorithm thus obtained tended to exhibit the usual cycling behavior which is characteristic of such a procedure, it was possible to locate the best method given in Table i  (This is particularly true near critical points.) It will of course be interesting to determine the extent to which it is possible to implement the present algorithm in a more sophisticated optimization procedure when higher order families are considered. It is remarked that, at least at this point, use of the algorithm of [8] (as modified in [i0] and herein) has permitted what seems to be the first attempt to seriously treat the determination of the relative stability characteristics of (3.1) strictly as a nonlinear optimization problem.
Several of the other programmatic aspects of the algorithm used to obtain the present results will now be discussed. When a good approximation to the principal root at the point e x 0 + i Y0 was not available, it was found that using exp(x + i y0 as an approximation worked well. (It is remarked that, in their definition of relative stability, many authors replace r0(e) with exp(e).) At any rate, after it was decided that a method in the plane a 3 was "optimal", the usual root-locus method was then also used to compute the corresponding stability region directly. This additional computation also serves as a safeguard against difficulties related to the numerical illconditioning of polynomials as functions of their coefficients.  [i0] for the values a -0.640 and a -0.675 as mentioned at the beginning of section 3.
The method corresponding to the entry in Table i   From the family of (5,4) iterative correctors, the method (5.2) which has optimal relative stability properties in the sense that its complex region of relative stability contains a maximal disk centered at the origin has been located. Besides being a feasible alternative to the standard fifth order methods currently being used, the method is also a natural candidate as a starting point for a detailed investigation of the relative stability characteristics of fifth order predictor-corrector combinations. Furthermore, the algorithm used to obtain this method may be generalized and applied to other classes of linear multistep methods. The results given in this paper indicate that it is possible to answer important questions, previously considered intractable,