An improved class of generalized Runge–Kutta methods for stiff problems. Part I: The scalar case

doi:10.1016/S0096-3003(01)00115-1

Applied Mathematics and Computation

Volume 130, Issues 2–3, 15 August 2002, Pages 537-560

https://doi.org/10.1016/S0096-3003(01)00115-1 Get rights and content

Abstract

A new family of p-stage methods for the numerical integration of some scalar equations and systems of ODEs is proposed. These methods can be seen as a generalization of the explicit p-stage Runge–Kutta ones, while providing better order and stability results. We will show in this first part that, at the cost of losing linearity in the formulas, it is possible to obtain explicit A-stable and L-stable methods for the numerical integration of scalar autonomous ODEs. Scalar autonomous ODEs are of very little interest in current applications. However, be begin studying this kind of problems because most of the work can be easily extended to a more general situation. In fact, we will show in a second part (entitled `The separated system case'), that it is possible to generalize our methods so that they can be applied to some non-autonomous scalar ODEs and systems. We will obtain linearly implicit L-stable methods which do not require Jacobian evaluations. In both parts, some numerical examples are discussed in order to show the good performance of the new schemes.

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 e^z.

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 problem $y^{′} (x)=f(y(x)), y(x_{0})=y_{0} .$ For this problem, let us consider the family of explicit p-stage methods defined by $y_{n+1} =y_{n} +hF_{p+1} (k_{1},k_{2},…,k_{p}),$ where the stages are given by $k_{1} = f(y_{n}),$ $k_{2} = f(y_{n} +hF_{2} (k_{1})),$ $k_{3} = f(y_{n} +hF_{3} (k_{1},k_{2})),$ $⋮$ $k_{p} = f(y_{n} +hF_{p} (k_{1},k_{2},…,k_{p−1})),$ and for each i with 2⩽i⩽p+1, F_i is any homogeneous function of degree 1, that is $F_{i} (αx_{1},αx_{2},…,αx_{i−1})=αF_{i} (x_{1},x_{2},…,x_{i−1}),$ holds for each $α∈ R$ and (x₁,x₂,…,x_i−1) in a subset

A useful formulation

We introduce the terms $s_{i} = k_{i} −k_{1} k_{1} = k_{i} k_{1} −1, 2⩽i⩽p,$ where the stages k_i are given by (3). From considerations that will be clear later, we take s_i=0 when k₁=0 in (6).

It is a simple task to show that s_i=O(h), and we will exploit this property in order to simplify our study. Moreover, the terms s_i can be seen as approximations to c_ihf_y(y_n). In fact s_i=c_ihf_y(y_n)+O(h²) (here we assume that f has a sufficient number of bounded derivatives), and so we can obtain approximations to the Jacobian f_y(y_n) by

Consistency and order of the methods

In what follows we will assume that $f: R → R$ is Lipschitz in $R$ , i.e. there exists a Lipschitz constant L such that $|f(y)−f(y^{∗})|⩽L|y−y^{∗} | for every y,y^{∗} ∈ R .$ With the previous assumptions it is well known that for any $y_{0} ∈ R$ , there exists a unique solution y(x) of problem (1) (throughout any interval [x₀,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 [x₀,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 $F_{i} (2⩽i⩽p+1)$ are assumed to be homogeneous functions (of degree 1) of the special form $F_{i} (x_{1},x_{2},…,x_{i−1})=∑_{j_{2}+j_{3}+⋯+j_{i−1}=0}^{r_{i}} A_{j_{2}j_{3}⋯j_{i−1}} x_{1} x_{2} x_{1}^{j_{2}} x_{3} x_{1}^{j_{3}} ⋯ x_{i−1} x_{1}^{j_{i−1}}$ with all r_i being non-negative integers.

The above family of methods, still contains all explicit p-stage Runge–Kutta methods. In fact, taking r_i=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 c₂, c₃ and a_i with 1⩽i⩽2 will appear in the order conditions, that is, it suffices to take in (15) r₃=2 (the other parameters can be arbitrarily chosen). This easily follows from the fact that s₂=O(h), and therefore s₂^k=O(h^k) for any given $k∈ N$ . 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 s₂ and s₃ are O(h). However, the resulting order conditions are still too cumbersome; hence we define a new term $s ̃_{3} =s_{3} −s_{2}$ 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 F_i with 2⩽i⩽p+1 are supposed to be homogeneous functions (of degree 1) of rational type. More precisely, we will consider functions F_i given in terms of the quotient of two homogeneous polynomials $N ̃_{i} (x_{1},x_{2},…,x_{i−1})$ and $D ̃_{i} (x_{1},x_{2},…,x_{i−1})$ with degrees r_i+1 and r_i, respectively, for some non-negative integer r_i, that is $F_{i} (x_{1},x_{2}$

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 G₂=c₂=2/3 and also $G_{3} (s_{2})=c_{3} 1+∑_{i=1}^{n_{3}^{∗}} n_{i} s_{2}^{i} 1+∑_{i=1}^{d_{3}^{∗}} d_{i} s_{2}^{i} =1+ 12 s_{2} + 16 s_{2}^{2} + O (s_{2}^{3})$ must hold. Obviously it is enough to consider $n_{3}^{∗} =d_{3}^{∗} =2$ in (46), obtaining the following conditions: $c_{2} =2/3,$ $c_{3} =1,$ $n_{1} =1/2+d_{1},$ $n_{2} =1/6+(1/2)d_{1} +d_{2} .$ If we want to minimize the principal part of the local truncation error, all we need is to take $n_{3}^{∗} =d_{3}^{∗} =3$ 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 by $y_{n+1} =y_{n} +hk_{1} G ̃_{4} (s_{2}, s ̃_{3})$ with $s ̃_{3} =s_{3} −s_{2}$ , and where as usually s₂ and s₃ are given through (13) in terms of the stages $k_{1} =f(y_{n}), k_{2} =f y_{n} +hk_{1} G_{2}, k_{3} =f y_{n} +hk_{1} G_{3} (s_{2}),$ Now the rational functions are given by $G_{2} =c_{2},$ $G_{3} (s_{2})=c_{3} 1+∑_{i=1}^{n_{3}^{∗}} n_{i} s_{2}^{i} 1+∑_{i=1}^{d_{3}^{∗}} d_{i} s_{2}^{i},$ $G ̃_{4} (s_{2}, s ̃_{3})=c_{4} 1+∑_{i+2j=1}^{ñ_{4}^{∗}} n_{ij} s_{2}^{i} s ̃_{3}^{j} 1+∑_{i+2j=1}^{d̃_{4}^{∗}} d_{ij} s_{2}^{i} s ̃_{3}^{j},$ and 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 equation $y^{′} =λy, λ∈ C,$ we get $y_{n+1} =R(z)y_{n},$ where R(z) is the associated stability function with z=hλ. Moreover, from , , , we obtain recursively $k_{1} =λy_{n},$ $s_{2} =zG_{2},$ $s_{3} =zG_{3} (zG_{2}),$ $⋮$ $s_{p} =zG_{p} (zG_{2},zG_{3} (zG_{2}),…,zG_{p−1} (zG_{2},…,zG_{p−2} (⋯))),$ from which we obtain the following expression for the stability function: $R(z)=1+zG_{p+1} (s_{2},s_{3},…,s_{p}),$ where the s_i 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) $y^{′} (x)= y(x)(1−y(x)) 2y(x)−1, y(0)= 56,$ for which the solution is $y(x)= 12 + 14 − 536 e^{−x} .$ With fixed step size h=2⁻ⁿ for various n=1,2,3,…,9 over 2ⁿ 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)

H. Zedan
Avoiding the exactness of the Jacobian matrix in Rosenbrock formulae
Comput. Math. Appl.
(1990)
J.D. Lambert
Numerical Methods for Ordinary Differential Systems. The Initial Value Problem
(1991)
G. Dahlquist
A special stability problem for linear multistep methods
BIT
(1963)
J.C. Butcher
The Numerical Analysis of Ordinary Differential Equations: Runge–Kutta and General Linear Methods
(1987)
J.C. Butcher
Implicit Runge–Kutta processes
Math. Comp.
(1964)
E. Hairer et al.
Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems
(1996)
R. Alexander
Diagonally implicit Runge–Kutta methods for stiff O.D.E.'s
SIAM J. Numer. Anal.
(1977)
H.H. Rosenbrock
Some general implicit processes for the numerical solution of differential equations
Comput. J.
(1963)
S.P. Nørsett et al.
Order conditions for Rosenbrock type methods
Numer. Math.
(1979)
P. Kaps et al.
Generalized Runge–Kutta methods of order four with stepsize control for stiff ordinary differential equations
Numer. Math.
(1979)

P. Kaps et al.

A study of Rosenbrock-type methods of high order

Numer. Math.

(1981)

Cited by (19)

Time discretization of the point kinetic equations using matrix exponential method and First-Order Hold
2013, Annals of Nuclear Energy
Citation 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.
A system of nonlinear differential equations is derived to model the dynamics of neutron density and the delayed neutron precursors within a point kinetics equation modeling framework for a nuclear reactor. The point kinetic equations are mathematically characterized as stiff, occasionally nonlinear, ordinary differential equations, posing significant challenges when numerical solutions are sought and traditionally resulting in the need for smaller time step intervals within various computational schemes. In light of the above realization, the present paper proposes a new discretization method inspired by system-theoretic notions and technically based on a combination of the matrix exponential method (MEM) and the First-Order Hold (FOH) assumption. Under the proposed time discretization structure, the sampled-data representation of the nonlinear point kinetic system of equations is derived. The performance of the proposed time discretization procedure is evaluated using several case studies with sinusoidal reactivity profiles and multiple input examples (reactivity and neutron source function). It is shown, that by applying the proposed method under a First-Order Hold for the neutron density and the precursor concentrations at each time step interval, the stiffness problem associated with the point kinetic equations can be adequately addressed and resolved. Finally, as evidenced by the aforementioned detailed simulation studies, the proposed method retains its validity and accuracy for a wide range of reactor operating conditions, including large sampling periods dictated by physical and/or technical limitations associated with the current state of sensor and digital reactor control system technology.
Convergence on error correction methods for solving initial value problems
2012, Journal of Computational and Applied Mathematics
Higher-order semi-explicit one-step error correction methods(ECM) for solving initial value problems are developed. ECM provides the excellent convergence $O (h^{2 p + 2})$ one wants to get without any iteration processes required by most implicit type methods. This is possible if one constructs a local approximation having a residual error $O (h^{p})$ on each time step. As a practical example, we construct a local quadratic approximation. Further, it is shown that special choices of parameters for the local quadratic polynomial lead to the known explicit second-order methods which can be improved into a semi-explicit type ECM of the order of accuracy $6$ . The stability function is also derived and numerical evidences are presented to support theoretical results with several stiff and non-stiff problems. It should be remarked that the ECM approach developed here does not yield explicit methods, but semi-implicit methods of the Rosenbrock type. Both ECM and Rosenbrock’s methods require to solve a few linear systems at each integration step, but the ECM approach involves $2 p + 2$ evaluations of the Jacobian matrix per integration step whereas the Rosenbrock method demands one evaluation only. However, it is much easier to get high order methods by using the ECM approach.
Generalized Runge-Kutta method for two- and three-dimensional space-time diffusion equations with a variable time step
2008, Annals of Nuclear Energy
An extensive knowledge of the spatial power distribution is required for the design and analysis of different types of current-generation reactors, and that requires the development of more sophisticated theoretical methods. Therefore, the need to develop new methods for multidimensional transient reactor analysis still exists. The objective of this paper is to develop a computationally efficient numerical method for solving the multigroup, multidimensional, static and transient neutron diffusion kinetics equations. A generalized Runge–Kutta method has been developed for the numerical integration of the stiff space–time diffusion equations. The method is fourth-order accurate, using an embedded third-order solution to arrive at an estimate of the truncation error for automatic time step control. In addition, the A(α)-stability properties of the method are investigated. The analyses of two- and three-dimensional benchmark problems as well as static and transient problems, demonstrate that very accurate solutions can be obtained with assembly-sized spatial meshes. Preliminary numerical evaluations using two- and three-dimensional finite difference codes showed that the presented generalized Runge–Kutta method is highly accurate and efficient when compared with other optimized iterative numerical and conventional finite difference methods.
An improved class of generalized Runge-Kutta methods for stiff problems. Part II: The separated system case
2004, Applied Mathematics and Computation
Citation 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.
In order to study a generalization of the explicit p-stage Runge–Kutta methods, we introduce a new family of p-stage formulas for the numerical integration of some special systems of ODEs that provides better order and stability results with the same number of stages. In our recent paper entitled `An improved class of generalized Runge–Kutta methods for stiff problems. Part I: The scalar case' we studied new schemes for the numerical integration of scalar autonomous ODEs. In this second part we will show that it is possible to generalize our GRK-methods so that they can be applied to some non-autonomous scalar ODEs and systems obtaining linearly implicit A-stable and L-stable methods. These methods do not require Jacobian evaluations in their implementation. Some numerical examples are discussed in order to show the good performance of the new schemes.
An improved class of generalized Runge-Kutta-Nyström methods for special second-order differential equations
2004, Communications in Nonlinear Science and Numerical Simulation
We propose a new family of explicit methods of order four with two evaluations per step, for the numerical integration of special second-order differential equations given by y^″=f(y). These two-stage formulas can be seen as a generalization of the explicit two-stage Runge–Kutta–Nyström methods, providing better order and stability results. We will show that it is possible to obtain methods that are more efficient than the classical Runge–Kutta–Nyström one-step methods with the same number of evaluations per step, specially when highly oscillatory problems are considered. Some numerical experiments are discussed in order to show the good performance of the new schemes.
On pade' approximations to the exponential function and application to the point kinetics equations
2004, Progress in Nuclear Energy
A general expression for a special type of functions has been introduced, which allows us to approximate the exponential function in an economical manner. Also, the different cases of Padé approximation are perturbed so that the resulting approximations have a smaller minimum maximum error on the desired interval especially at large transient time steps. In addition, the Padé approximations and the resulting perturbed ones are applied to the solution of the point kinetics equations of nuclear reactor with different types of reactivity. The results of the approximations developed show large correction effects at small transients and perform quite well at large transients. Also, the factorization methods based on using temporarily the complex plane with the analytical inversion are applied to obtain the general solution of the point kinetics equations.

View all citing articles on Scopus

¹: 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.

²: 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.

View full text

An improved class of generalized Runge–Kutta methods for stiff problems. Part I: The scalar case

Abstract

Introduction

Section snippets

The new family of methods

A useful formulation

Consistency and order of the methods

Methods of polynomial type

Two-stage methods of polynomial type: attainable order

Three-stage methods of polynomial type: attainable order

Methods of rational type

Two-stage methods of rational type

Three-stage methods of rational type

Linear stability properties of the methods

Three-stage methods being A-stable

Numerical experiments (the scalar case)

Conclusions

Comput. Math. Appl.

Numerical Methods for Ordinary Differential Systems. The Initial Value Problem

A special stability problem for linear multistep methods

BIT

The Numerical Analysis of Ordinary Differential Equations: Runge–Kutta and General Linear Methods

Implicit Runge–Kutta processes

Math. Comp.

Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems

Diagonally implicit Runge–Kutta methods for stiff O.D.E.'s

SIAM J. Numer. Anal.

Some general implicit processes for the numerical solution of differential equations

Comput. J.

Order conditions for Rosenbrock type methods

Numer. Math.

Generalized Runge–Kutta methods of order four with stepsize control for stiff ordinary differential equations

Numer. Math.

A study of Rosenbrock-type methods of high order

Numer. Math.