Exponential fitting BDF algorithms and their properties

doi:10.1016/j.amc.2007.01.008

Applied Mathematics and Computation

Volume 190, Issue 1, 1 July 2007, Pages 80-110

https://doi.org/10.1016/j.amc.2007.01.008 Get rights and content

Abstract

We present two families of explicit and implicit BDF formulas, capable of the exact integration (with only round-off errors) of differential equations whose solutions belong to the space generated by the linear combinations of exponential of matrices, products of the exponentials by polynomials and products of those matrices by ordinary polynomials. Those methods are suitable for stiff and highly oscillatory problems, then we will study their properties of local truncation error and stability. Plots of their 0-stability regions in terms of Ah are provided. Plots of their regions of absolute stability that include all the negative real axis in the implicit case are provided. Exponential fitting algorithms depend on a parameter Ah, how can we find this parameter is a big question, here we show different ways to find a good parameter. Finally, numerical examples show the efficiency of the proposed codes, specially when we are integrating stiff problems where the jacobian of the function has complex eigenvalues or problems where the jacobian has positive eigenvalues but the solutions of the problems have not positive exponentials.

Introduction

In this paper, we are going to derive explicit and implicit BDF methods, that solve the IVP problem $y^{'} (x) = G (x, y (x)), y (x_{0}) = y_{0},$ where $y = [y^{1}, \dots, y^{m}]^{t}$ , and $G = [g^{1}, \dots, g^{m}]^{t}$ .

The most important improve of these new methods shall be that they are able to integrate problems where the jacobian has positive eigenvalues but the solutions of the problems have not positive exponentials. We will see that these problems are very hard, and most of the classical methods (including classical methods for stiff problems) are not stable with these problems. We are going to show very good results in problems that appear in a recent bibliography (see [11] or [25]).

The new methods are very efficient with other two kind of numerical examples: stiff and highly oscillatory problems. The most common definition of stiffness (see Definition 2.3 in [39], and for example, [11] or [25]) is: stiffness occurs if for most classical explicit methods, the largest step size $h_{n}^{*}$ , guaranteeing numerical stability is much smaller than the largest step h_n for which the local discretization error is still sufficiently small (in norm), i.e., $h_{n}^{*} ≪ h_{n}$ . With this definition stiff problems include highly oscillatory problems, however, many scientists do not consider highly oscillatory problems as stiff problems.

Stiff problems (and highly oscillatory problems) are very common problems in many fields of the applied sciences (see [1]): atmosphere, biology, combustion, control theory, dynamics of missile guidance, dispersed phases, electronic circuit theory, fluids, heat transfer, chemical kinetics (this is the most common area), lasers, mechanics, molecular dynamics, nuclear, process industries, process vessels, reactor kinetics, … Then, stiff computation is a very interesting area and it is impossible to include all the contributions.

Many methods have been considered to integrate stiff problems, but implicit Runge–Kutta (Radau, see [17]; STRIDE, see [2], [10] or [32]; DIRK, [2] or [7]; SDIRK, Gauss, Rosembrock, …) and implicit BDF (MEBDF, see [5]; DASSL, see [34] or [35]; LSODE, see [18], …) codes have frequently been used.

The idea of using exponential fitting/adapted formulae with the numerical integration of stiff systems has recently received considerable attention (see [6], [22], [26], …).

The basic idea behind exponential fitting is to integrate exactly (1) when y(x) belongs to a large space where linear combinations of exponential of matrices are included, and we will denote exponential fitting algorithms, in this paper, with the capital letters EF.

While we will use the name adapted methods with the methods that solve the IVP problem $y^{'} (x) - Ay (x) = F (x, y (x)), y (x_{0}) = y_{0}$ (where $y = [y^{1}, \dots, y^{m}]^{t}$ , A is a m × m constant matrix and $F = [f^{1}, \dots, f^{m}]^{t}$ ) exactly when y(x) belongs to a space where linear combinations of exponential of matrices are included, here, we will denote those algorithms with the capital letter A.

This is, exponentially fitting methods will depend on G(x, y(x)), while adapted methods (see [27], [42]) depend on $F (x, y (x)) = G (x, y (x)) - Ay (x)$ . In this case, implicit adapted methods function as explicit methods with linear problems.

Exponential fitting and adapted methods are very similar and we will show a relation between them in the following sections. In Section 2 we will derive two families of exponential fitting methods, while in Section 3 we will provide the local truncation error of those methods. In Section 4 we will provide plots of the 0-stability of the methods, and then, in Section 5 the definition of absolute stability of the methods and some plots of those regions are provided, in those Sections 4 0-Stability, 5 Absolute stability regions, we will see that adapted and exponential fitting algorithms have different properties of stability and we will explained this apparent contradiction. Finally in Section 6 we will compare our algorithms with well known ODE solvers for the numerical solution of many different kind of problems.

Section snippets

Derivation of the exponential fitting/adapted methods

Through the basic theory of difference equations (see [41]), we are going to build two families of exponential fitting/adapted methods. The first one integrate with no truncation error the problem (1) when y(x) is a vector of functions that belongs to the space generated by the linear combinations of $〈 e^{Ax}, 1, x, \dots, x^{r - 1} 〉$ (we will denote them with EF-I-r or A-I-r), while the second family integrate exactly the problem (1) when y(x) is a vector of functions that belongs to the space generated by the

Local truncation error

Theorem 3

The multistep implicit adapted method of type I $\nabla_{P} \sum_{j = 0}^{r - 1} β_{j}^{0} \nabla^{j} y_{n} = {hf}_{n}$ of r steps is convergent of order r. Its local truncation error can be expressed in the form: $h^{r + 1} β_{r}^{0} (D - A) D^{r} y (x) + O (h^{r + 2}) .$

Proof

As we have seen in [27] the operator providing us with the local truncation error of (21) is written in the form $\nabla_{P} β_{r}^{0} \nabla^{r} y (x) + O (h^{r + 2}),$ and the differences satisfy $\nabla^{r} y (x) = h^{r} D^{r} y (x) + O (h^{r + 1})$ , and $\nabla_{P} y (x) = h (D - A) y (x) + O (h^{2})$ , then the local truncation error of (21) is written in the form $h^{r + 1} β_{r}^{0} (D - A) D^{r} y (x) + O (h^{r + 2}) . □$

Theorem 4

The

0-Stability

It is common to assume that if the method is exponential fitting or adapted then the method is 0-stable, but this fact is not true in the case of BDF type formulas. There are some regions where the method is not 0-stable although it is really easy to select a step-size that avoid those regions.

In attempting to set up zero-stability, we consider the linear constant homogeneous system $y^{'} (x) = Ay (x),$ where the eigenvalues $λ_{i}, i = 1, \dots, m$ of the m × m matrix A satisfy $Re (λ_{i}) < 0$ .

It follows that solutions of the

Absolute stability regions

The classical definitions of regions of absolute stability and A-stability were designed for linear multistep methods with constant coefficients. In this section those definitions are modified so as to provide a basis for linear stability analysis of exponential fitting methods for ordinary differential equations of type (1). The stability properties of proposed methods are analyzed, to demonstrate its relevance to stiff oscillatory problems. The way is very similar to that used in [8] to

Numerical examples

We have explained that the improve of these methods is not only the fact that one is explicit, but the behavior with problems where the jacobian of the function has a positive eigenvalue (but the solution of the ODE has not positive exponentials) or with stiff oscillatory problems, too.

We will say that one ODE is stiff oscillatory if $y^{'} = f (x, y)$ and the jacobian f_y has eigenvalues $λ_{i}, λ_{j}$ , where $Re (λ_{i}) ≪ 0$ comparing with the interval of integration and $| \frac{Im (λ_{j})}{Re (λ_{j})} | ≫ 1$ .

In this kind of problems we

Conclusions

We have shown different types of exponential fitting methods that integrates (1) when y(x) (the solution of the problem) belongs to the space generated by the linear combinations of $〈 e^{Ax}, \dots, x^{k - 1} e^{Ax}, 1 〉$ or the space generated by $〈 e^{Ax}, 1, x, \dots, x^{r - 1} 〉$ .

In [28] the A-I-r algorithms are shown. The A-I-r methods and the EF-I-r algorithms really verify the same methods, the only one difference is the presentation of the ODE solver.

We have shown that the absolute stability regions of the EF-I-r and the A-I-r

References (43)

W.S. Edwards et al.
Krylov methods for the incompressible Navier–Stokes equations
J. Comput. Phys.
(1994)
N. Guzel et al.
On the numerical solution of stiff systems
Appl. Math. Comput.
(2005)
L.Gr. Ixaru et al.
Frequency evaluation in exponential fitting multistep algorithms for ODEs
J. Comput. Appl. Math.
(2002)
J. Martín-Vaquero et al.
Exponential fitting BDF algorithms: explicit and implicit 0-stable methods
J. Comp. Appl. Math.
(2006)
M.N. Spijker
Stiffness in the numerical initial-value problems
J. Comput. Appl. Math.
(1996)
J. Vigo-Aguiar et al.
On the stability of exponential fitting BDF algorithms
J. Comput. Appl. Math.
(2005)
R.C. Aiken
Stiff Computation
(1985)
R. Alexander
Diagonally implicit Runge–Kutta methods for stiff ODEs
SIAM J. Numer. Anal.
(1977)
D. Barton
On Taylor series and stiff equations
ACM Trans. Math. Software
(1980)
N.N. Bogoliubov et al.
Asymptotic Methods in the Theory of Non-linear Oscillations
(1985)

J.R. Cash et al.

An MEBDF code for stiff initial value problems

ACM Trans. Math. Soft.

(1992)

J.R. Cash

On the exponential fitting of composite multiderivative linear multistep methods

SIAM J. Numer. Anal.

(1981)

J.R. Cash

Diagonally implicit Runge–Kutta formulae with error estimates

J. IMA

(1979)

J.P. Coleman et al.

P-stability and exponential-fitting methods for y″ = f(x, y)

IMA J. Numer. Anal.

(1996)

G.F. Corliss et al.

Hihg-order ODE solvers via automatic differentiation and rational prediction

(1997)

M. Crouzeix, Sur l’approximation des equations differentiales operationelles lineaires par des methodes Runge–Kutta,...

K. Dekker et al.

Stability of Runge–Kutta methods for stiff nonlinear differential equations

(1984)

W.H. Enright et al.

Comparing numerical methods for stiff systems of ODE’s

BIT

(1975)

R.A. Friesner et al.

A method for exponential propagation of large systems of stiff nonlinear differential equations

J. Sci. Comput.

(1989)

E. Gallopoulos et al.

Efficient solution of parabolic equations by Krylov approximation methods

SIAM J. Sci. Stat. Comput.

(1992)

E. Hairer et al.

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

(1991)

Cited by (14)

Present state-of-the-art in exponential fitting. A contribution dedicated to Liviu Ixaru on his 70th birthday
2012, Computer Physics Communications
The standard monograph in this area is the book Exponential fitting by Ixaru and Vanden Berghe (Kluwer, Boston - Dordrecht - London, 2004) but a fresh look on things is necessary because many new contributions have been accumulated in the meantime. With no claim that our investigation is exhaustive we consider various directions of interest, try to integrate the new contributions in a natural, easy to follow way, and also detect some open problems of acute interest.
On the numerical solution of the heat conduction equations subject to nonlocal conditions
2009, Applied Numerical Mathematics
Many physical phenomena are modelled by nonclassical parabolic boundary value problems with nonlocal boundary conditions. Many different papers studied the second-order parabolic equation, particularly the heat equation subject to the specifications of mass. In this paper, we provide a whole family of new algorithms that improve the CPU time and accuracy of Crandall's formula shown in [J. Martin-Vaquero, J. Vigo-Aguiar, A note on efficient techniques for the second-order parabolic equation subject to non-local conditions, Appl. Numer. Math. 59 (6) (2009) 1258–1264] (and this algorithm improved the results obtained with BTCS, FTCS or Dufort–Frankel three-level techniques previously used in other works, see [M. Dehghan, Efficient techniques for the second-order parabolic equation subject to nonlocal specifications, Appl. Numer. Math. 52 (2005) 39–62]) with this kind of problems. Other methods got second or fourth order only when $k = s h^{2}$ , while the new codes got nth order for $k = h$ ; therefore, the new schemes require a smaller storage and CPU time employed than other algorithms. We will study the convergence of the new algorithms and finally we will compare the efficiency of the new methods with some well-known numerical examples.
Exponential fitting BDF-Runge-Kutta algorithms
2008, Computer Physics Communications
Citation Excerpt :
Several papers (see [40], for example) showed, using perturbation theory, a way to get linear equations from nonlinear ones. The solutions of both IVP's are very similar, then we can choose the parameter Ah from the linear problem and Ah allows us to get good results as we did in [28] with several stiff oscillatory problems. Now, in the numerical example 3, we shall show that with the new EF-BDF-RK methods we get better results than with the EF-BDF algorithms that we showed in [28].
In other papers, the authors presented exponential fitting methods of BDF type. Now, these methods are used to derive some BDF–Runge–Kutta type formulas (of second-, third- and fourth-order), capable of the exact integration (with only round-off errors) of differential equations whose solutions are linear combinations of an exponential with parameter A and ordinary polynomials. Theorems of the truncation error reveal the good behavior of the new methods for stiff problems. Plots of their absolute stability regions that include the whole of the negative real axis are provided. Different procedures to find the parameter of the method are proposed, using these techniques there will not be necessary to compute the exponential matrix at each step, even when nonlinear problems are integrated. Numerical examples underscore the efficiency of the proposed codes, especially when they are integrating stiff problems.
Fractional exponential fitting backward differential formulas for solving differential equations of fractional order
2023, Computational and Applied Mathematics
An explicit trigonometrically fitted Runge–Kutta method for stiff and oscillatory problems with two frequencies
2020, International Journal of Computer Mathematics
An adapted explicit hybrid four-step method for the numerical solution of perturbed oscillators
2018, Journal of Mathematical Chemistry

View all citing articles on Scopus

View full text

Exponential fitting BDF algorithms and their properties

Abstract

Introduction

Section snippets

Derivation of the exponential fitting/adapted methods

Local truncation error

0-Stability

Absolute stability regions

Numerical examples

Conclusions

J. Comput. Phys.

Appl. Math. Comput.

J. Comput. Appl. Math.

J. Comp. Appl. Math.

J. Comput. Appl. Math.

J. Comput. Appl. Math.

Stiff Computation

Diagonally implicit Runge–Kutta methods for stiff ODEs

SIAM J. Numer. Anal.

On Taylor series and stiff equations

ACM Trans. Math. Software

Asymptotic Methods in the Theory of Non-linear Oscillations

An MEBDF code for stiff initial value problems

ACM Trans. Math. Soft.

On the exponential fitting of composite multiderivative linear multistep methods

SIAM J. Numer. Anal.

Diagonally implicit Runge–Kutta formulae with error estimates

J. IMA

P-stability and exponential-fitting methods for y″ = f(x, y)

IMA J. Numer. Anal.

Hihg-order ODE solvers via automatic differentiation and rational prediction

Stability of Runge–Kutta methods for stiff nonlinear differential equations

Comparing numerical methods for stiff systems of ODE’s

BIT

A method for exponential propagation of large systems of stiff nonlinear differential equations

J. Sci. Comput.

Efficient solution of parabolic equations by Krylov approximation methods

SIAM J. Sci. Stat. Comput.

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