ARTICLES

A new numerical method for stiff differential equations

G. Boroni^†, P. Lotito^† and A. Clausse^†‡

^† CONICET and Universidad Nacional del Centro, Tandil, Argentina.
^‡ Comisión Nacional de Energía Atómica, Argentina.
{gboroni, plotito, clausse}@exa.unicen.edu.ar

Abstract -A new class of multistep methods for stiff ordinary differential equations is presented. The method is based in the application of estimation functions for the derivatives and the state variables, allowing the transformation of the original system in a purely algebraic system using the solutions of previous steps. From this point of view these methods adopt a semi-implicit scheme. The novelty introduced is an adaptive formula for the estimation function coefficients, deduced from a combined analysis of stability and convergence order. That is, the estimation function coefficients are recomputed at each time step. The convergence order of the resulting scheme is better than the equivalent linear multistep methods, while preserving the properties of stability.

Keywords -Multistep Methods. A-Stability. Convergence Order.

I. INTRODUCTION

Solving differential equations numerically is, even today, still a great challenge. This applies especially to stiff differential equations and to closely related problems involving differential equations. Linear multistep methods and, in particular, Backward Differentiation Formulae (BDF) (Gear, 1971; and Ascher and Petzold, 1998) are regularly used for the numerical solution of stiff initial value problems.

A high-quality numerical method for the solution of stiff Ordinary Differential Equations (ODE) must have good accuracy and some wide region of absolute stability (Dahlquist, 1963; Enright, 1974; Ismail and Ibrahim, 1999; and Lambert, 1972). The latter imposes a strong limitation on the choice of suitable methods for stiff problems.

In general, the search for high-order A-stable methods has been carried out in the two main directions:

• Use of higher derivatives of the solutions

• Introduction of additional stages, off-step points, and super-future points (this leads into the large field of general linear methods; Hairer and Wanner 1996).

In the present paper a new class of multistep methods is derived, having good stability properties and improvements in the convergence order compared with equivalent linear schemes. The difference with the others multistep methods is the application of estimation functions not only of the derivatives but also of the state variables, which leads to the transformation of the original system into a purely algebraic system. The method proposed is an improvement over the linear multistep schemes in the sense that it incorporates nonlinear correction coefficients while keeping A-stability and high-order of convergence. In the last section, numerical experiments are presented comparing the new method with BDF.

II. LINEAR MULTISTEP METHOD

Let us consider the initial value problem

(1)

where t ∈ [t₀, t₀ + Nh] (N being a natural number and h a constant time step), y : [t₀, t₀ + Nh] → R^m, y⁽¹⁾ stands for the first temporal derivative, and f : R^m → R^m is continuous and differentiable.

The general multistep method can be written in the form (Ascher and Petzold, 1998)

(2)

where α_i,β_i are parameters to be determined, y_n = y(t₀ + nh). A multistep method is of order p if and only if (Butcher 2003)

(3)

with 0 ≤ q ≤ p.

The well-known multistep scheme BDF is given by

(4)

This scheme is a class of k-step formula of order k, and for k = 2, the BDF coefficients are

(5)

III. NEW MULTISTEP METHOD

The general multistep formula (i.e. Eq. 2) is essentially a transformation of the differential Eq. 1 into a purely algebraic equation by means of the estimators

(6)

In this work, let us propose the following set of transformations for , 1 ≤ j ≤ m, and

(7)

where A_i,j and B_i,j are coefficients satisfying and .

Equations 7 lead to the following alternative multistep algebraic equation

(8)

To make clear the deduction of the coefficients, let us consider a two-step instance of Eqs. 8 (k = l = 2). The general method to determine the coefficients A_i,j and B_i,j, can be easily deduced from this particular case.

A. Convergence order

Expanding y_n-2,j and y_n,j in Taylor series around (t-h) leads to

(9)

where stands for the k-th derivative of y respect to time. Combining Eqs. 7 and 9 yields

(10)

Expanding f_j around (t-h), gives

(11)

where and H_j(y_n) stand for the transpose gradient and the Hessian of f_j(y) evaluated at y = y_n, respectively.

Likewise, expanding the right term of the Eq. 8, gives

(12)

The relation between f_j(y_n), and H_j(y_n) with can be found by successive differentiations of Eq. 1. Combining Eqs. 10 to 12 yields

(13)

which leads to the following set of algebraic equations for the coefficients A_i,j and B_i,j

(14)

where J(y_n) stands for the Jacobian of f(y) evaluated at y = y_n. Eqs. 14 is a set of 5×m algebraic equations with 6×m unknowns. Therefore, there is a family of coefficients A_i,j and B_i,j that ensures O(h³) convergence

(15)

In Eqs. 15, f, H_j and are evaluated at y_n-1. Hence, every term c_j can be seen as a non-linear correction applied each step to the estimator coefficients of y_j in order to increase the convergence order of the scheme. From the Eq. 15 show that all the coefficients B_i,j are the same for all derivative, while the coefficients A_i,j are different for each variable through c_j. Moreover, if f is linear (i.e. H_j = 0), c_j = 0, and all the coefficients A_i,j are the same for all variables y_j.

B. A-stability

In order to choose one particular set of coefficients from the O(h³) convergence families, A-stability is required. A method is A-stable if, applied to a stable linear set of differential equations, the resulting iterative scheme is also stable independently of h. In that way, ensuring A-stability, h is determined just for precision purposes, without restrictions on linear stability. Such methods are considered good candidates to solve stiff problems (Ascher and Petzold, 1998).

If f is linear, all the coefficients of one variable are independent of those of the other variables, and A-stability can be analyzed with a single variable. Applying Eqs. 2 to the linear test equation y⁽¹⁾ = λy yields

(16)

The stability of Eq. 16 is ensured if the real part of the roots of the characteristic polynomial

(17)

are negative. A-stability is then given by (Ascher and Petzold, 1998) Re(z) ≥ 0, where

(18)

for all (unitary) complex numbers q = cosθ + isinθ, for θ ∈ [0,2π]. Next, the scheme k = l = 2 is A-stable if (Fig. 1)

(19)

Fig. 1. Absolute stability region of Hybrid Multistep Method (HMM) and BDF methods.

IV. NUMERICAL CONSIDERATIONS

The term c_j in the Eq. 15 has a singularity if . In fact, whenever the magnitude becomes very small, the coefficient c_j becomes very large, and the numerical solution of the implicit algebraic Eq. 8 fails. In those singular time steps, an ad-hoc solution can be chosen, either temporary compromising the stability or temporary reducing the convergence order. A good practical alternative is to switch to the BDF method -which is A-stable and O(h³) (Ascher and Petzold, 1998)- whenever any of the absolute values of the non-linear terms, |c_j|, exceeds a critical threshold.

V. NUMERICAL EXPERIMENTS

In order to assess the performance of the new multistep method, it was applied to the integration of specific equations, and compared to the BDF method. The studies were carried out considering B₁=0.001.

A. Ricatti equation (m = 1)

Let us consider the following Ricatti equation (Abramowitz and Stegun, 1972)

(20)

with initial value y₀ = 1.8. The exact solution is given by

(21)

Figure 2 shows the temporal evolution of f⁽⁰⁾, f⁽¹⁾ and f⁽²⁾, the coefficients A_i and the non-linear correction term c. BDF is activated when |c| (between 0.4 and 1.4). Figure 3 shows the absolute difference between the analytic and the numerical solutions. It can observe that the numerical solution obtained with the new method is always better than the BDF solution.

Fig. 2. Temporary evolution of f⁽⁰⁾, f⁽¹⁾ and f⁽²⁾, the A-coefficients and the non-linear correction term c.

Fig. 3. Calculation of the absolute difference |e| between the analytic and the numerical solutions, with h = 10^-3.

B. Van der Pol oscillator

The Van der Pol (VDP) oscillator is a second order system that can be derived from the Rayleigh equation (Thompson, 1986). The system is dissipative and leads to limit cycles. A simple form of VDP in terms of first order ODE is

(22)

where μ is a constant parameter which determines the size of the limit cycle. The corresponding Jacobian and Hessians of f₁ and f₂ are

(23)

The non-linear correction terms are given by

(24)

Numerical experiments were performed for the case μ = 3.5. Figure 4 shows the temporal evolution of the state variables, and the A-coefficients of y₂. It can be seen that the coefficients vary significantly, not only at the "elbows" of the cycle but also during the apparently smooth periods of the phase-space trajectory.

Fig. 4. Temporary evolution of the state variables, and the coefficients A_0,2, A_1,2, A_2,2.

VI. CONCLUSIONS

A 2-step method O(h³) by means of non-linear corrections of the estimation-function coefficients was presented. This scheme is A-stable, and it is a good candi- date method for the solution of stiff problems. The new method requires an update of the coefficients at each time step, and, as a consequence, an increase in the computational cost, however the new method benefits of better precision and a larger stability region, as it is shown in the proposed examples. As future work of the new method was compared with other methods for stiff problems, as for example the so-called exponential integrator method (De la Cruz et al., 2007).

REFERENCIAS

1. Abramowitz, M. and I.A. Stegun, "Riccati-Bessel Functions." Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th edition. New York, Dover (1972).         [ Links ]

2. Ascher, U. and L. Petzold, Computer Method for Ordinary Differential Equations and Differential Algebraic Equations, Siam (1998).         [ Links ]

3. Butcher, J., Numerical Method for Ordinary Differential Equations, Wiley (2003).         [ Links ]

4. Dahlquist, G., "A special stability problem for linear multistep methods," BIT Numerical Mathematics, 3, 27-43 (1963).         [ Links ]

5. De la Cruz, H., R.J. Biscay, F. Carbonell, T. Ozaki and J.C. Jimenez, "a higher order local linearization method for solving ordinary differential equations," Applied Mathematics and Computation, 185, 197-212 (2007).         [ Links ]

6. Enright, W.H., "Second derivative multistep methods for stiff ordinary differential equations," SIAM J. Numer. Anal., 11, 321-331 (1974).         [ Links ]

7. Gear, W., "Simultaneous Numerical Solution of Differential Algebraic Equation," IEEE Transaction on Circuit Theory, 18, 89-95 (1971).         [ Links ]

8. Hairer, E. and G. Wanner, Solving Ordinary Differential Equation II: Stiff and Differential-Algebraic Problems, Springer, Berlin (1996).         [ Links ]

9. Ismail, G. and I. Ibrahim, "New efficient second derivative multistep methods for stiff systems," Appl. Math. Model, 23, 279-288 (1999).         [ Links ]

10. Lambert, J.D., Computational Methods in Ordinary Differential Equations, John Wiley and Sons (1972).         [ Links ]

11. Thompson, J. and H. Stewart, Nonlinear Dynamics and Chaos. Wiley (1987).         [ Links ]

Received: December 20, 2007.
Accepted: March 6, 2008.
Recommended by Subject Editor José Guivant.