J.R. Cash

The integration of stiff initial value problems in ODEs using modified extended backward differentiation formulae

Abstract A class of modified extended backward differentiation formulae (MEBDF) suitable for the approximate numerical integration of first order systems of stiff ordinary differential equations is introduced. The computational aspects of this new class of formulae are examined in detail. In particular, algorithms for varying both order and stepsize are given and this leads to a variable step/variable order process using highly stable formulae of order 2–8. Extensive numerical results for the well known DETEST set are given and on this basis a comparison is made between a code incorporating MEBDF, a code due to Hindmarsh (based on conventional backward differentiation formulae) and a code due to Skeel and Kong (based on blended linear multistep methods). It is shown that on this test set the MEBDF code is very reliable and is competitive with the other two codes on a significant class of problems. A computer listing of the MEBDF code used to obtain the results presented in this paper is available from the author.

Lobatto Deferred Correction for Stiff Two-Point Boundary Value Problems

Abstract An iterated deferred correction algorithm based on Lobatto Runge-Kutta formulae is developed for the efficient numerical solution of nonlinear stiff two-point boundary value problems. An analysis of the stability properties of general deferred correction schemes which are based on implicit Runge-Kutta methods is given and results which are analogous to those obtained for initial value problems are derived. A revised definition of symmetry is presented and this ensures that each deferred correction produces an optimal increase in order. Finally, some numerical results are given to demonstrate the superior performance of Lobatto formulae compared with mono-implicit formulae on stiff two-point boundary value problems.

A variable order deferred correction algorithm for the numerical solution of nonlinear two point boundary value problems

Abstract We develop two variable order deferred correction algorithms, based on finite difference methods, for the approximate numerical solution of nonlinear two point boundary value problems. The deferred corrections provide a powerful local error estimator and this estimate can be used to refine the mesh, if need be, in order to obtain a prescribed tolerance. Mesh refinement is particularly straightforward since our algorithms are truly one step in nature. The approach described is particularly appropriate for problems where only low, or “engineering”, precision is required and for very large systems of equations where storage space is an extremely important consideration. Some numerical results are given and a comparison with the method developed by Lentini and Pereyra and with Richardson extrapolation is made.

The MOL solution of time dependent partial differential equations

Abstract A modified BDF scheme is proposed for the numerical integration of the ordinary differential equations that arise in the method of lines solution of time dependent partial differential equations. It is to be expected that this new approach will be superior to the use of BDF for at least three important classes of problems, namely for advection diffusion problems where advection dominates and here we expect very large gains in efficiency, for problems where high accuracy is required and for problems where function evaluations are very expensive. Some analytic and numerical results are given to confirm these expectations.

An MEBDF package for the numerical solution of large sparse systems of stiff initial value problems

Abstract An efficient algorithm for the numerical integration of large sparse systems of stiff initial value ordinary differential equations (ODEs) and differential-algebraic equations (DAEs) is described. The algorithm is constructed by embedding a standard sparse linear algebraic equation solver into a suitably modified MEBDF code. An important practical application of this algorithm is in the numerical solution of time dependent partial differential equations (PDEs), particularly in two or more space dimensions, using the method of lines (MOL). A code based on this algorithm is illustrated by application to several problems of practical interest and its performance is compared to that of the standard code LSODES.

Adaptive Runge-Kutta methods for nonlinear two-point boundary value problems with mild boundary layers

Abstract A variable mesh deferred correction algorithm based on implicit Runge-Kutta formulae is described for the approximate numerical solution of nonlinear two-point boundary value problems. A strategy for automatically choosing the variable mesh spacing is described and this seeks to equidistribute an approximation to the global truncation error of the Runge-Kutta formula. The facility of being able to use a variable mesh extends considerably the range of applicability of Runge-Kutta methods and in particular allows the possibility of solving rather more difficult problems, such as those with mild boundary layers or turning points, in an efficient manner. An extensive set of numerical results is given to illustrate the various algorithms described and a comparison is made with the deferred correction code of Lentini and Pereyra appearing in the NAG library.

Diagonally implicit runge-kutta formulae for the numerical integration of nonlinear two-point boundary value problems

Abstract Fully implicit Runge-Kutta formulae, based on interpolatory quadrature schemes, for the approximate numerical solution of nonlinear two-point boundary value problems have been investigated by Weiss. Such formulae are not used very often in practice, however, because they generally require such a large computational effort as to make them uncompetitive with, for example, integration schemes based on the trapezoidal or implicit mid-point rules. In this paper we consider an alternative class of Runge-Kutta formulae, namely diagonally implicit Runge-Kutta (DIRK) formulae, which can be implemented more efficiently than the fully implicit formulae considered by Weiss. We also consider how these DIRK formulae can be implemented in a defect or deferred correction framework and we give some numerical results to illustrate the algorithms derived. One particular formula belonging to the DIRK class is the implicit mid-point rule. In this paper we derive an efficient implementation of this formula which is applicable when the given boundary conditions are non-separated.

A comparison of some codes for the stiff oscillatory problem

Abstract In 1984, Gaffney [1] carried out a survey of codes for the numerical solution of the stiff oscillatory problem. This was prompted by a desire to derive an efficient method of lines algorithm for the reduced resistive MHD equations. His main conclusion was that none of the codes available at that time was entirely satisfactory for dealing with this problem. In the present paper, we outline the considerable progress that has been made since Gaffneys survey and show that there now exist much more satisfactory codes for dealing with these problems. Two of the three FORTRAN codes that we recommend are both freely available from NETLIB, as well as from certain web pages, and so it is straightforward to reproduce all the numerical results given in the present paper.

Block embedded explicit runge-kutta methods

Abstract Embedded explicit Runge—Kutta formulae are amongst the most popular methods currently in use for the approximate numerical integration of nonstiff systems of ordinary differential equations. In particular the fourth-order methods RKF45 of Shampine and Watts and the Runge—Kutta—Merson code of the NAG Library have been particularly popular for some time. In this paper we consider the derivation of a block embedded explicit Runge—Kutta (BERK) formula of order 4. BERK formulae have all the characteristics of standard explicit Runge—Kutta formulae except that they are no longer single step in nature in the sense that a p th order BERK formula produces p th order approximations to the solution at several step points instead of at one point only. The results of some fairly extensive numerical testing are presented and these indicate that the new formula is very competitive, both in terms of efficiency and reliability, with the codes R—K Merson of NAG and RKF45. Also considered in detail are the problems of computing solutions at “off-step” points and of efficiently computing low accuracy solutions. In particular we show that BERK formulae allow intermediate solutions at “off-step” points to be computed with relatively little additional computational effort. This makes BERK formulae particularly attractive for problems where output is requested at many off-step points since it is for this class of problems that conventional Runge—Kutta formulae can become very inefficient. Also considered is the question of computing low accuracy solutions. It is shown that, as is to be expected, such problems are solved more efficiently using low order BERK formulae. In particular, BERK formulae of orders 1 and 3 are derived and their performance is compared with that of the fourth-order BERK formula on a large set of test problems.

Stability concepts in the numerical solution of difference and differential equations

Abstract Several algorithms have been proposed for the stable numerical computation of non-dominant solutions of linear difference equations. The unifying approach of these algorithms is that they replace the task of solving the given initial value problem, by that of solving a boundary value problem, which is in some sense, equivalent. A similar task is involved in the numerical solution of initial value problems for ordinary differential equations using linear multistep methods. However, in this case, it is always assumed that the numerical solution will be found by forward recurrence, and this assumption imposes severe restrictions on the numerical methods that can be used. By adapting some of the ideas that have been proposed for linear recurrence relations, new algorithms are derived for the numerical solution of initial value problems, and it is shown that these new methods can be competitive with the best methods currently available.