Mathematical Problems in Engineering

Volume 2013 (2013), Article ID 989381, 11 pages

http://dx.doi.org/10.1155/2013/989381

Research Article

A Numerical Algorithm for Solving Stiff Ordinary Differential Equations

S. A. M. Yatim,¹ Z. B. Ibrahim,² K. I. Othman,³ and M. B. Suleiman²

¹School of Distance Education, Universiti Sains Malaysia, 11800 Penang, Malaysia

²Department of Mathematics, Faculty of Science, Universiti Putra Malaysia, Selangor, 43400 Serdang, Malaysia

³Department of Mathematics, Faculty of Computer and Mathematical Sciences, Universiti Teknologi MARA, Selangor, 40450 Shah Alam, Malaysia

Received 6 August 2012; Revised 24 October 2012; Accepted 7 November 2012

Academic Editor: J. Rodellar

Copyright © 2013 S. A. M. Yatim et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

An advanced method using block backward differentiation formula (BBDF) is introduced with efficient strategy in choosing the step size and order of the method. Variable step and variable order block backward differentiation formula (VSVO-BBDF) approach is applied throughout the numerical computation. The stability regions of the VSVO-BBDF method are investigated and presented in distinct graphs. The improved performances in terms of accuracy and computation time are presented in the numerical results with different sets of test problems. Comparisons are made between the proposed method and MATLAB’s suite of ordinary differential equations (ODEs) solvers, namely, ode15s and ode23s.

1. Introduction

Many studies on solving the equations of stiff ordinary differential equations (ODEs) have been done by researchers or mathematicians specifically. With the numbers of numerical methods that currently exist in the literature, extensive research has been done to unveil the comparison between their rate of convergence, number of computations, accuracy, and capability to solve certain type of test problems [1]. The well-known numerical methods that are used widely are from the class of BDFs or commonly understood as Gear’s Method [2]. However, many other methods that have evolved to this date are for solving stiff ODEs which arise in many fields of the applied sciences [3–6]. The problems considered in this paper are for the numerical solution of the initial value problem,

with given initial values in the given interval .

The existing numerical methods have often been compared to one another to find the best approximation and the best method. MATLAB ODE suite is one of the most preferred solvers to be used for comparison purposes [7, 8]. Stiff ODE solvers that are available in MATLAB ODE suite are ode15s and ode23s which are based on the numerical differentiation formulas [8, 9] and the modified Rosenbrock formula of order 2, respectively [8].

An improved method has been identified in [10], which was expanded from the method incorporated with BDF proposed by Gear. Since then, the study on producing block approximations also known as Block Backward Differentiation Formulae (BBDF) [11] has attracted much attention. BBDF method by using variable step size in [12] demonstrates the competency of computing concurrent solution values at different points. Consequently, study in [13, 14] is an extension of a previous study in a way that the accuracy is improved by increasing the order of the method up to order 5. Thus, these studies lead to enhancing the existing method to become Variable Order Variable Step Block Backward Differentiation Formula of order 3 until 5 (VSVO-BBDF).

In Section 2, we introduce the general formulation of VSVO-BBDF of order 3 to 5. Order conditions are listed in Section 3 while in Section 4, we consider the implementation of VSVO-BBDF. The analysis of stability regions is illustrated in Section 5. This is followed by the strategy in choosing the step size and order of VSVO-BBDF in Section 6. Numerical results are presented in Section 7. The appendix describes the algorithm applied for VSVO-BBDF method.

2. VSVO-BBDF Method Formulation

Two values of and were computed simultaneously in block by using earlier blocks with each block containing a maximum of two points (Figure 1). The orders of the method (,  , and  ) are distinguished by the number of back values contained in total blocks. The ratio distance between current and previous step is represented as and in Figure 1. In this paper, the step size is given selection to decrease to half of the previous steps or increase up to a factor of 1.9. For simplicity, is assigned as 1, 2, and 10/19 for the case of constant, halving and increasing the step size, respectively. The zero stability is achieved for each of these cases and explained in the next section.

Figure 1: VSVO-BBDF method of order (–).

We find approximating polynomials , by means of a  -degree polynomial interpolating the values of at given points, that are :

where

Predictors for the first point and second point were computed using back values as the interpolating points. The resulting Lagrange polynomial for each order was given as follows.

For VSVO-BBDF of order

For VSVO-BBDF of order

For VSVO-BBDF of order

Substituting and gives the predictor for the first and second points, respectively. Therefore by letting , and . This produced the following coefficients (Tables 1, 2, and 3) for the first and second points of predictor formulae for VSVO-BBDF method.

Table 1: Coefficients for the first and second point of predictor formulae for VSVO-BBDF when , and .

Order

Table 2: Coefficients for the first and second point of predictor formulae for VSVO-BBDF when , and .

Order

Table 3: Coefficients for the first and second point of predictor formulae for VSVO-BBDF when , and .

Order

The interpolating polynomial of the function using Lagrange polynomial in (2) gives the following corrector for the first point and second point . The resulting Lagrange polynomial for each order was given as follows.

For VSVO-BBDF of order

For VSVO-BBDF of order

For VSVO-BBDF of order

Linear Multistep Method (LMM) given in [15] is given in the definition below.

Definition 1. The linear  -step method can be represented in standard form by an equation , where and , coefficients ,   are suitably chosen constants subject to conditions , , and is defined as the order of the particular method employed. This method is said to be explicit if and implicit otherwise.

Substituting and gives the corrector for the first and second points, respectively. Therefore by letting , and . This produced the following coefficients (Tables 4, 5, and 6), as in Definition 1 for the first and second points of VSVO-BBDF method.

Table 4: Coefficients for the first and second point of VSVO-BBDF when , and .

Order

Table 5: Coefficients for the first and second point of VSVO-BBDF when , and .

Order

Table 6: Coefficients for the first and second point of VSVO-BBDF when , and .

Order

3. Order Conditions for General VSVO-BBDF

As similar to analysis for order of Linear Multistep Method (LMM) given in [15], we use the following to determine the order of VSVO-BBDF method.

Definition 2. The LMM [15] and the associated difference operator defined by

are said to be of order if , . The general form for the constant is defined as

Consequently, BBDF method can be represented in standard form by an equation , where and are by matrices with elements and for . Since VSVO-BBDF for variable order is a block method, we extend the Definition 2 in the form of

And the general form for the constant is defined as

is equal to the coefficients of , where and .

4. Implementation of VSVO-BBDF Method

Throughout this section, we illustrate the effect of Newton-type scheme which in a general form of

The general form of VSVO-BBDF method is

with and are the back values. By setting,

Equation (15) in matrix-vector form is equivalent to

Equation (17) is simplified as

Newton iteration is performed to the system ; by taking the analogous form of (14) where is the Jacobian matrix of with respect to . Equation (14) is separated to three different matrices denoted as

Two-stage Newton iteration works to find the approximating solution to (1) with two simplified strategies based on evaluating the Jacobian () and LU factorization of [16]. Two-stage Newton iteration is carried out as follows.

Step 1. Compute the values for , where

Step 2. Calculate the corrected value for with the value from  (1).

Step 3. Solve for the second iteration.

Step 4. Let be equal to the second iteration of .

The error is defined as

and the maximum error is defined as

where TS is the total steps and is the number of equations. Consequently, the rest of the method is carried out as the appendix.

5. Stability Conditions for VSVO-BBDF Method

To begin this section, we provide some definitions for the stability conditions of VSVO-BBDF method as in (15).

Definition 3. A method is said to be zero stable if the roots of the polynomial satisfy the condition ,  .

Definition 4. A method is said to be absolutely stable in a region for a given if for that all the roots of the stability polynomial satisfy ,  .

Definition 5. A method is said to be stiffly stable if and are contained in the absolute stability region, and it is accurate for all when applied to the scalar test equation ; is a complex constant with , where , , and ,  , and are positive constants.

The stability polynomial, associated with the method of (15), is given by , while the absolute stability region of this method in the plane is determined by solving . Consequently, to determine zero stable, we substitute to for each order of VSVO-BBDF method.

Hence, the roots for the three different step size and order selections obtained by using Maple are listed in Table 7.

Table 7: Roots for different step sizes and orders.

Roots for different step size

, and
, and
,
,

,
,

and
and

, and
, and
,
,

,
,

and
and

, and
, and
,
,

,
,

and
and

The stability region was given by the set of points determined by the boundary ,  . We obtained the stability region by finding the region for which . Since all of the roots in Table 7 have modulus less than or equal to 1, the method (15) is said to be zero stable. Figure 2 shows the stability for orders  ,  , and    of VSVO-BBDF method, respectively. The stability regions lie outside the closed region for each case.

(a)

(b)

(c)

Figure 2: Stability regions for order (a)  , (b)  , and (c)  .

Based on Figure 2, VSVO-BBDF method possesses the region absolute stability, which contains almost the whole of the half-plane .

6. Choosing Order and Step Size

To increase the efficiency in BBDF algorithm, VSVO-BBDF algorithm is designed to have the capacity to vary not only the step size, but also the order of the method employed. The importance of choosing the step size is to achieve reduction in computation time and number of iterations. Meanwhile changing the order of the method is designed for finding the best approximation. The essential components of VSVO-BBDF algorithm are the local truncation error (LTE), strategies for deciding when to change step size and order, and techniques for changing the step size and order. Strategies proposed in [17] are applied in this study for choosing the step size and order. The strategy is to estimate the maximum step size for the following step. Methods of order ,  , and    are selected depending on the occurrence of every successful step. Consequently, the new step size is obtained from which order produces the maximum step size.

The user initially will have to provide an error tolerance limit, TOL, on any given step and obtain the LTE for each iteration. The LTE is obtained from

where is the th order method and is the th order method. By finding the LTEs, the maximum step size is defined as

where is the step size from previous block and is obtained from the maximum step size given in the above equations.

There are 3 cases of the possibilities in choosing the step size.

Case 1. From order 3 to order 4 and from order 4 to order 3, the step size is allowed to increase, decrease, or be maintained, that is, , or .

Case 2. From order 3 to order 5 and from order 4 to order 5, the step size is allowed to decrease, or be maintained, that is, , or .

Case 3. From order 5 to order 3 and from order 5 to order 4, the step size is allowed to decrease, or be maintained, that is, or .

The successful step is dependent on the condition . If this condition fails, the values of  ,   are rejected, and the current step is reiterated with step size selection . On the contrary, the step size increment for each successful step is defined as

  and if  then ,

where is the safety factor, is the order of the method while and is the step size from previous and current block, respectively. In this paper, is set to be 0.8 so as to make sure the rejected step is being reduced.

7. Numerical Results

We carry out numerical experiments arising from problems in physics to compare the performance of VSVO-BBDF method with stiff ODE solvers in MATLAB version 7 in Windows XP or later. The syntax for using ode15s and ode23s is similar for every test problem. For example,

where odefun is a function stored that evaluates the right side of the differential equations, int is a vector of the integration interval, , init is a column vector of initial conditions, , and options is an adjustable format of optional parameters to change the default integration properties.

We created the above option structure to vary the relative tolerance (RelTol) and absolute tolerance (AbsTol) according to specific error tolerance. These test problems are performed under different conditions of error tolerances—(a) 10⁻², (b) 10⁻⁴, and (c) 10⁻⁶.

The test problems and solutions are listed below.

Problem 1. Consider,

with solution

Problem 2. Consider,

with solution

See Kaps [18].

Problem 3. Consider,

with solution

See Lambert [15].

This paper considers the comparison of four different factors, namely, number of steps taken, average error, maximum error, and computation time. From Table 8, among the three methods tested, our method, VSVO-BBDF method, requires the shortest execution time, smallest maximum error, and average error for each given tolerance level. From Figure 3, we can see more clearly that VSVO-BBDF gives the lowest maximum error for every tolerance level.

Table 8: Numerical results for Problem 1.

TOL
MTD
TS
AVEE
MAXE
TIME

10⁻²
VSVO-BBDF
21

0.0103

ode15s
28

0.0313

ode23s
19

0.1406

10⁻⁴
VSVO-BBDF
48

0.0105

ode15s
60

0.0156

ode32s
42

0.0313

10⁻⁶
VSVO-BBDF
164

0.0115

ode15s
100

0.0313

ode23s
143

0.0469

(a)

(b)

Figure 3: Efficiency curves for Problem 1.

Again by comparing the four factors mentioned earlier, we can see VSVO-BBDF in Table 9 gives the minimum maximum error for every tolerance level except for TOL  . However, our method prevails in terms of average error for each given tolerance level. VSVO-BBDF once again requires the shortest execution time for each given tolerance level. Figure 4 illustrates the efficiency of VSVO-BBDF as compared to ode15s and ode23s.

Table 9: Numerical results for Problem 2.

TOL
MTD
STs
AVEE
MAXE
TIME

10⁻²
VSVO-BBDF
22

0.0106

ode15s
29

0.0781

ode23s
25

0.0781

10⁻⁴
VSVO-BBDF
54

0.0111

ode15s
55

0.1250

ode32s
118

0.2031

10⁻⁶
VSVO-BBDF
194

0.0143

ode15s
197

0.2344

ode23s
773

0.4531

(a)

(b)

Figure 4: Efficiency curves for Problem 2.

From Table 10, VSVO-BBDF method gives the shortest execution time for every given tolerance level. While in terms of average error and maximum error, VSVO-BBDF method once again gives the best result as compared to ode15s and ode23s. Figure 5 clarifies the efficiency of the proposed method based on its total steps and tolerance level.

Table 10: Numerical results for Problem 3.

TOL
MTD
TSs
AVEE
MAXE
TIME

10⁻²
VSVO-BBDF
35

0.0111

ode15s
45

0.0156

ode23s
137

0.2031

10⁻⁴
VSVO-BBDF
84

0.3281

ode15s
93

0.0469

ode32s
1211

0.0156

10⁻⁶
VSVO-BBDF
380

0.0199

ode15s
186

0.0781

ode23s
2829

2.3281

(a)

(b)

Figure 5: Efficiency curves for Problem 3.

8. Conclusion

The objective is met when VSVO-BBDF method outperformed ode15s and ode23s in terms of average error as well as maximum error. In most of the cases, VSVO-BBDF has successfully managed to reduce the number of total steps taken. As for the computation time wise, it gave lesser values for all cases. Therefore, we can conclude that VSVO-BBDF can serve as an alternative solver for solving stiff ordinary differential equations which arise in engineering and applied sciences.

Appendix

The following notation is used in the algorithm of VSVO-BBDF:

 : step size, : tolerance, TS: total steps, Order: order of the method, JCBN: jacobian, E: increment.

The algorithm for VSVO-BBDF code is given as follows.

Step 1. Let order = 3.

Step 2. Calculate initial step size and .

Step 3. Calculate initial array ,  .

Step 4. Calculate the predictor values for ,  ,  .

Step 5. Calculate for .

Step 6. Calculate the corrector values for ,  ,  .

Step 7. Calculate for .

Implementation of VSVO-BBDF.

Step 8. Calculate Jacobian (JACBN).

Step 9. Perform LU-Decomposition to calculate the increment .

Consider,

Test for convergence.

Step 10. If , convergence is false. Else go to Step 11.

Step 11. Calculate ,  .

Change order.

Step 12. Let order equal 4 and order equal 5.

Step 13. Repeat from Step 2 until Step 11 for each order correspondingly.

Step size and order control.

Step 14. Calculate for choosing the new step size.

Consider,

The final step size after a successful step () is given by

The final step size after a failure step () is given by

Step 15. If then .

Else .

Step 16. Update values:

Step 17. Repeat Step 1 to Step 16 until    = end point.

The Abbreviations Used in Figures  3–5 and Tables  8–10

TSs: | The total number of steps taken
TOL: | The initial value for the local error estimate
MAXE: | The maximum error
AVEE: | The average error
MTD: | The method used
TIME: | The total execution time (seconds).

References

W. H. Enright, T. E. Hull, and B. Lindberg, “Comparing numerical methods for stiff systems of O.D.E:s,” BIT, vol. 15, no. 1, pp. 10–48, 1975.

G. D. Byrne, A. C. Hindmarsh, K. R. Jackson, and H. G. Brown, “A comparison of two ode codes: gear and episode,” Computers and Chemical Engineering, vol. 1, no. 2, pp. 125–131, 1977.

R. Sacks-Davis, “Fixed leading coefficient implementation of sd-formulas for stiff ODES,” ACM Transactions on Mathematical Software, vol. 6, no. 4, pp. 540–562, 1980.

A. S. Mahmood, L. Casasús, and W. Al-Hayani, “The decomposition method for stiff systems of ordinary differential equations,” Applied Mathematics and Computation, vol. 167, no. 2, pp. 964–975, 2005.

J. Ibáñez, V. Hernández, E. Arias, and P. A. Ruiz, “Solving initial value problems for ordinary differential equations by two approaches: BDF and piecewise-linearized methods,” Computer Physics Communications, vol. 180, no. 5, pp. 712–723, 2009.

H. Aminikhah and J. Biazar, “A new analytical method for system of ODEs,” Numerical Methods for Partial Differential Equations, vol. 26, no. 5, pp. 1115–1124, 2010.

S. Abelman and K. C. Patidar, “Comparison of some recent numerical methods for initial-value problems for stiff ordinary differential equations,” Computers and Mathematics with Applications, vol. 55, no. 4, pp. 733–744, 2008.

L. F. Shampine and M. W. Reichelt, “The MATLAB ode suite,” SIAM Journal on Scientific Computing, vol. 18, no. 1, pp. 1–22, 1997.

L. F. Shampine, M. W. Reichelt, and J. A. Kierzenka, “Solving index-I DAEs in MATLAB and simulink,” SIAM Review, vol. 41, no. 3, pp. 538–552, 1999.

T. Kato and K. Ikeuchi, “Variable order and variable step-size integration method for transient analysis programs,” IEEE Transactions on Power Systems, vol. 6, no. 2, pp. 206–213, 1991.

Z. B. Ibrahim, M. Suleiman, and K. I. Othman, “Fixed coefficients block backward differentiation formulas for the numerical solution of stiff ordinary differential equations,” European Journal of Scientific Research, vol. 21, no. 3, pp. 508–520, 2008.

Z. B. Ibrahim, K. I. Othman, and M. B. Suleiman, “Variable step size block backward differentiation formula for solving stiff odes,” in Proceedings of the World Congress on Engineering, vol. 2, pp. 785–789, London, UK, 2007.

S. A. M. Yatim, Z. B. Ibrahim, K. I. Othman, and F. Ismail, “Fifth order variable step block backward differentiation formulae for solving stiff ODEs,” Proceedings of World Academy of Science, Engineering and Technology, vol. 62, pp. 998–1000, 2010.

S. A. M. Yatim, Z. B. Ibrahim, K. I. Othman, and M. B. Suleiman, “Quantitative comparison of numerical method for solving stiff ordinary differential equations,” Mathematical Problems in Engineering, vol. 2011, Article ID 193691, 12 pages, 2011.

J. D. Lambert, Computational Methods in Ordinary Differential Equation, John Wiley & Sons, New York, NY, USA, 1991.

K. R. Jackson, “The numerical solution of large systems of stiff IVPs for ODEs,” Applied Numerical Mathematics, vol. 20, no. 1-2, pp. 5–20, 1996.

G. Hall and J. M. Watt, Modern Numerical Methods for Ordinary Differential Equations, Clarendon Press, Oxford, UK, 1976.

P. Kaps and G. Wanner, “A study of Rosenbrock-type methods of high order,” Numerische Mathematik, vol. 38, no. 2, pp. 279–298, 1981.