Software based on explicit RK formulas

MATHEMATICS ELSEVIER

Applied Numerical Mathematics 22 (1996) 293-308

Software based on explicit RK formulas L.F. Shampine,

I. Gladwell *

Department of Mathematics, Southern Methodist University, Dallas, TX 75275, USA

Abstract

A personal view of the historical development of software based on explicit Runge-Kutta formulas is presented, and some future developments are suggested. Keywords: History; Software; Explicit RK formulas; Step size control; Global error; Stiffness; Continuous

extension; Evaluation

1. I n t r o d u c t i o n

This paper presents a personal view of the historical developments we have observed in three decades of work on software for ordinary differential equations (ODEs) based on explicit RungeKutta (RK) methods. Continuing this historical treatment, we describe some current work and trends for the future. Our main concern is with general-purpose software, though we make some observations about solving ODEs arising in special circumstances. Our work with the SLATEC and NAG libraries and specific codes written by ourselves and others serve as mileposts. This survey is necessarily brief; basic definitions and fundamental algorithms can be found in other papers of this volume and some issues are treated in more detail but from a similar viewpoint in [43]. The history of developments that led to features found in the suite of codes RKSUITE [7,8] provides the skeleton of this survey. It is fleshed out with specific developments in formulas, initial step size selection, the selection of step sizes for efficiency and control of the local error, global error assessment, diagnosis of stiffness, continuous extensions, and evaluation of software. We conclude with some "things to do", useful capabilities that could be added at this time, and with observations about the influence of new computing environments on the development of ODE solvers based on explicit RK methods. * Corresponding author. E-mail: [email protected].

0168-9274/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved PH S0168-9274(96)00039-6

294

L.F. Shampine, L Gladwell / Applied Numerical Mathematics 22 (1996) 293-308

2. Formulas In deriving RK formulas, at first the game was to get the highest order possible with a given number of stages. Since much of this work was done by hand, the preference for coefficients that are the ratio of small integers or a root of such a ratio is understandable. From merely finding a formula, attention turned to finding an accurate, efficient formula (one with "small" truncation error coefficients) that had a large absolute stability region. When the advantages of estimating and controlling the error were recognized, the game changed in a major way. Because the first formulas were derived without consideration of error estimation, the error was estimated by means of a general principle exploiting the order of consistency of the integration formula, usually "doubling" (see below). As error estimation is viewed now, each step is taken with a pair of formulas and the error in the lower order formula is estimated by comparing its result to that of the higher. The cost of a step is then not just the number of stages in the formula used to compute the solution, but also those needed for the error estimate. Though some workers, notably Sarafyan, developed competitive formulas contemporaneously and others, notably Merson [33], had developed them much earlier, a pair of formulas of orders 4 and 5 due to Fehlberg [15] became generally accepted as the "best" at this popular order. This was largely due to the paper of Hull et al. [28] comparing methods for solving ODEs and to the subsequent appearance of a production-grade FORTRAN 77 code, RKF45, of Watts and Shampine [49]. RKF45 was written for a software library that evolved into the SLATEC library, where it is still available. This code has appeared in a number of textbooks, papers, and packages (sometimes without attribution). It has been translated into other languages, including C, C++, and MATLAB and been embedded in application programs. At roughly the same time, important codes were developed for the two most widely distributed libraries of mathematical software, namely D02PAF, based on the Merson (3,4) formulas [33], for the NAG FORTRAN 77 library [17] and DVERK, based on a (5,6) pair due to Verner, cf. [54], for the IMSL Math library [29]. The latter code is extant, and available from other sources. Both the NAG and IMSL libraries now include RKSUITE. The collection of codes comprising RKSUITE is based on two pairs of formulas by Bogacki and Shampine and one by Prince and Dormand [1,2,36]. (A unique, undocumented aspect of the implementation of the (4,5) pair in RKSUITE is that, as a double check, two independent estimates of the local error are formed at each step.) The three pairs of formulas in RKSUITE were developed in the light of experience with the previous generation of formulas and an improved understanding of the issues. The search for still "better" formulas continues, see for example [52]. An issue arises when each step is taken with two formulas. It is the error of the lower order formula that is estimated, but why discard a result that is believed to be more accurate? Advancing the integration with the higher order formula was called local extrapolation by Shampine [37]. Local extrapolation raises the order of the integration at no extra cost, but the stability of the integration is that of the higher order formula and the error estimated is not that of the formula actually used. After a long polemic, practice settled on the use of local extrapolation. This decision influences the derivation of a pair because of its significance for stability, but also for another reason we take up now. The first stage of a RK formula is the evaluation of the slope at the beginning of the step. An idea called FSAL (First Same As Last) by Dormand and Prince [1 I] is to use this stage from the next step in the companion formula for error estimation in the current step. If the current step is a success, and most are if "safety factors" in the code are chosen appropriately, the stage will be reused in the

L.F. Shampine, I. Gladwell / Applied Numerical Mathematics 22 (1996) 293-308

295

next step and may fairly be regarded as "free". The FSAL idea is independent of the idea of local extrapolation, but it requires a decision about whether local extrapolation is to be employed with the pair. The additional flexibility provided by a free stage from FSAL is valuable, but at first attention was directed to minimal stage pairs. More recently, Bogacki and Shampine [2] and Sharp and Smart [51] observe that the additional flexibility of extra stages can more than pay for themselves on an equal cost basis. There are, however, practical objections to the large step sizes needed for the efficient use of a "large" number of stages, so this approach cannot be taken too far. Nowadays the derivation of formulas relies heavily upon computer assistance, both symbolic and numeric. With computer tools and a better understanding of the issues that govern reliability, efficiency, and stability, we are able to search for "optimal" pairs. The search is further complicated by the need for a sufficiently accurate, yet efficient, continuous extension, a capability we take up below. So far we have been discussing formulas to be used in general-purpose software. There is always a question as to whether a particular class of problems is sufficiently important to warrant developing software that exploits its particularities. Problems arising naturally as second order equations without first derivatives present in the definition of the ODE are of sufficient importance that they received special attention right from the beginning. It is possible to derive formulas that solve problems in this important class more effectively than by the standard device of reducing them to first order systems. A modern package of codes for this problem using explicit Runge-Kutta-Nystr6m formulas and with design similar to RKSUITE was developed by Brankin et al. [5]. A related area of contemporary interest is the derivation of formulas that are symplectic for long-time integration of Hamiltonian (symplectic) problems. Much remains to be understood about long-time integration of ODEs, but in some applications it seems important to integrate with formulas that preserve certain qualitative features of solutions of the ODE. Though no general-purpose software for this problem is at present available, there are now formulas with error estimation [22] that could replace a conventional RK formula pair. In some contexts it is still common to use fixed step size. For instance, in the method of lines (MOL) for approximating the solution of time-dependent systems of partial differential equations (PDEs), the spatial variables are discretized, with or without attempting to control their error, then the resulting system of ODEs is integrated in time. Often there is little interest in controlling the temporal accuracy, or it is just assumed that the temporal error will be sufficiently small it can safely be ignored, so a fixed time step is used, usually chosen to ensure stability for all time by satisfying a global CFL condition. Because the systems arising in the MOL are typically large, there is great interest in methods with low memory requirements. Although low memory methods with special structure have been proposed that allow a moderately high order, the RK formulas in wide use in the MOL literature are of low order because they involve only a few stages and they are efficient for the modest accuracies that are typical. Also, these formulas provide the large scaled stability intervals that are often needed. More specific to the context is the search for RK formulas with special properties that hold even at the longest step sizes permitted by the global CFL condition. For example, there has been a search for formulas that are variation-diminishing [32,53], a property that helps assure an approximation with appropriate qualitative behavior to certain hyperbolic PDEs in conservation form. A related question [10] is an appropriate way to treat boundary conditions when evaluating the RK stages. Other qualitative properties that may be important when integrating time discretizations of PDEs are the dissipation and/or dispersion of the formula [27].

296

L.F. Shampine, I. Gladwell / Applied Numerical Mathematics 22 (1996) 293-308

3. Choosing the step size At first codes integrated with a fixed step size that the user would guess. Later codes adapted the step size to the problem, but they still asked the user to provide the initial step size. This is at best a nuisance. Generally the user does not have the information needed to make a good guess. If the initial step size is too big, the results may not be reliable, and if it is too small, the computation may be inefficient. Today production-grade codes select automatically the initial step size, though they do allow users to influence the choice. Because most codes adjust the step size quickly, the practical issue is recognizing quickly the scale of the problem. The first non-trivial scheme for selecting automatically an initial step size appears to have been that of [45]. The solver has more information at its disposal than the user does, but still very little. Shampine's rule of thumb uses all the data available from the problem definition and the order of consistency of the method to produce a step size with the properties expected for a typical problem. Variations on the theme are used now in many codes. However, more can be done. The somewhat elaborate scheme of RKSUITE obtains a tentative initial step size with the rule of thumb and then cautiously takes a first step. Whilst computing the stages of the first step, a Lipschitz condition is monitored and the step size reduced as necessary to make the product of the step size and a Lipschitz constant of modest size. The two devices inspire some confidence in the error estimate computed for this first step. In the usual fashion this estimate is used to adjust the step size. However, rather than go on to the next step, the step is repeated until an on-scale step size is found. Only then does the integration proper begin. This starting step size algorithm is essentially that proposed in [19]. It combines the best of corresponding procedures in RKF45 and the NAG RK-Merson code, D02PAE Providing a good guess for the initial step size when using a variable-order Adams or BDF code is especially difficult for users because the codes start with a surprisingly low order, typically one. Several authors, e.g., [6,16], have suggested taking a first step with an explicit RK pair. A continuous extension associated with the pair could be used first to select an appropriate step size and order and then to generate the starting values needed to begin the integration with the multistep method at this order. Although this idea is attractive, its realization is somewhat complicated and is not yet seen in production codes. Early codes estimated the local error by doubling. If the step size is halved when the step is a failure, it is possible to reuse a number of stages. Perhaps by analogy with concerns special to multistep methods and perhaps by symmetry, doubling the step size when possible was common. One of the attractions of RK formulas is that it is easy to adjust the step size at every step. Restricting the choice to halving and doubling means that the integration does not use the most efficient step size. It also means that the global error does not behave as smoothly as we would like. Furthermore, a number of workers found in substantial experiments that the savings obtained by halving on a step failure were illusory. This is because halving the step size is sometimes too much, which is inefficient, and is sometimes not enough, leading to subsequent step failures, which is very inefficient. The authors of production-grade codes quickly shifted from the tactic of halving to that of an "optimal" reduction, optimal in the sense that the new step size is expected to result in a norm of the local error just less than the accuracy tolerance. As the general principle for error estimation was abandoned in favor of schemes tailored to the basic formula, the issue became moot. There was a lengthy polemic about whether local error ought to be controlled by error per unit step (EPUS) or error per step (EPS). The principal attraction of EPUS was that there is a simple and

L.E Shampine, 1. Gladwell / Applied Numerical Mathematics 22 (1996) 293-308

297

general proof of convergence and a bound on the error that is proportional to the tolerance. In this generality the results are not true for EPS. Still, the production codes of the time used EPS and they were clearly more effective than codes using EPUS. It was some time before it was appreciated that the crux of the matter is that the step size chosen in the production codes is efficient, which is to say that it is about as big as possible for the tolerance. With this additional assumption, the desired behavior can be shown to hold for EPS, too. A number of the important production codes based on EPS also used local extrapolation, and it was eventually appreciated that EPS with local extrapolation corresponds to a generalized EPUS for which a modification of the analysis for EPUS establishes the desired tolerance proportionality results. In the event, the production-grade explicit RK codes in wide use today control the local error by EPS and use local extrapolation. The results discussed in the previous paragraph are bounds. The earliest investigations establishing an asymptotic expansion for the error of integration with an explicit RK formula are for constant step size. It is easy to extend the analysis to variable step size when the step size at a given point in the integration is specified by a step size selection function. However, for a long time no one asked if the step sizes selected automatically in the codes could be described adequately by such a function. This turns out to be the case [39], and the expansions can be used to understand the implications of controlling the local error by EPS and EPUS and of doing local extrapolation. With reasonable assumptions it is possible to analyze the behavior of the global error to leading order in codes that select the step size automatically. For instance, when the local error is controlled by EPS and local extrapolation is used, the global error is proportional to the tolerance. Another useful conclusion is that local extrapolation improves the efficiency of the computation. The principal assumption is that the step size selected is efficient. A less critical, but still important, assumption is that the initial step size is onscale. In part this explains our interest in these issues when developing step size selection algorithms. The earliest step size selection algorithms assumed that the behavior of the local error on repeating a step with a smaller step size would be qualitatively the same. A robust algorithm must always be concerned with the possibility that basic assumptions are not true and a number of tactics have been used for identifying and dealing with this. The same assumption about the behavior of the local error was made when selecting the step size for a step following a successful one, though this is rather less plausible and greater caution is appropriate. Zonneveld [59] was probably the first to exploit information garnered at previous steps in the step size selection algorithm. The idea is to approximate the behavior of the leading term in an expansion of the local error using estimates of the local error at one or more previous steps. Watts [56] proceeded similarly and developed an algorithm that is rather more efficient than the traditional one, an algorithm that is the basis for the one used in RKSUITE. An interesting approach to the task taken by Gustafsson et al. [20,21] is to regard step size selection as a (feedback) control problem. This is mainly a difference in viewpoint, for the actual algorithms of Watts and of Gustafsson are remarkably similar. Finally, the traditional step size adjustment algorithm is notably inefficient when an explicit RK method is applied to a stiff problem. Interesting work presented in [24] examines the roles of the method and the error estimator in this situation.

4. G l o b a l e r r o r a s s e s s m e n t

In modem codes step sizes are chosen to control the local error, but it is the global error that concerns the user. Generally, it is not possible to control the global error directly in step-by-step

298

L.E Shampine, L Gladwell / Applied Numerical Mathematics 22 (1996) 293-308

methods. Usually, there is no way of knowing how much error can be allowed at the current step and still achieve the desired global error at future points of interest, because this depends on the stability of the differential equation in the neighborhood of the (unknown) solution. The stability of routine problems is such that, in carefully crafted codes, the global error is often comparable in size to the local error tolerance. When it is not, it usually behaves roughly proportionally to changes in the local error tolerance. Most users are content with this, but it would be useful to verify that the global error is behaving as expected, so some attention has been given to assessing it. A major reason for an early general-purpose code that estimates global error, GERK [48], was to use it for the solution of Sturm-Liouville eigenvalue problems by shooting. Competing methods estimate the errors in the computed eigenvalues. To do the same with shooting, the code for initial value problems must be able to estimate the global error of its solution. GERK estimates global error by global extrapolation. A problem is integrated with the Fehlberg (4,5) pair much as in RKF45. A parallel integration is carried out simultaneously using step sizes that are half those selected in the primary integration. In principle the parallel integration might be done with double the step size of the primary integration, but this is not possible in practice because often the integration would not be stable. The error of either integration can be estimated. The object of this code was to solve the problem and provide an estimate of the error in its solution, so it reports the more accurate result. Accordingly, the local error tolerances supplied to the code are adjusted internally so that the accuracy of the parallel integration is comparable to the accuracy desired by the user. This is different from wanting to check the accuracy of the results of an integration by setting a switch, as is done in a similar code in the RK-Merson package in the NAG library [17]. For such a code the integration at half the step size is purely for estimating the global error, making global error estimation comparatively expensive. Shampine and Watts [48] pay particular attention to the reliability of global error estimation. It is difficult to estimate the global error reliably at crude tolerances when the asymptotic approximations underlying basic algorithms can be poor because the step size is "large", and at stringent tolerances when they can be poor because of finite precision arithmetic. There are important theoretical limitations on the whole approach of global extrapolation. For instance, the occurrence of a single instant in the interval of integration when the solution is not sufficiently smooth for the postulated asymptotic expansion ruins the estimate at all subsequent times. There have been a number of investigations of less costly ways of estimating the global error, but reliability remains a critical issue. In the RK context, see especially the work of Dormand and Prince [ 12,13] who, for smooth problems, endeavor to compute an "accurate" global error estimate throughout the integration, not just an order of magnitude estimate. There is a serious conceptual difficulty in deciding just what to ask of global error estimation. In our experience, what is usually desired is an assessment of the general size of the error throughout the interval of interest rather than estimates at specific points. The global error at a given point may be atypical and it is obviously difficult to estimate even crudely near a point where it vanishes. The eigenfunctions of a Sturm-Liouville problem are illuminating. A small phase shift in the computed eigenfunction is generally not of practical importance as long as it is not so large as to affect integrals involving the eigenfunction or to affect the computed eigenfunction's appearance. Yet a small phase shift can imply a "large" error near points where either the eigenfunction or its numerical approximation vanishes.

L.E Shampine, 1. Gladwell / Applied Numerical Mathematics 22 (1996) 293-308

299

In RKSUITE we prefer to a s s e s s the global error rather than estimate it. Although the code estimates the global error at each step internally, what it reports is the RMS error made in the integration from the start to the current time. The code reports the RMS error rather than estimates at specific points because the average over the interval of interest provides what we think most users want and it should be estimated far more reliably than errors at specific points. The capability is provided as a check, so the numerical solution of the ODE is exactly the same whether or not the assessment is switched on. Because it is a check, a way of assessing global error was to be grafted on to a method for solving the problem that was chosen for high performance. The possibilities were limited by the fact that RKSUITE implements three pairs of formulas, one of which is not accompanied by a continuous extension. Global extrapolation is an obvious possibility, but this code does not use global extrapolation. What it does might be described as a reliable version of the method of reintegration. Specifically, a parallel integration is done with a step size that is a fraction of that used in the primary integration and the accuracy of the primary integration is estimated by comparison. This scheme is applicable to any RK method and produces a good estimate provided only that integration with a much smaller step size gives a somewhat more accurate solution. In RKSUITE half the step size is used in the parallel integration in the cases of the (4,5) and (7,8) pairs and one third in the case of the (2,3) pair. The low order formula is treated differently because it is used most often at crude tolerances which can imply step sizes for which asymptotic approximations may be poor. Also, in this range of tolerances the assumption that a parallel integration with half the step size is a better approximation to the true solution may not be true. This is a relatively expensive way to estimate the global error, but it is reliable because it requires so little of the behavior of the secondary integration. Local errors are estimated in both the primary and secondary integrations. They are compared so as 1o determine whether the two solutions being tracked have diverged to the point that the global error estimate might be unreliable. An assessment of the global error is useless as a check if you cannot count on it, so if there is any doubt, the code reports that it no longer trusts the assessment and terminates the run. The solvers in RKSUITE also provide the maximum global error seen in any component of the solution and the point in the integration where it occurred. This can help identify difficulties with the problem and its solution. An attractive general approach to global error assessment is to monitor the residual (defect). For methods that have an inexpensive, accurate continuous extension, this is appealing because the size of the residual might be determined reliably even at crude tolerances. Current research has focussed on estimating the residual inexpensively [14], but the fruits of this research need to be combined with reliable tactics for computing a credible assessment when asymptotic results are of limited value. An important step towards credibility was taken by Higham [25]. It is to be appreciated that the residual is intermediate between the local error that the codes control and the global error that users want to control. Numerical analysts with an understanding of backward error analysis might be reassured by a report that the residual is "small", but it is not clear that a "typical" user concerned only about the solution would share this feeling.

5. Diagnosing stiffness Explicit RK methods are not effective for the solution of stiff problems. Because a user often cannot tell in advance whether a problem is stiff, it is useful to be able to diagnose stiffness. Stiffness is a

300

L.E Shampine, I. Gladwell / Applied Numerical Mathematics 22 (1996) 293-308

murky concept. In practice what is meant is that a good code intended for stiff problems, typically based on the backward differentiation formulas (BDFs), is much more efficient than a good explicit non-stiff code, typically using RK or Adams formulas. There are many factors that play a role in this and, correspondingly, it is difficult to decide automatically whether a problem is stiff. If a code based on an explicit RK formula is applied to a stiff problem using a constant step size that is too large, the integration "blows up". As first pointed out in [38], automatic selection of the step size so as to control the local error stabilizes the integration. The argument shows that there are many step failures and the average step size corresponds roughly to being on the boundary of the stability region of the method. These have been accepted as the facts though the theory remains incomplete. These facts can be used to obtain an inexpensive indicator of stiffness, namely an expensive integration with a high percentage of step failures and a step size that is roughly constant. More direct attempts to recognize stiffness rely mainly on comparing the behavior of two formulas with different stability regions in a way rather like local error is estimated. Combinations of these and other tactics are used in RKF45 and in the NAG RK-Merson stiffness detection code, see [46] and [17], respectively, for details. A primary factor in stiffness is the (absolute) stability of the formula. The eigenvalues of local Jacobians are critical to recognizing whether stability is holding back the step size. Although it is out of the question for general-purpose software to consider all the eigenvalues, it is practical to look at a couple of the largest in magnitude. This is not as simple as the textbooks suggest because the stability regions of popular methods are not (half) disks, meaning that the precise location of the eigenvalues in the complex plane matters. Problems with eigenvalues near the imaginary axis are especially hard to classify because for many formulas the boundary of the stability region changes sharply there. It is not even clear what stiffness means when the dominant eigenvalues are in this area because perturbations to the solution are not strongly damped. Problems with eigenvalues with "small" positive real part are also hard to classify. Despite these qualifications, in RKSUITE we found that a reliable diagnostic of stiffness [42] can be made at reasonable cost by making the relatively expensive computation of the dominant eigenvalues of a local Jacobian (using a nonlinear power method) periodically and then only when very inexpensive tests suggest that stiffness is possible.

6. Continuous extension

Adams formulas and the BDFs are based on interpolating polynomials [43], and this makes it clear how to approximate the solution between steps. Modem codes based on these formulas integrate with about the largest step size that will provide the specified local accuracy and keep the computation stable. They obtain output at specific points by evaluating the interpolant that underlies the discrete variable method. Accordingly the number and placement of output points has little effect on the cost of the integration. An important application of interpolation is to event location. An event is said to occur at time t* if the solution y(t) of the ODE satisfies the algebraic equation g(t*, y(t*)) -- 0 for a given event function g(t, y). There may be many event functions associated with an initial value problem. Event location is difficult, both in theory and practice. Among the difficulties are that the task may not be well-posed and there is no way to be sure of even realizing the existence of a zero of a general algebraic equation, much less of computing the first zero. When it is realized that an event occurs in the course of a step, it is often necessary to locate the event very accurately. This can

L.E Shampine, I. Gladwell / Applied Numerical Mathematics 22 (1996) 293-308

301

be found efficiently using values for y(t) computed inexpensively by evaluating an interpolant when solving the algebraic equation for t* and y(t*). For explicit RK methods there is no natural way of approximating solutions between steps. However, because they provide both the solution and the slope at each step, it is natural to use an Hermite interpolating polynomial as a continuous extension. This was found to work well enough at low orders, and it provides the continuous extension of the (2,3) pair [1] of RKSUITE. Unfortunately, there are difficulties at even moderate orders because it is necessary to interpolate the results of several steps to preserve the accuracy of the method. Gladwell [17] discusses the difficulties that arise when interpolating over just two steps. One is that to get an accurate interpolant, it is necessary to control the rate of increase of step size, vitiating an important advantage of RK methods. Also, there is an obvious difficulty at the start of the integration. A more fundamental difficulty arose with the exploitation of higher order formulas. RK formulas of even moderate order make many more function evaluations per step than Adams formulas and BDFs. To be competitive they must take a step that is much longer. The difficulty is that interpolation over several steps can be justified in a limit of the step sizes tending to zero, but at practical tolerances it is often found that the step sizes appropriate to the RK formula are too large to achieve the same accuracy with the interpolant. A paper of Horn [26] was the first to present a modem approach to continuous extensions. This matter can be viewed in two distinct ways. In the view taken by Horn, a family of formulas is created, one for each point in the span of a step. Each member of the family is a linear combination of the stages used in the basic formula plus other stages as necessary. By virtue of reusing stages, it is possible to approximate the solution anywhere in the span of the step with a small number of extra stages. Another view is based directly on interpolation [40]. The intermediate values of a RK formula provide approximate solutions in the span of a step that are generally, but not always, of low order. Extra stages can be used to obtain intermediate solutions that are of the desired order, along with approximations to the slope there, and these values can be used along with values at both ends of the step to form an interpolant within the span of a single step. It is of practical importance that a continuous extension is required only after the step is judged successful. Further, the extra stages for forming a continuous extension are required only when intermediate output is needed in the span of the step. Once these extra stages have been formed, the number of intermediate output points in the step has little effect on the cost. These considerations are crucial to event location where, normally, events occur on only a very small proportion of the steps but many evaluations of the interpolant in the step may be needed to locate the event accurately. The accuracy of the intermediate results can be analyzed using either interpolation or truncation error theory and related to the accuracy of one of the pair of formulas. There is a conceptual difficulty-what is an appropriate order for the extension? To be concrete, let us think of a (4,5) pair with local extrapolation. The actual error of the numerical solution is composed of the global error of the integration to the current point, which is of order 5, and the local error of the current step. If the extension is of order 4, meaning that its local error is of order 5, the intermediate result will be of the same order of accuracy as the result at the end of the current step. However, it might be argued that the error being controlled is of order 4 and that is the order that ought to be matched rather than the order achieved by local extrapolation. In this view, an extension of (global) order 3 would suffice. Another view is that the integration is advanced with the formula of order 5, so the order of the extension should also be 5 if we want extensions that are truly comparable in accuracy to the formula used to compute the numerical solution at mesh points.

302

L.F. Shampine, L Gladwell / Applied Numerical Mathematics 22 (1996) 293-308

The BS(4,5) pair [2] of RKSUITE involves more than the minimal number of stages for its order and results in a very accurate step. Because of this a relatively large number of extra stages is needed to get interpolation accuracy consistent with the fifth order formula. Bogacki and Shampine considered not only the order of the error, but also its size relative to the error at the end of the step and its distribution within the span of the step. A remarkable property of the interpolant they develop is that, to leading order, the distribution of the local error is independent of the ODE. Indeed, this property also holds for the Hermite interpolant associated with their (2,3) pair [1]. This property implies that the error in the interpolant everywhere in each interval is controlled by the local error control of the solution value at the end of the interval. The Prince-Dormand (7,8) pair [36] of RKSUITE is not accompanied by a continuous extension. Interpolants proposed for pairs of this order are rather expensive, and when RKSUITE was written we did not think that a continuous extension was as important in the range of tolerances that are typical at this order. For these reasons RKSUITE does not at present have a continuous extension for this pair though extensions of appropriate order are now available, e.g., [55]. Graphical display of solutions has become the standard in many contexts, e.g., in a computing environment like MATLAB and on graphing calculators such as the Texas Instruments TI-85. In such contexts a modest accuracy is appropriate. With a relaxed tolerance, modem formulas take steps too long for the "natural output", the values computed at each step, to result in a smooth graph formed by joining output points with straight lines, as is typical of plotting packages. Polking [35] makes this observation about the MATLAB (4,5) pair. The new generation of MATLAB ODE solvers [47] deals with this difficulty by evaluating a continuous extension at a number of points in the span of a step. This raises questions of how many points are needed for a smooth graph and where they should be placed. Another issue that becomes relevant in such a computing environment is the shape of the solution. Because derivatives of interpolants converge to derivatives of the solution, qualitative properties of the solution like monotonicity and convexity are true of the interpolant as the step size tends to zero. However, this need not be the case at relaxed tolerances and it is then necessary to impose these properties on the interpolant, see, e.g., Brankin and Gladwell's work [4] where they propose a combination of Hermite polynomial and rational interpolants, fitted to the solution and derivative values typically available from an ODE solver, to preserve shape, and they calculate how frequently the interpolants must be evaluated in a graph plotting environment.

7. Evaluating methods Typically a user has a choice of solvers and first must decide which to use. When writing a solver similar choices must be made about the choice of formulas and the details of the algorithm. Testing and analysis help with these choices. Early on, testing involved solving a few problems with several codes and comparing the accuracy of the numerical solution and/or the cost in terms of CPU time and/or function evaluations. As the limitations of such comparisons became better understood, reports of more systematic testing began to appear in the literature. A landmark was the paper of Hull et al. [28]. A set of problems was assembled for their tests that still serves as a benchmark. A severe difficulty in evaluating software is that two codes can be rather similar, but apply to different classes of problems. Besides, when both are applicable, they may not attempt to solve the problem in the same sense or attack them in the same way. Hull et al. aimed to compare methods, in contrast to implementations. Because of this they were willing to alter the codes so as to make them solve problems in the same

L.F. Shampine, I. Gladwell / Applied Numerical Mathematics 22 (1996) 293--308

303

sense, alterations that in some cases had serious implications for the effectiveness of the codes. An interesting theoretical approach to comparing RK formulas is the use of effectiveness theorems based on theoretical cost criteria. This was developed in a number of papers published by the Toronto group including one specifically on RK formulas [31]. See [30] in this volume for more details. Gladwell et al. [18] took the opposite approach and compared codes as "black boxes", recognizing that this limits the codes that can be tested and the conclusions that can be drawn. One of the difficulties they recognize is, how do you account for codes failing to approximate the same solution, which sometimes happens at relaxed tolerances, or failing entirely for some reason? An alternative approach [3,34] based on measuring the behavior of codes on parametrized families of ODEs confronts many of the same difficulties and others, including step size strategies that do not behave smoothly with variation of the parameters. As a concrete illustration of the distinction between the approaches of testing the method in comparison with testing the code, a code based on Fehlberg's (4,5) pair showed up well in the tests [28]. However, RKF45 is much more efficient because it does local extrapolation and the code tested does not. Also, it controls the local error according to EPS and the codes tested in [28] were altered if necessary to use EPUS. Further, how can you compare two codes fairly when one determines an initial step size automatically and the other requires the user to supply a guess? The approach reported in [28] was to force the same starting step size on all its methods, an approach consistent with the stated goal of testing methods. The evaluation of [50] attempts to compare codes and recognizes that simplifying testing by providing the same "reasonably good" initial step size to all codes favors the less convenient codes that normally require the user to specify the initial step size. Two approaches to measuring efficiency have been used. In one the cost of the integration is plotted against the tolerance specified. In the other the cost is plotted against the accuracy of the numerical solution. The first is attractive to users because it corresponds exactly to the task set. It is implicit that the code actually solves the problem in some acceptable way. No credit is given for getting more accuracy than requested because the user has no way of knowing if any extra accuracy was achieved. The second measure is favored by those writing solvers because it exposes the effects of algorithmic decisions. It is important to appreciate the distinction, but the approaches are not incompatible and they are often mixed in practice. Now that issues of accuracy and the quality of error estimation are better understood, it is possible to analyze some important properties by studying the coefficients of terms in a Taylor series expansion of the local error. Similarly, stability is studied by computing and plotting stability regions. Rather than attempt to study the properties of a solver as a whole, particular properties are studied in a quantitative way. This is important because it separates characteristics of the method from the quality of the implementation. Equally important is that it is applicable to a broad class of problems. Numerical testing is always liable to a peculiar correlation of method, implementation, and test problem. Classic examples of formulas with such a correlation are the error estimate of the Merson pair behaving quite differently for the subclass comprising linear differential equations and the error estimates of high order pairs derived by Fehlberg behaving quite differently for the subclass comprising quadrature problems. When proceeding carefully, predictions of relative effectiveness made theoretically are confirmed in experiments with solvers. Checking the properties of formulas is relatively straightforward. This is not true of the goal of a smooth behavior of the error with respect to a change of tolerance because this behavior depends in a complicated way on the step size selection algorithm. Nonetheless, an understanding of the various parts of the algorithm and the implications for the behavior of the error can be used to predict the

304

L.E Shampine, I. Gladwell / Applied Numerical Mathematics 22 (1996) 293-308

behavior observed in practice. Observe that none of the popular codes goes as far as it might to secure smooth behavior; certainly none goes so far as to reject a step with an acceptably small local error estimate in order to take a larger step that is predicted will yield a local error closer to, but still smaller than, the tolerance. (True, RKSUITE does this on the first step, but only the first step.) A compromise is made with efficiency taking priority. This compromise weakens the link between local and global error, indeed so much so that it is quite possible for the global error to increase when the local error tolerance is reduced. This is one reason why reintegration using a smaller tolerance is not a truly reliable way to assess the global error.

8. Things to do In our experience the capabilities of software have been driven by the needs of users. There is always a question as to whether a particular need is of sufficiently general interest that it is worth providing in a general-purpose code. Here we mention a few capabilities that seem to us worth adding to RKSUITE. One already mentioned is a continuous extension for the (7,8) pair. Recently one of us had to consider the needs of the dynamic system simulation package SIMULINK in the course of developing the codes of the MATLAB ODE Suite [47]. Continuous simulation is an important application of ODE solvers that calls for capabilities not found in the popular general-purpose codes. Simulation packages often provide a solver based on an explicit RK formula as the default, and they must provide the capabilities somehow. For instance, event location is so important in the context that it is often provided. However, the implementations lack generality and lack the robustness of the quality stand-alone codes for initial value problems. Event location is not available in general-purpose RK codes, and it would be useful to add it to RKSUITE. The robust techniques in [44] are designed for this type of implementation, but adding event location in a user-friendly way is a non-trivial task. The techniques in [44] are robust partly because they are designed for a special, but frequently occurring, subclass of event functions 9 (see Section 6). A very careful attempt to deal with general g is that of Watts [57,58]. This work is in the context of Adams and BDF codes but could be transferred in a straightforward way to a RK code with a continuous extension. Another capability needed for continuous simulation is the solution of differential-algebraic equations (DAEs). We believe it is practical to provide for the solution of non-stiff, semi-explicit DAEs of index one in a general-purpose RK code like RKSUITE. Conservation laws express properties of the solution that are implied by the ODEs, properties that do not automatically hold for the numerical solution. It is shown in [41] how to impose them on numerical solutions computed with explicit RK methods. This is quite useful in certain contexts, and it is practical to make the capability available unobtrusively in general-purpose codes. It might provide an alternative to using symplectic formulas in some long-time integrations where a conservation law, say of energy or angular momentum, and symplecticness are both desirable but cannot be achieved simultaneously. Other matters being equal, the user might well be more concerned about conserving important physical quantities than symplecticness. The problem of parameter determination by means of initial value problems for ODEs is important in many areas. This is an optimization problem (normally constrained nonlinear least squares) and it requires the solution of an initial value problem to provide the data. It has long been recognised that the data must be smooth for the optimization to be efficient. Smoothness can be achieved by integrating the ODEs to a much greater accuracy than the accuracy tolerance in the optimization or by

L.E Shampine, L Gladwell /Applied Numerical Mathematics 22 (1996) 293-308

305

integrating using a "frozen" step size sequence, i.e., all integrations use the step size sequence chosen automatically in the first integration. The former approach is fine as long as the derivatives required for the optimization algorithm can be calculated to reasonable accuracy. Possibly, with the proliferation of symbolic and automatic differentiation packages available for calculating the requisite "variational" equations from the user's definition of the ODEs, this will become the method of choice because then the derivatives can be calculated as accurately as the data. However, if numerical differentiation is to be used by the optimization package to compute the derivatives it requires, it is unlikely that sufficient accuracy will be achieved in the ODE solution. An interesting illustration of this is found in [23]. In such circumstances the approach of freezing the step sequence is a viable alternative. It would be straightforward to add to RKSUITE the ability to write the "frozen" step size sequence and appropriate solution and derivative information to a sequential access file for retrieval for subsequent integrations and for evaluations of the data. The capability for RKSUITE to integrate with this "frozen" step size sequence could also be added easily. The difficulty here is that the first integration (the one written to a file) determines the step size sequence using a local error norm and the specified accuracy tolerance. Are we to monitor the local error norm in subsequent integrations using the "frozen" step size sequence and if so, how are we to define and to respond to estimated local errors that are "excessive"? Finally, consider the impact of new computing environments. The microcomputers that sit on our desks are entirely comparable in power to the mainframes we used when we first started solving ODEs numerically. The hand-held calculators of today were beyond our imagination of the time. These tools have made possible computer experiments in the teaching of ODEs, which has evoked specialized software for the purpose. These new computing environments are qualitatively different and the implications are still to be understood fully. An obvious difference is graphical display of results. This has implications for what codes compute and return to users. For instance, in the past few codes wrote results to an output file at every step because of issues of speed and portability. Now the equivalent is normal output when solutions are to be plotted routinely. There is much more emphasis in this context on software being easy to use. Reliance upon default values of parameters is not as dangerous as in the "batch processing" environments of the recent past, for in some ways these new environments resemble those of the early days when users monitored computations very closely for indications of trouble. In this context when a solution does not appear to behave as expected, adjustments are made and a new solution computed immediately. Still, there is an obvious need for codes to monitor computations and to assist users who do not know what to expect by providing good advice. This is especially true because of another implication of graphical output, namely relaxed tolerances. When only graphical output is desired, it is natural to specify nominal accuracies, either directly or by default in the code. This may be imprudent because at graphical accuracies the asymptotic approximations that underlie many algorithms in the codes may be poor. More research will be required to meet the challenge of reliable computation at nominal accuracies. In the case of hand-held calculators there remains a premium on memory, so short programs implementing explicit RK formula pairs with few stages are the methods of choice. There is an interesting analogy here with the RK formulas used in the massive MOL calculations in three space dimensions where memory constraints continue to dominate to the extent that few stages can be kept simultaneously. The same constraints influenced the early teaching packages, but today memory and speed of desktop computers are so great that general-purpose codes can be used. (Will this be the case soon for the calculators or the MOL calculations?) For instance, a multmedia teaching tool being developed by C.ODE-E (Consortium for Ordinary Differential Equations Experiments) is based on RKSUITE and

306

L.E Shampine, 1. Gladwell / Applied Numerical Mathematics 22 (1996) 293-308

VODE [9] with a thin interface to make the underlying codes solve problems in the same sense and "look" the same to the graphical user interface. This raises another point, the user interface. The new environments call for very different ways of specifying problems and interacting with solvers. In particular, teaching packages solve problems of restricted type and this simplifies development of easy-to-use schemes for defining problems and displaying results. In the future this might also influence the choice of formulas.

Acknowledgements Much of the work described in this paper was supported by NATO through grants 898/83 and CRG920037.

References [1] P. Bogacki and L.E Shampine, A 3(2) pair of Runge-Kutta formulas, Appl. Math. Lett. 2 (1989) 1-9. [2] P. Bogacki and L.E Shampine, An efficient Runge-Kutta (4,5) pair, Report 89-20, Mathematics Department, SMU, Dallas (1989), Comput. Math. Appl., to appear. [3] E Biddle, A technique for evaluating the performance of initial value solvers for non-stiff ordinary differential equations, M.Sc. Thesis, Department of Mathematics, University of Manchester, UK (1977). [4] R.W. Brankin and I. Gladwell, Using shape preserving local interpolants for plotting solutions of ordinary differential equations, IMA J. Numer. Anal 9 (]989) 555-566. [5] R.W. Brankin, I. Gladwell, J.R. Dormand, EJ. Prince and W.L. Seward, ALGORITHM 670: A RungeKutta-Nystr/Sm code, ACM Trans. Math. Software 15 (1989) 31-40. [6] R.W. Brankin, I. Gladwell and L.E Shampine, Starting BDF and Adams codes at optimal order, J. Comput. Appl. Math. 21 (1988) 357-368. [7] R.W. Brankin, I. Gladwell and L.E Shampine, RKSUITE: a suite of Runge-Kutta codes for the initial value problem for ODEs, Softreport 91-1, Mathematics Department, SMU, Dallas (1991). [8] R.W. Brankin, I. Gladwell and L.F. Shampine, RKSUITE: A suite of explicit Runge-Kutta codes, in: R. Agarwal, ed., Contributions in Numerical Mathematics, World Scientific Series in Applicable Analysis 2 (World Scientific, Singapore, 1993) 85-98. [9] EN. Brown, G.D. Byrne and A.C. Hindmarsh, VODE: a variable-coefficient ODE solver, SIAM J. Sci. Statist. Comput. 10 (1989) 1038-1051. [10] M.H. Carpenter, D. Gottlieb, S. Abarbanel and W.S. Don, The theoretical accuracy of Runge-Kutta discretizations for the initial boundary value problem: A study of the boundary error, SIAM J. Sci. Comput., to appear. [I1] J.R. Dormand and P.J. Prince, A family of embedded Runge-Kutta formulae, J. Comput. Appl. Math. 6 (1980) 19-26. [12] J.R. Dormand and EJ. Prince, Global error estimation with Runge-Kutta methods, IMA J. Numer. Anal. 4 (1984) 169-184. [13] J.R. Dormand and EJ. Prince, Global error estimation with Runge-Kutta methods II, 1MA J. Numer. Anal. 5 (1985) 481-497. [14] W.H. Enright, Analysis of error control strategies for continuous Runge-Kutta methods, SIAM J. Numer. Anal. 26 (1989) 588-599. [15] E. Fehlberg, Klassische Runge-Kutta-Formeln vierter und niedrigerer Ordnung mit Schrittweiten-Kontrolle und ihre Anwendung auf W~irmeleitungsprobleme, Computing 6 (1970) 61-71.

L.E Shampine, L Gladwell / Applied Numerical Mathematics 22 (1996) 293-308

307

[16] C.W. Gear, Runge-Kutta starters for multistep methods, ACM Trans. Math. Software 6 (1980) 263-279. [17] I. Gladwell, Initial value routines in the NAG library, ACM Trans. Math. Software 5 (1979) 386-400. [18] I. Gladwell, J.A.I. Craigie and C.R. Crowther, Testing initial value routines as black boxes, Numerical Analysis Report 34, Department of Mathematics, University of Manchester (1979). [19] I. Gladwell, L.E Shampine and R.W. Brankin, Automatic selection of the initial step size for an ODE solver, J. Comput. Appl. Math. 18 (1987) 175-192. [20] K. Gustafsson, Control theoretic techniques for stepsize selection in explicit Runge-Kutta methods, ACM Trans. Math. Software 17 (1991)533-554. [21] K. Gustafsson, M. Lundh and G. S6derlind, A PI stepsize control for the numerical solution of ordinary differential equations, BIT 18 (1988) 270-287. [22] E. Hairer and D. Stoffer, Reversible long-time integration with variable step sizes, SIAM J. Sci. Comput., to appear. [123] E. Hairer, S.R N6rsett and G. Wanner, Solving Ordinary Differential Equations L Nonstiff Problems (Springer, Berlin, 1987). [24] E. Hairer and G. Wanner, Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems (Springer, Berlin, 1991). [25] D.J. Higham, Robust defect control with Runge-Kutta schemes, SIAM J. Numer. Anal. 26 (1989) 1175-1183. [26] M.K. Horn, Fourth and fifth-order scaled Runge-Kutta algorithms for treating dense output, SIAM J. Numer. Anal. 20 (1983) 558-568. [27] F.Q. Hu, M.Y. Hussaini and J. Manthey, Low dissipation and dispersion Runge-Kutta schemes for computational acoustics, ICASE Report 94-102, NASA Langley Research Center (1994). [28] T.E. Hull, W.H. Enright, B.M. Fellen and A.E. Sedgwick, Comparing numerical methods for ordinary differential equations, SIAM J. Numer. Anal. 9 (1972) 603-637. [29] T.E. Hull, W.H. Enright and K.R. Jackson, User's guide for DVERK--a subroutine for solving non-stiff ODEs, Technical Report 100, Department of Computer Science, University of Toronto, Ontario (1976). [30] T.E. Hull, W.H. Enright and K.R. Jackson, Runge-Kutta research at Toronto, Appl. Numer. Math. 22 (1996) 225-236. [31] K.R. Jackson, W.H. Enright and T.E. Hull, A theoretical criterion for comparing Runge-Kutta formulas, SlAM J. Numer. Anal. 15 (1978) 618-641. [32] G.-S. Jiang and C.-W. Shu, Efficient implementation of weighted ENO schemes, ICASE Report 95-73, NASA Langley Research Center (1995). [33] R.H. Merson, An operational method for the study of integration processes, in: Proc. Sympos. Data Ptwcessing, Weapons Research Establishment, Salisbury, Australia (1957) 110-1-110-25. [34] M.L. Parker, Performance profile testing of ordinary differential equation routines, M.Sc. Thesis, Department of Mathematics, University of Manchester, UK (1979). [35] J.C. Polking, MATLABManual for Ordinary Differential Equations (Prentice-Hall, Englewood Cliffs, NJ, 1995). [36] EJ. Prince and J.R. Dormand, High order embedded Runge-Kutta formulae, J. Comput. Appl. Math. 7 (1981) 67-75. [37] L.F. Shampine, Local extrapolation in the solution of ordinary differential equations, Math. Comp. 27 (1973) 91-97. [38] L.F. Shampine, Stiffness and non-stiff differential equation solvers, in: L. Collatz et al., eds., Numerische Behandlung von Differentialgleichungen, International Series of Numerical Mathematics 27 (Birkhauser, Basel, 1975) 287-301. [39] L.F. Shampine, The step sizes used by one-step codes for ODEs, Appl. Numer. Math. 1 (1985) 95-106. [40] L.F. Shampine, Interpolation for Runge-Kutta methods, SlAM J. Numer. Anal. 22 (1985) 1014-1027.

308

L.E Shampine, 1. Gladwell/ Applied Numerical Mathematics 22 (1996) 293-308

[41] L.E Shampine, Conservation laws and the numerical solution of ODEs, Comput. Math. Appl. 12B (1986) 1287-1296. [42] L.F. Shampine, Diagnosing stiffness for Runge-Kutta methods, SlAM J. Sci. Statist. Comput. 12 (1991) 260-272. [43] L.F. Shampine, Numerical Solution of Ordinary Differential Equations (Chapman and Hall, New York, 1994). [44] L.E Shampine, I. Gladwell and R.W. Brankin, Reliable solution of event location problems for ODEs, ACM Trans. Math. Software 17 (1991) 11-25. [45] L.E Shampine and M.K. Gordon, Some numerical experiments with DIFSUB, SIGNUM Newsletter 7 (1972) 24-26. [46] L.E Shampine and K.L. Hiebert, Detecting stiffness with the Fehlberg (4,5) formulas, Comput. Math. Appl. 3 (1977) 41-46. [47] L.E Shampine and M.W. Reichelt, The MATLABODE suite, SlAM J. Sci. Comput., to appear. [48] L.F. Shampine and H.A. Watts, Global error estimation for ordinary differential equations and Algorithm 504, GERK: global error estimation for ordinary differential equations, ACM Trans. Math. Software 2 (1976) 172-186 and 200-203. [49] L.F. Shampine and H.A. Watts, The art of writing a Runge-Kutta code I, in: J.R. Rice, ed., Mathematical Software III (Academic Press, New York, 1977) 257-275; and The art of writing a Runge-Kutta code II, Appl. Math. Comput. 5 (1979) 93-121. [50] L.E Shampine, H.A. Watts and S.M. Davenport, Solving nonstiff ordinary differential equations--the state of the art, SlAM Rev. 18 (1976) 376-411. [51 ] P.W. Sharp and E. Smart, Explicit Runge-Kutta pairs with one more derivative evaluation than the minimum, SlAM J. Sci. Comput. 14 (1993) 338-348. [52] P.W. Sharp and J.H. Verner, Explicit Runge-Kutta 4-5 pairs with interpolants, Math. Preprint 1995-03, Queen's University, Kingston, Ontario (1995). [53] C.-W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes I, J. Comput. Phys 77 (1989) 439-471. [54] J.H. Verner, Explicit Runge-Kutta methods with estimates of the local truncation error, SlAM J. Numer. Anal. 15 (1978) 772-790. [55] J.H. Verner, Private communication (1992). [56] H.A. Watts, Step size control in ordinary differential equation solvers, Trans. Soc. Computer Simulation 1 (1984) 15-25. [57] H.A. Watts, RDEAM--an Adams ODE code with root solving capability, Report SAND85-1595, Sandia National Laboratories, Albuquerque, NM (1985). [58] H.A. Watts, Backward differentiation formulas revisited: improvements in DEBDF and a new root solving code RDEBD, Report SAND85-2676, Sandia National Laboratories, Albuquerque, NM (1986). [59] J.A. Zonneveld, Automatic Numerical Integration (Mathematisch Centrum, Amsterdam, 1964).