Table of Contents Author Guidelines Submit a Manuscript
Mathematical Problems in Engineering
Volume 2011, Article ID 193691, 12 pages
Research Article

A Quantitative Comparison of Numerical Method for Solving Stiff Ordinary Differential Equations

1Department of Mathematics, Faculty of Science, Universiti Putra Malaysia (UPM), Selangor, 43400 Serdang, Malaysia
2Department of Mathematics, Faculty of Computer and Mathematical Sciences, Universiti Teknologi MARA, Selangor, 40450 Shah Alam, Malaysia

Received 17 April 2011; Accepted 28 June 2011

Academic Editor: Claude Lamarque

Copyright © 2011 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.


We derive a variable step of the implicit block methods based on the backward differentiation formulae (BDF) for solving stiff initial value problems (IVPs). A simplified strategy in controlling the step size is proposed with the aim of optimizing the performance in terms of precision and computation time. The numerical results obtained support the enhancement of the method proposed as compared to MATLAB's suite of ordinary differential equations (ODEs) solvers, namely, ode15s and ode23s.

1. Introduction

Consider the first-order ordinary differential equations in the form of 𝑦=𝑓(𝑥,𝑦)(1.1) with given initial values 𝑦(𝑎)=𝑦0 in the given interval 𝑥[𝑎,𝑏]. The system of (1.1) is said to be stiff if the eigenvalues of the matrix 𝜕𝑓(𝑥,𝑦)/𝜕𝑦 have negative real parts at every time 𝑥 and varies greatly in magnitude. Solving stiff ODEs which aroused from mainly applied problems, for example, from physical, chemical, and biological phenomena is made easy with several methods of interest appeared in subroutine libraries [1]. Some examples can be seen in [24]. From the code pioneered by Gear called DIFSUB, the wide interest in stiff ODEs has caused many other codes evolved to meet the same objective of finding the most accurate approximation for IVPs [5]. Some renowned codes are EPISODE and LSODE [5]. With the advancements of the existing methods of solving ODEs, many of these codes are preinstalled in MATLAB to work as numerical solvers. An overview of MATLAB’s numerical solvers presented in [2] shows a wide choice of solvers according to the problem type, step, and order to deal with stiff and nonstiff problems. Since then, several measures have been brought up by many researchers to evaluate which method turns out to be the most efficient method [6]. Despite that there exist some limitations in each method, these methods proven to have their own strength in solving certain initial value problems [1].

The need to solve large systems of stiff IVPs has also influenced the development of other existing method. These include iterative linear equation solvers and sparse direct linear equation solvers. Further discussion and additional references on this method are discussed in more detail in [2]. The method that incorporated with BDF proposed by Gear has been expanded gradually to become an improved method [7]. The study on producing block approximations 𝑦𝑛+1,𝑦𝑛+2,𝑦𝑛+3,,𝑦𝑛+𝑘 also known as block backward differentiation formulae (BBDF) [7] was inspired by Gear’s method. BBDF method in [7] has verified the competency of computing concurrent solution values at different points. Apart from that, block approximations have been used in different methods to improve the methods to give a better accuracy and computation time. Consequently, the study by Ibrahim et al. in [8] is extended in a way that the accuracy is improved with reduction of total steps and lesser computational time.

In the next section, we discuss the derivation of the 5th-order variable step block backward differentiation formulae (5oBBDF). The strategies and implementations are presented next, followed by the results and discussion. It is also one of the main aims of this paper to prove that the method proposed acts as an alternative way in solving IVPs which arise in engineering and applied sciences. We are interested to compare the numerical results obtained with stiff ODEs solvers provided by MATLAB, namely, ode15s and ode23s.

2. General 5oBBDF Formulation

The ratio distance between current (𝑥𝑛) and previous step (𝑥𝑛1) is represented as 𝑞 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: 5th order block method of variable step size.

We find approximating polynomials 𝑃𝑘(𝑥) by means of a 𝑘-degree polynomial interpolating the values of 𝑦 at given points are (𝑥𝑛2,𝑦𝑛2),(𝑥𝑛1,𝑦𝑛1),,(𝑥𝑛+2,𝑦𝑛+2), 𝑃𝑘=𝑘𝑗=0𝑦𝑥𝑛+1𝑗𝐿𝑘,𝑗(𝑥),(2.1) where 𝐿𝑘,𝑗(𝑥)=𝑘𝑖=0𝑖𝑗𝑥𝑥𝑛+1𝑖𝑥𝑛+1𝑗𝑥𝑛+1𝑖foreach𝑗=0,1,,𝑘.(2.2) Define 𝑠=(𝑥𝑥𝑛+1)/ to find the 4th-order interpolating polynomial for (2.1), 𝑃𝑥(𝑥)=𝑃𝑛+1=+𝑠(2𝑞+1+𝑠)(𝑞+1+𝑠)(1+𝑠)(𝑠)𝑦4(𝑞+1)(𝑞+2)𝑛+2+(2𝑞+1+𝑠)(𝑞+1+𝑠)(1+𝑠)(𝑠1)𝑦(𝑞+1)(2𝑞+1)𝑛+1+(2𝑞+1+𝑠)(𝑞+1+𝑠)(𝑠)(𝑠1)4𝑞2𝑦𝑛+(2𝑞+1+𝑠)(1+𝑠)(𝑠)(𝑠1)𝑞2𝑦(𝑞+1)(𝑞+2)𝑛1+(𝑞+1+𝑠)(1+𝑠)(𝑠)(𝑠1)4𝑞2𝑦(2𝑞+1)(𝑞+1)𝑛2(2.3)𝑠 is equated to 0 and 1 after the polynomial is differentiated with respect to 𝑠, when 𝑥=𝑥𝑛+1and𝑥=𝑥𝑛+2, respectively, 𝑓𝑛+1=1+2𝑞2+3𝑞𝑦4(𝑞+1)(𝑞+2)𝑛+2+2+3𝑞𝑦(𝑞+1)(2𝑞+1)𝑛+1+12𝑞23𝑞4𝑞2𝑦𝑛+1+2𝑞𝑞2𝑦(𝑞+1)(𝑞+2)𝑛1+1𝑞4𝑞2𝑦(2𝑞+1)(𝑞+1)𝑛2,𝑓𝑛+2=10+3𝑞2+12𝑞𝑦2(𝑞+1)(𝑞+2)𝑛+2+4𝑞212𝑞8𝑦(𝑞+1)(2𝑞+1)𝑛+1+2+𝑞2+3𝑞2𝑞2𝑦𝑛+4𝑞4𝑞2𝑦(𝑞+1)(𝑞+2)𝑛1+𝑞+22𝑞2𝑦(2𝑞+1)(𝑞+1)𝑛2.(2.4) Upon substituting 𝑞=1,2 and 10/19 into 𝑃(𝑥𝑛+1) and 𝑃(𝑥𝑛+2), we obtained the coefficients for points 𝑦𝑛+1 and 𝑦𝑛+2 as follows: 𝑦for𝑞=1,𝑛+1=65𝑓𝑛+13𝑦10𝑛+2+95𝑦𝑛35𝑦𝑛1+1𝑦10𝑛2,𝑦𝑛+2=1225𝑓𝑛+1+48𝑦25𝑛+236𝑦25𝑛+16𝑦25𝑛13𝑦25𝑛2,𝑦for𝑞=2,𝑛+1=158𝑓𝑛+175𝑦128𝑛+2+225𝑦128𝑛25𝑦128𝑛1+3𝑦128𝑛2,𝑦𝑛+2=1223𝑓𝑛+2+192𝑦115𝑛+118𝑦23𝑛+3𝑦23𝑛12𝑦115𝑛2,𝑦for𝑞=10/19,𝑛+1=11311292𝑓𝑛+114703𝑦82688𝑛+2+1279161𝑦516800𝑛183027𝑦108800𝑛1+10469𝑦27200𝑛2,𝑦𝑛+2=13923095𝑓𝑛+2+89088𝑦40235𝑛+1242208𝑦77375𝑛+198911𝑦77375𝑛1658464𝑦1005875𝑛2.(2.5)

3. Implementation of 5oBBDF Method

It is commonly known that stiff ODEs codes must solve both nonlinear and linear systems at each step of differentiation. This is due to implicity of the formulae for solving stiff IVPs. Throughout this section, we illustrate the effect of Newton-type scheme to find the approximation solutions of 𝑦𝑛+1 and 𝑦𝑛+2 simultaneously in every step. The general forms of the 5th-order BBDF method are 𝑦𝑛+1=𝛼1𝑓𝑛+1+𝜃1𝑦𝑛+2+𝜓1,𝑦𝑛+2=𝛼1𝑓𝑛+2+𝜃1𝑦𝑛+1+𝜓2(3.1) with 𝜓1 and 𝜓2 are the back values. Equation (3.1) in matrix-vector form is equivalent to (𝐼𝐴)𝑌𝑛+1,𝑛+2=𝐵𝐹𝑛+1,𝑛+2+𝜉𝑛+1,𝑛+2.(3.2) By setting 𝐼=1001, 𝑌𝑛+1,𝑛+2=𝑦𝑛+1𝑦𝑛+2, 𝐴=0𝜃1𝜃20, 𝐵=𝛼100𝛼2, 𝐹𝑛+1,𝑛+2=𝑓𝑛+1𝑓𝑛+2, and 𝜉𝑛+1,𝑛+2=𝜓1𝜓2, (3.1) is simplified as 𝑓𝑛+1,𝑛+2=(𝐼𝐴)𝑌𝑛+1,𝑛+2𝐵𝐹𝑛+1,𝑛+2𝜉𝑛+1,𝑛+2=0.(3.3) Newton iteration is performed to the system 𝑓𝑛+1,𝑛+2=0, by taking the analogous form 𝑌(𝑖+1)𝑛+1,𝑛+2𝑌(𝑖)𝑛+1,𝑛+2=(𝐼𝐴)𝐵𝜕𝐹𝑌𝜕𝑌(𝑖)𝑛+1,𝑛+21(𝐼𝐴)𝑌(𝑖)𝑛+1,𝑛+2𝑌𝐵𝐹(𝑖)𝑛+1,𝑛+2𝜉𝑛+1,𝑛+2,(3.4) where 𝐽𝑛+1,𝑛+2=(𝜕𝐹/𝜕𝑌)(𝑌(𝑖)𝑛+1,𝑛+2) is the Jacobian matrix of 𝐹 with respect to 𝑌. Equation (3.4) is separated to three different matrices denoted as 𝐸(𝑖+1)1,2=𝑌(𝑖+1)𝑛+1,𝑛+2𝑌(𝑖)𝑛+1,𝑛+2,(3.5)𝐴=(𝐼𝐴)𝐵𝜕𝐹𝑌𝜕𝑌(𝑖)𝑛+1,𝑛+2,(3.6)𝐵=(𝐼𝐴)𝑌(𝑖)𝑛+1,𝑛+2𝑌𝐵𝐹(𝑖)𝑛+1,𝑛+2𝜉𝑛+1,𝑛+2.(3.7)

Two-stage Newton iteration is now introduced in order to find the approximate solution to (1.1). Thus, the corresponding linear system to be solved is 𝐴𝐸(𝑖+1)1,2=𝐵. According to Jackson in [9], evaluating the Jacobian 𝐽𝑛+1,𝑛+2 and LU factorization of 𝐴 require the most computation time during the operations. As a result, two strategies which are based on the step size are applied to ensure the efficiency of the method:(i)no calculation of new matrices 𝐴 and 𝐵 if the step size remains as previous step size. Hence, there is no new calculation of the Jacobian matrix 𝐽𝑛+1,𝑛+2,(ii)the matrices 𝐴 and 𝐵 are updated with new evaluation of the Jacobian matrix 𝐽𝑛+1,𝑛+2 for every occurrence of changing step size.

As a result, this paper agrees with [10] that this method has a higher tendency in saving computational time. Unlike most modern solvers that incorporated full Newton iterations, 5oBBDF method offers refined strategy for reevaluating the Jacobian matrix 𝐽𝑛+1,𝑛+2.

3.1. Stability Conditions for 5oBBDF

In this section, we provide the conditions for the stability of 5oBBDF method as in (2.5). To begin with, we include some definitions to support the practical criterion for a method to be useful in solving ODEs.

Definition 3.1. A method is said to be zero stable if the roots of the polynomial 𝜌(𝑧)=𝜋(𝑧,0) satisfy the condition 𝑧1=1|𝑧2|,𝑧21.

Definition 3.2. A method is said to be absolute stable in a region 𝑅 for a given 𝜆 if for that 𝜆, all the roots 𝑟𝑠 of the stability polynomial 𝜋(𝑟,𝜆)=𝜌(𝑟)𝜆𝜎(𝑟)=0 satisfy |𝑟𝑠|<1, 𝑠=1,2,,𝑘.

The stability polynomial, 𝑅(𝑡,), associated with the method of (2.5) is given by det(𝐴𝑡2𝐵𝑡𝐶), while the absolute stability region of this method in the 𝜆 plane is determined by solving det(𝐴𝑡2𝐵𝑡𝐶=0). Below are the absolute stability regions in the method of (2.5) for different step size selections: (𝑞=1)(𝑞=2) and (𝑞=10/19), respectively, 𝑅=𝑡,197𝑡125442𝑡254153𝑡12539𝑡252+72𝑡12542252𝑡125318𝑡12521+𝑅=125𝑡=0,𝑡,91𝑡464441𝑡1844173𝑡923289𝑡29442+45𝑡46421155𝑡73633𝑡9221+𝑅=2944𝑡=0,𝑡,1393273𝑡99968545298909𝑡3998740442149421𝑡1999370003544951521𝑡4209200002+393588𝑡999685421398292623𝑡3998740003753768𝑡13153752+47045881420920000𝑡=0.(3.8) To determine for zero stable, we substitute =𝜆=0 to (3.8). We will have the following equations for each step size selection mentioned previously: 197𝑡1254153𝑡12539𝑡252+1125𝑡=0,91𝑡464173𝑡923289𝑡29442+12944𝑡=0,1393273𝑡999685442149421𝑡1999370003544951521𝑡4209200002+47045881420920000𝑡=0.(3.9) Hence, the roots for three different step size selections obtained by using Maple are listed below(1)𝑡=.24414201370,𝑡=0.02079175991, and 𝑡=1,(2)𝑡=0.052708171410,𝑡=0.003257621961, and 𝑡=1,(3)𝑡=0.93455113330,𝑡=0.08581158625, and 𝑡=1.

Since all of the roots have modulus less than or equal to 1, the method (2.5) when (𝑞=1)(𝑞=2) and (𝑞=10/19) is zero stable.

The stability region was given by the set of points determined by the boundary 𝑡=𝑒𝑖𝜃,0𝜃2𝜋. The stability region is obtained by finding the region for which |𝑡|<1. Figure 2 shows the stability for the cases (𝑞=1)(𝑞=2) and (𝑞=10/19), respectively. The stability regions lie outside the closed region for each case.

Figure 2: Stability regions when 𝑞=1, 𝑞=2, and 𝑞=10/19.

Based on Figure 2, the 5th-order BBDF possesses the region absolute stability, which contains almost whole of the half-plane Re(𝜆)<0.

3.2. Choosing Step Size

The importance in choosing the step size is to achieve reduction in computation time and number of iterations. Therefore, this paper proposed three basic strategies for the step size selections. For each successful step, the step size remains constant (𝑞=1) or increased by a factor of 1.9 (𝑞=10/19). Inversely, when a fail step occurs, the next step size will be halved of the previous step size (𝑞=2). The user initially will have to provide an error tolerance limit, TOL on any given step, and obtain the local truncation error (LTE) for each iteration. The LTE is obtained from LTE=𝑦(𝑘+1)𝑛+2𝑦(𝑘)𝑛+2,𝑘=4,(3.10) where 𝑦(𝑘+1)𝑛+2 is the (𝑘+1)th order method, and 𝑦(𝑘)𝑛+2 is the 𝑘th order method.

The successful step is dependent on the condition LTE<TOL. If this condition fails, the values of 𝑦𝑛+1, 𝑦𝑛+2 are rejected, and the current step is reiterated with step size selection (𝑞=2). On the contrary, the step size increment for each successful step is defined as new=𝑐×old×TOLLTE1/𝑝.(3.11)

And if new>1.9×old, then new=1.9×old. Where 𝑐 is the safety factor, 𝑝 is the order of the method, while old and new are the step size from previous and current blocks, respectively. In this paper,𝑐 is set to be 0.8.

4. MATLAB’s Stiff Numerical Solvers

MATLAB offers several numerical solvers to solve either stiff or nonstiff ordinary differential equations. These built-in numerical solvers are capable in approximating solutions to almost any system of differential equations [2]. As far as this paper is concerned, we are interested in solving stiff ODEs. For comparison purposes, this paper considers only ode15s and ode23s. This is because both of the methods deal with stiff differential equations based on backward differentiation formulas (BDFs) and on a modified Rosenbrock formula of order 2, respectively. Therefore, a fair comparison is obtained, and the discussions are made easy.

5. Numerical Results

We carry out numerical experiments to compare the performance of the 5th-order BBDF method with stiff ODE solvers in MATLAB mentioned above. Three sets of stiff test problems aroused from problems in physics are tested for the purpose of elucidating the difference in numerical results obtained. These test problems are performed under different conditions of error tolerances—(a) 10−2, (b) 10−4, and (c) 10−6, and the error trend is defined by ||𝑦Error=approximate𝑦actual||.(5.1) The test problems and solution are listed below (1)basic circuit problem is 𝑦=20𝑦+24,𝑦(0)=0,0𝑥10,(5.2) with solution 𝑦(𝑥)=6/5(6/5)𝑒20𝑥,(2) oscillatory problems in physics: Torsion spring oscillator with dry and viscous friction [7] 𝑦1=998𝑦1+1998𝑦2,𝑦1𝑦(0)=1,0𝑥10,2=999𝑦11999𝑦2,𝑦2(0)=2,(5.3) with solution 𝑦1(𝑥)=2𝑒𝑥𝑒1000𝑥,𝑦2(𝑥)=𝑒𝑥+𝑒1000𝑥,(5.4) eigenvalues, 𝜆=1,1000,(3)Physical problem [11] is 𝑦1=20𝑦10.25𝑦219.75𝑦3,𝑦1𝑦(0)=1,0𝑥10,2=20𝑦120.25𝑦2+0.25𝑦3,𝑦2𝑦(0)=0,3=20𝑦119.75𝑦20.25𝑦3,𝑦3(0)=1,(5.5) with solution 𝑦1𝑒(𝑥)=0.50.5𝑥+𝑒20𝑥,𝑦(cos20𝑥+sin20𝑥)2𝑒(𝑥)=0.50.5𝑥𝑒20𝑥,𝑦(cos20𝑥sin20𝑥)3𝑒(𝑥)=0.50.5𝑥+𝑒20𝑥,(cos20𝑥sin20𝑥)(5.6) eigenvalues, 𝜆=0.5,20±20𝑖.

From Table 1, among the three methods tested, our method, 5oBBDF, requires the shortest execution time for each given tolerance level. On top of that, this method also gives the smallest maximum error and average error. From Figure 3, we can see that 5oBBDF gives the lowest maximum error for every tolerance level. When we look at the convergence of the methods compared in Figure 3, our method comparably converges faster with minimum number of steps taken.

Table 1: Numerical results for problem (1.1).
Figure 3: Efficiency curves for problem (1.1).

From Table 2, once again, 5oBBDF requires the shortest execution time for each given tolerance level. From the error trends described in Figure 4, we can see that the maximum errors for 5oBBDF are the smallest for each TOL as compared to ode15s and ode23s. So the convergence of the method, 5oBBDF is found to converge fastest among the methods tested.

Table 2: Numerical results for problem (2.1).
Figure 4: Efficiency curves for problem (2.1).

From Table 3, similarly, 5oBBDF gives the shortest execution time for each given tolerance level. This method is also completed with lesser total steps except when TOL is 10−2. The error trends for all of the methods tested are depicted in Figure 5. The error trend shows that 5oBBDF converged fastest with smallest maximum error for each TOL as compared to ode15s and ode23s.

Table 3: Numerical results for problem (2.2).
Figure 5: Efficiency curves for problem (2.2).

6. Conclusion

For all problems tested, it is proven that, the 5th-order variable step BBDF has outperformed the ode15s and ode23s in terms of average errors as well as maximum errors. It also managed to reduce the number of total steps taken in most of the cases especially when TOL is less than 10−2. When we tested the convergence of the methods compared in this paper, the error trend worked as the identifier for the fastest error of the method approaches zero. Figures 3, 4, and 5 illustrate how fast the error of 5th-order variable step BBDF converged as compared to the other two methods used. As a result, we found that this method converges fastest for all problems tested. In terms of computation timewise, it gave lesser values for all of the test problems. Therefore, we can conclude that the 5th-order variable step BBDF is one compatible alternative solver for solving stiff ordinary differential equations which mostly arise in engineering and applied sciences.


TS:The total number of steps taken
TOL:The initial value for the local error estimate
MAXE:The maximum error
AVEE:The average error
METHOD:The method used
TIME:The total execution time (seconds)
5oBBDF:5th-order variable step BBDF
𝑒:Error for each point of approximations.


  1. G. D. Byrne and A. C. Hindmarsh, “Stiff ODE solvers: a review of current and coming attractions,” Journal of Computational Physics, vol. 70, no. 1, pp. 1–62, 1987. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  2. C. F. Quo and M. D. Wang, “Quantitative comparison of numerical solvers for models of oscillatory biochemical systems,” in Proceedings of the 2nd International Multi-Symposiums on Computer and Computational Sciences (IMSCCS '07), pp. 44–51, University of Iowa, Iowa City, IA, USA, August 2007.
  3. C. W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice-Hall, Englewood Cliffs, NJ, USA, 1971.
  4. J. G. Verwer and D. Simpson, “Explicit methods for stiff ODEs from atmospheric chemistry,” Applied Numerical Mathematics, vol. 18, no. 1–3, pp. 413–430, 1995. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  5. 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. View at Google Scholar
  6. W. H. Enright, T. E. Hull, and B. Lindberg, “Comparing numerical methods for stiff systems of ODEs,” BIT Numerical Mathematics, vol. 15, no. 1, pp. 10–48, 1975. View at Publisher · View at Google Scholar
  7. Z. B. Ibrahim, M. B. 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. View at Google Scholar
  8. 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.
  9. 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. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  10. 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. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  11. J. C. Polking and D. Arnold, Ordinary Differential Equation Using MATLAB, Upper Saddle River, NJ, USA, Prentice Hall, 2nd edition, 1999.