Research Article  Open Access
H. Musa, M. B. Suleiman, F. Ismail, N. Senu, Z. B. Ibrahim, "An Accurate Block Solver for Stiff Initial Value Problems", International Scholarly Research Notices, vol. 2013, Article ID 567451, 10 pages, 2013. https://doi.org/10.1155/2013/567451
An Accurate Block Solver for Stiff Initial Value Problems
Abstract
New implicit block formulae that compute solution of stiff initial value problems at two points simultaneously are derived and implemented in a variable step size mode. The strategy for changing the step size for optimum performance involves halving, increasing by a multiple of 1.7, or maintaining the current step size. The stability analysis of the methods indicates their suitability for solving stiff problems. Numerical results are given and compared with some existing backward differentiation formula algorithms. The results indicate an improvement in terms of accuracy.
1. Introduction
Stiff problems occur in many fields of engineering science and they represent coupled physical systems having components varying with very different time scales [1]. Some physical examples include chemical kinetics where reactions go at very different speeds and in circuit simulation where components respond at widely different time scales. A considerable research has been done on numerical methods for solving stiff initial value problems (IVPs) of the following form: The most successful ones are implicit in nature. Examples include single step methods like the diagonally implicit Runge Kutta (DIRK) methods, single diagonally implicit RungeKutta (SDIRK), and singly implicit RungeKutta methods [2β4]. Also well known are methods based on Richardson extrapolation and deferred correction schemes [5β8]. Others are multistep methods like the backward differentiation formula (BDF) (see, e.g., [4]). Initially, the BDF methods are known to be stable only up to order 2 [9]. However, Cash in [10] proved that the barrier in [9] can be circumvented by adding extra future point to produce higher order stable methods. Methods of this type are known as extended backward differentiation formula [10β12]. Some block methods for stiff problems are also known to be Astable and have order greater than two (see, e.g., [13β18]).
One of the motivation of this paper is to develop a new block method based on the block backward differentiation formula [13, 14, 17]. The method is similar in fashion with that in [15] but differs greatly by the choice of the constant . Furthermore, our choice of gives rise to different formulae with better stability regions than in [15], even with the choice of the same step ratio . In the following sections of the paper, we will present the derivation of the method, its stability analysis, and implementation. Standard test problems will also be presented and a comparison of their numerical results with some known stiff solvers will be given.
2. Formulation of the Method
Consider a special case of the following BBDF given in [14, 17]: where = the variable step size used; ; = computed solution; , ; ; represents the formula for the first point; represents the formula for the second point; and = the value used in the step size changing strategy.
In this paper, a nonzero coefficient is introduced in (2) to have a formula of the following form: where .
Define the linear difference operator associated with the method (3) by where stands for the first point and stands for the second point. Expanding (4) using Taylorβs series and collecting like terms give a set of equations to be solved for and . For the derivation of the first point, is normalized to 1 and for the second point, is normalized to 1; the following 2point implicit block method is obtained asFor optimal performance, we choose different parameters from those used in [15]. We choose and and where the values of correspond to the step size control strategy of maintaining constant, halving, and increasing the step size by a multiple of 1.7.
When we have When we have And when we have Formulae (6), (7), and (8) compute two solution values concurrently at each step of the integration thereby having the advantage of reducing the total number of steps when a conventional nonblock method is used. Also, the formulae obtained are entirely different from those in [15]. It should be noted that when in (5), formula (3) reduces to (2) and we say that formula (2) is contained in (3). Hence, (3) is a superclass of (2) and contains (2) as a subclass. We will refer to our method as variable step size superclass of block backward differentiation formula (VSSBBDF). In the following section, the stability analysis of the method will be presented.
3. Stability of the Method
This section investigates the stability regions of the methods (6), (7), and (8). When , the formulae in (6) can be represented as: Collecting like terms we have
Let . Applying the scalar test equation, to (10) gives
Equation (12) is equivalent to where The stability polynomial is obtained by evaluating the determinant of to obtain To show that the method is zero stable, we let in (15) and the following equation for the first characteristics polynomial is obtained as Solving (16) for , we have The values of above indicate that the method is zero stable since no magnitude of the root is greater than one and the root is simple.
The plot of the stability polynomial (15) is given in Figure 1. The region of stability is the region outside the circular shape. It indicates that the method (6) is stable since the region covers the entire negative half plane.
The same procedure is applied to obtain the stability polynomial of the method (7) (corresponding to ) as When in (18), the first characteristic polynomial of the method (7) becomes and solving (19) for gives the following roots:
Again, it shows that the method (7) is zero stable.
The stability plot of the polynomial (18) is given in Figure 2 and the region outside the circle indicates that the method is also stable.
Similarly, the stability polynomial for the method (8) (corresponding to ) is given by and when in (21), the first characteristics polynomial is given by When (22) is solved for , the following values are obtained: Hence, the method (8) is also zero stable.
The stability plot for (21) is given in Figure 3. The method (8) is therefore almost stable.
Astability is a desirable property for numerical methods of solving stiff problems. From the forgoing analysis, the methods developed are suitable for solving stiff initial value problems. In the following section, the procedure involved in the implementation of the method will be explained.
4. Implementation of the Method
The implementation involves Newtonβs iteration. To simplify the code, methods (6), (7), and (8) are written in the following form: which is equivalent to the following matrix equation: Equation (25) can be represented as where
Let Newtonβs iteration takes the form Equation (32) is equivalent to We show the accuracy of the method by finding the maximum error in the following way. Let and be the computed and the exact solutions of (1), respectively. The absolute error is defined by
The maximum error is given by: where is the total number of steps and is the number of equations.
Let denote the th iterate and Then (33) can be written as which is equivalent to where LU decomposition method is used to solve the system (35). For the different values of , When When and when
4.1. Step Size Selection
The step size increase in this paper differs with that in [15] in the sense that an increase is made by a multiple of 1.7 of the previous step size. This is done in order to optimize total number of steps and computation time. Three techniques are employed in controlling the step size. At the beginning of the integration, an initial step size is determined and the local truncation error (LTE) is computed as follows: where refers to the formula of order and refers to the formula of order .
If the LTE is less than or equal to a specified error tolerance (TOL) (i.e., ), the step size (from the previous block) is maintained as constant (equivalent to the formula when ) and a new step size is computed as where is the order of the method and is a safety factor. Again we choose a safety factor () different from that in [15].
If , then the step size becomes (equivalent to the formula when ). Otherwise if , the step size is halved (equivalent to the formula when ).
5. Test Problems
The performance of the method is tested on the following stiff problems that are commonly found in engineering and physical sciences. (1) The circuit problem taken from [17] is as follows: βThe exact solution is as follows: (2) The cosine problem taken from [8, 15, 19] is as follows: βThe problem becomes increasingly stiff as . We choose .βThe exact solution is as follows: (3) We consider the system of nonlinear stiff IVP taken from [18, 20] as follows: βThe solution, which is independent of , is βThe eigenvalues are and . This allows the degree of stiffness to be regulated. We choose . (4) This physical problem is taken from [21] as follows: βThe exact solution is as follows: βThe eigenvalues are , , and .
6. Numerical Results
In the following are tables of results indicating the maximum error and the total steps taken to solve each of the problems given in the previous section. Each problem is solved at a prescribed tolerance, that is, , , and . The results are compared with some known successful BDF based algorithms like the nonblock BDF method (Gearβs method), the VSBBDF [14] method, and the NVSBBDF [15] method. Also given are the graphs of the performance profile for the methods.
The following notations are used in presenting the results in Tables 1, 2, 3, and 4:βTOL: Tolerance value usedβMethod: Name of the methodβIST: Number of successful stepsβFS: Number of failed stepsβTS: Total stepsβMAXE: Maximum error.




From Tables 1, 2, 3, and 4, the maximum error for the VSSBBDF method is less when compared to the NBDF, VSBBDF, and the NVSBBDF methods. This is also depicted clearly in the graphs shown in Figures 4, 5, 6, and 7. To compare between the methods, we can therefore conclude that the present scheme is more accurate. This is achieved by an optimal choice of the constant that gives better stability region and ensures accuracy at the same time. Furthermore, the total number of steps to complete the integration is competitive for all the problems tested.
7. Conclusion
A new variable step size algorithm based on block backward differentiation formula is developed for the integration of stiff problems occurring in the fields of engineering and physical sciences. The step size changing algorithm involves a strategy of increase by a multiple of 1.7, halving, or maintaining the current step size. By choosing , stable methods are obtained corresponding to the step size ratio of 1 and 2. This is an added advantage to the method because in most BDF algorithms, the method is not stable, but almost Astable when the step ratio is 1. In addition, choosing reduced the algorithm developed to the BBDF, thus making it its superclass. The results obtained show the suitability of the method for solving stiff physical problems and indicate better accuracy when compared with some existing BDF algorithms. Thus, the algorithm developed can be an alternative solver for initial value problems arising in engineering and physical sciences.
Acknowledgments
The authors thank the anonymous reviewers for their comments which improved the paper and the Institute for Mathematical Research, Universiti Putra Malaysia, for funding the research.
References
 L. Brugnano, F. Mazzia, and D. Trigiante, βFifty years of stiffness,β in Recent Advances in Computational and Applied Mathematics, pp. 1β21, 2011. View at: Google Scholar
 J. Álvarez and J. Rojo, βAn improved class of generalized RungeKutta methods for stiff problems. I. The scalar case,β Applied Mathematics and Computation, vol. 130, no. 23, pp. 537β560, 2002. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 K. Dekker and J. G. Verwer, Stability of RungeKutta Methods for Stiff Nonlinear Differential Equations, vol. 984, NorthHolland, Amsterdam, The Netherlands, 1984. View at: MathSciNet
 E. Hairer and G. Wanner, Solving Ordinary Differential Equations II: Stiff and DifferentialAlgebraic Problems, vol. 2, Springer, New York, NY, USA, 2004.
 K. Bohmer and H. J. Stetter, Defect Correction Methods: Theory and Applications, vol. 5, Springer, New York, NY, USA, 1984.
 R. Frank and C. W. Ueberhuber, βIterated defect correction for the efficient solution of stiff systems of ordinary differential equations,β BIT Numerical Mathematics, vol. 17, no. 2, pp. 146β159, 1977. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 B. Gustafsson and W. Kress, βDeferred correction methods for initial value problems,β BIT Numerical Mathematics, vol. 41, no. 5, pp. 986β995, 2001. View at: Publisher Site  Google Scholar  MathSciNet
 A. T. Layton and M. L. Minion, βImplications of the choice of quadrature nodes for Picard integral deferred corrections methods for ordinary differential equations,β BIT Numerical Mathematics, vol. 45, no. 2, pp. 341β373, 2005. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 G. G. Dahlquist, βA special stability problem for linear multistep methods,β BIT Numerical Mathematics, vol. 3, pp. 27β43, 1963. View at: Google Scholar  Zentralblatt MATH  MathSciNet
 J. R. Cash, βOn the integration of stiff systems of O.D.E.s using extended backward differentiation formulae,β Numerische Mathematik, vol. 34, no. 3, pp. 235β246, 1980. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 J. R. Cash, βModified extended backward differentiation formulae for the numerical solution of stiff initial value problems in ODEs and DAEs,β Journal of Computational and Applied Mathematics, vol. 125, no. 12, pp. 117β130, 2000. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 H. Musa, M. B. Suleiman, and F. Ismail, βAstable 2point block extended backward differentiation formula for solving stiff ordinary differential equations,β AIP Conference Proceedings, vol. 1450, pp. 254β258, 2012. View at: Publisher Site  Google Scholar
 Z. B. Ibrahim, K. I. Othman, and M. Suleiman, βImplicit $r$point block backward differentiation formula for solving firstorder stiff ODEs,β Applied Mathematics and Computation, vol. 186, no. 1, pp. 558β565, 2007. View at: Publisher Site  Google Scholar  MathSciNet
 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. View at: Google Scholar
 M. B. Suleiman, H. Musa, F. Ismail, and N. Senu, βA new variable step size block backward differentiation formula for solving stiff IVPs,β International Journal of Computer Mathematics, 2013. View at: Publisher Site  Google Scholar
 J. Williams and F. de Hoog, βA class of $A$stable advanced multistep methods,β Mathematics of Computation, vol. 28, pp. 163β177, 1974. View at: Google Scholar  MathSciNet
 S. A. M. Yatim, Z. B. Ibrahim, K. I. Othman, and M. B. Suleiman, βA quantitative comparison of numerical method for solving stiff ordinary differential equations,β Mathematical Problems in Engineering, vol. 2011, Article ID 193691, 12 pages, 2011. View at: Publisher Site  Google Scholar
 S. A. M. Yatim, Z. B. Ibrahim, K. I. Othman, and M. B. Suleiman, βA numerical algorithm for solving stiff ordinary differential equations,β Mathematical Problems in Engineering, vol. 2013, Article ID 989381, 11 pages, 2013. View at: Publisher Site  Google Scholar
 D. Kushnir and V. Rokhlin, βA highly accurate solver for stiff ordinary differential equations,β SIAM Journal on Scientific Computing, vol. 34, no. 3, pp. A1296βA1315, 2012. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 H. Aminikhah and M. Hemmatnezhad, βAn effective modification of the homotopy perturbation method for stiff systems of ordinary differential equations,β Applied Mathematics Letters, vol. 24, no. 9, pp. 1502β1508, 2011. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 J. D. Lambert, Computational Methods in Ordinary Differential Equations, John Wiley & Sons, New York, NY, USA, 1973. View at: MathSciNet
Copyright
Copyright © 2013 H. Musa 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.