An improved class of generalized Runge–Kutta methods for stiff problems. Part I: The scalar case
Introduction
During the last decades there has been a considerable amount of research on methods for numerical integration of stiff systems of ODEs, usually looking for better stability properties. Nearly all such methods are implicit in character.
The most widely used algorithms are those based on linear multistep formulas like the BDF methods (see e.g. [1]), because they are very efficient for general stiff problems. However, from a result of Dahlquist it is known that no linear multistep method of order greater than 2 can be A-stable [2], and so these formulas are not suited to some stiff problems (for example, those with Jacobians whose eigenvalues have large imaginary parts).
Implicit Runge–Kutta formulas [3], [4] have been widely used because of their excellent stability properties (such as A-stability, L-stability and B-stability), but the need for solving non-linear algebraic equations at each step makes these formulas generally too costly when considering some huge systems of ODEs.
To reduce the amount of computational effort required to solve the non-linear equations (by Newton-type iterations) when integrating with a fully implicit Runge–Kutta method, some classes of Runge–Kutta formulas have been developed. With the class of diagonally implicit Runge–Kutta (DIRK) methods (also called semi-implicit or semi-explicit Runge–Kutta methods) the algebraic cost of the LU factorization is reduced. By considering the class of singly diagonally implicit Runge–Kutta (SDIRK) methods and also the class of singly implicit Runge–Kutta (SIRK) methods, it is possible to reduce even more the algebraic cost of the LU factorization (see e.g. [3], [5], [6] for more details).
Many other attempts have been made in order to reduce the computation cost per step by considering linearly implicit methods, in this way eliminating the need for solving non-linear systems which, as pointed before, usually are solved by Newton-type iteration (and hence require additional function evaluations for every iteration at every step). Such formulas have the computational advantage that it is necessary to solve only linear systems of algebraic equations at each step.
Among the many different RK-like methods of this type we have the Rosenbrock methods [7] and the ROW-methods (also called Rosenbrock–Wanner methods and modified Rosenbrock methods) [8], [9], [10]. These formulas, however, require the exact Jacobian at every step. Therefore the computations are costly when the Jacobian matrix is expensive to evaluate. For this reason, extensions of Rosenbrock methods have been considered in which the exact Jacobian is fixed for some number of steps so that the computation cost is reduced (see e.g. [11], [12], [13], [14]). Moreover, Rosenbrock-type methods in which the exact Jacobian is no longer needed have been considered. The so-called W-methods [15], the MROW-methods [16] and the generalized Runge–Kutta methods [17] (see also [18]) fall into this class. For an excellent survey of some of these methods the reader is referred to [5].
In our recent papers [19], [20], examples are shown of explicit and linearly implicit two-stage methods of order 3 for the numerical integration of scalar autonomous ODEs, which do not require Jacobian evaluations, some of them being as well A-stable and L-stable. Some comparisons with Runge–Kutta methods as well as numerical experiments are also reported. In [21] we describe how to construct from a given function R a one-parameter family of explicit (or linearly implicit) two-stage methods, having R as the associated stability function, and illustrate this fact by obtaining a two-stage third-order formula whose associated stability function is given by ez.
Our first aim in the present paper is to introduce the general form of the new explicit p-stage methods for the numerical integration of scalar autonomous ODEs. We begin studying these methods for scalar problems because, as we will see in a second part, most of the work can be easily extended in order to integrate separated systems of ODEs. These methods can be seen as a generalization of the explicit Runge–Kutta methods providing better order and stability results with the same number of stages. In fact, from Butcher's theory we know that an p-stage explicit Runge–Kutta method cannot have order greater than p. Moreover, the stability function of such methods is a polynomial, and so none of them is A-stable. We will show that it is possible to obtain A-stable explicit formulas for scalar autonomous problems of order 3 and 5, with only two and three stages, respectively, from our class of methods. By a further generalization of our schemes, we will show in the second part that it is possible to obtain A-stable linearly implicit formulas for some non-autonomous scalar ODEs and systems, which do not require Jacobian evaluations. For example, in [22] we present a two-stage third-order method for separated systems of ODEs being L-stable, as well as preliminary results using our method to integrate systems that arise when solving some non-linear parabolic PDEs by the method of lines approach.
Finally we illustrate the efficiency of those schemes by carrying out some numerical experiments.
Section snippets
The new family of methods
In this first part, we will restrict our attention to the scalar autonomous initial value problemFor this problem, let us consider the family of explicit p-stage methods defined bywhere the stages are given byand for each i with 2⩽i⩽p+1, Fi is any homogeneous function of degree 1, that isholds for each and (x1,x2,…,xi−1) in a subset
A useful formulation
We introduce the termswhere the stages ki are given by (3). From considerations that will be clear later, we take si=0 when k1=0 in (6).
It is a simple task to show that si=O(h), and we will exploit this property in order to simplify our study. Moreover, the terms si can be seen as approximations to cihfy(yn). In fact si=cihfy(yn)+O(h2) (here we assume that f has a sufficient number of bounded derivatives), and so we can obtain approximations to the Jacobian fy(yn) by
Consistency and order of the methods
In what follows we will assume that is Lipschitz in , i.e. there exists a Lipschitz constant L such thatWith the previous assumptions it is well known that for any , there exists a unique solution y(x) of problem (1) (throughout any interval [x0,b]), where y(x) is continuous and differentiable.
We will investigate the consistency of any p-stage method given by (7). When doing so, we will assume the existence (and continuity) of y′(x) in [x0,b], but
Methods of polynomial type
Now, for every fixed p, we restrict our attention to the family of methods given by , , where now all the are assumed to be homogeneous functions (of degree 1) of the special formwith all ri being non-negative integers.
The above family of methods, still contains all explicit p-stage Runge–Kutta methods. In fact, taking ri=1 in (12) we get from , all explicit p-stage Runge–Kutta methods.
In what follows, we will
Two-stage methods of polynomial type: attainable order
To show that we can obtain explicit two-stage methods of polynomial type from , (together with , ) of order 3, it is enough to consider Taylor's expansion of the associated local truncation error. Note that when doing so, only the parameters c2, c3 and ai with 1⩽i⩽2 will appear in the order conditions, that is, it suffices to take in (15) r3=2 (the other parameters can be arbitrarily chosen). This easily follows from the fact that s2=O(h), and therefore s2k=O(hk) for any given . We obtain
Three-stage methods of polynomial type: attainable order
When considering three-stage methods of polynomial type from , together with , , it is possible to get a family of fifth-order formulas, depending on many free parameters. As in the two-stage case, only some of this parameters will appear in the order conditions, due to the fact that both s2 and s3 are O(h). However, the resulting order conditions are still too cumbersome; hence we define a new term to make our study easier. Note that from our previous considerations it is clear that
Methods of rational type
In the last three sections we have only considered methods of polynomial type. Now we will consider methods of rational type, that is, methods given by , , where now all the Fi with 2⩽i⩽p+1 are supposed to be homogeneous functions (of degree 1) of rational type. More precisely, we will consider functions Fi given in terms of the quotient of two homogeneous polynomials and with degrees ri+1 and ri, respectively, for some non-negative integer ri, that is
Two-stage methods of rational type
To obtain all the explicit two-stage methods (of rational type) of order 3, it suffices to note that from , we have that G2=c2=2/3 and alsomust hold. Obviously it is enough to consider in (46), obtaining the following conditions:If we want to minimize the principal part of the local truncation error, all we need is to take in (46) and expand to higher order the second term,
Three-stage methods of rational type
Now, following our notations in Section 7, we can describe any three-stage method of rational type bywith , and where as usually s2 and s3 are given through (13) in terms of the stagesNow the rational functions are given byand it is easily seen that when looking for fifth-order formulas, it is enough to consider the
Linear stability properties of the methods
Now we are going to study the linear stability properties of the methods. When we apply a p-stage method (7) of our family to the scalar test equationwe getwhere R(z) is the associated stability function with z=hλ. Moreover, from , , , we obtain recursivelyfrom which we obtain the following expression for the stability function:where the si are given in terms of z by the
Three-stage methods being A-stable
From the above section it is clear that the stability function of a method of polynomial type is always a polynomial function. Therefore it is not possible to obtain formulas of polynomial type with good linear stability properties such as A-stability or L-stability. However, it is also clear that we can obtain methods of rational type whose associated stability function is a rational function. Moreover, we can obtain A-stable and L-stable methods of rational type, without losing the highest
Numerical experiments (the scalar case)
In order to show the behaviour as h→0 for the methods explained in the last section, we will consider the following simple problem (taken from [3], p. 134)for which the solution isWith fixed step size h=2−n for various n=1,2,3,…,9 over 2n steps, the value of y(1) was computed using our methods and two Runge–Kutta methods. The magnitude of the error E for different h and for each of these methods is shown in Fig. 1 in logarithmic scale. The
Conclusions
In this first part we have introduced and studied a new family of methods for the numerical integration of scalar autonomous ODEs. This work will help us in a second part because most of the results and properties can be easily extended to system problems. These new methods seem quite promising, for instance in the context of solving some non-linear parabolic equations (by the method of lines). We will investigate this and other questions in part II: `The separated system case'.
References (22)
Avoiding the exactness of the Jacobian matrix in Rosenbrock formulae
Comput. Math. Appl.
(1990)Numerical Methods for Ordinary Differential Systems. The Initial Value Problem
(1991)A special stability problem for linear multistep methods
BIT
(1963)The Numerical Analysis of Ordinary Differential Equations: Runge–Kutta and General Linear Methods
(1987)Implicit Runge–Kutta processes
Math. Comp.
(1964)- et al.
Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems
(1996) Diagonally implicit Runge–Kutta methods for stiff O.D.E.'s
SIAM J. Numer. Anal.
(1977)Some general implicit processes for the numerical solution of differential equations
Comput. J.
(1963)- et al.
Order conditions for Rosenbrock type methods
Numer. Math.
(1979) - et al.
Generalized Runge–Kutta methods of order four with stepsize control for stiff ordinary differential equations
Numer. Math.
(1979)
A study of Rosenbrock-type methods of high order
Numer. Math.
Cited by (19)
Time discretization of the point kinetic equations using matrix exponential method and First-Order Hold
2013, Annals of Nuclear EnergyCitation Excerpt :Li et al., developed a numerical integral method for solving the point kinetics equations using the better bias function (BBF) to approximate the neutron density in on-time step integrations (Nahla, 2011). Also, Alvarez proposed a new family of p-stage methods for the numerical integration of some scalar equations and systems of ODEs (Alvarez and Rojo, 2002). The above given methods are unconventional numerical techniques and require a “small” time step in order to be deemed accurate.
Convergence on error correction methods for solving initial value problems
2012, Journal of Computational and Applied MathematicsGeneralized Runge-Kutta method for two- and three-dimensional space-time diffusion equations with a variable time step
2008, Annals of Nuclear EnergyAn improved class of generalized Runge-Kutta methods for stiff problems. Part II: The separated system case
2004, Applied Mathematics and ComputationCitation Excerpt :From the preceding comments it is now clear that matrices Si can be seen as a generalization (together with a minor modification) of the terms si considered in the scalar case. In fact, it is not difficult to see that the family of methods introduced in [7] for scalar problems can also be formulated in this manner, but we will not do that here. From our last observation we can conclude that our GRK-methods are similar to the Rosenbrock methods [19] and other related formulae such as the W-methods [22], the MROW-methods [24] and the generalized Runge–Kutta methods [23] for which the exact Jacobian matrix must not be computed when methods are implemented.
An improved class of generalized Runge-Kutta-Nyström methods for special second-order differential equations
2004, Communications in Nonlinear Science and Numerical SimulationOn pade' approximations to the exponential function and application to the point kinetics equations
2004, Progress in Nuclear Energy
- 1
This work was partially supported by Consejerı́a de Educación y Cultura de la Junta de Castilla y León and Unión Europea (Fondo Social Europeo) under proyect VA19/00B and by Junta de Castilla y León under proyect VA11/99.
- 2
This work was partially supported by Consejerı́a de Educación y Cultura de la Junta de Castilla y León and Unión Europea (Fondo Social Europeo) under proyect VA19/00B and by Spanish DGES under project PB98-0359.