ISRN Applied MathematicsVolume 2013 (2013), Article ID 567451, 10 pageshttp://dx.doi.org/10.1155/2013/567451

Research Article

An Accurate Block Solver for Stiff Initial Value Problems

H. Musa,¹ M. B. Suleiman,² F. Ismail,¹ N. Senu,¹ and Z. B. Ibrahim¹

¹Department of Mathematics, Faculty of Science, Universiti Putra Malaysia, 43400 Serdang, Selangor, Malaysia²Institute for Mathematical Research, Universiti Putra Malaysia, 43400 Serdang, Selangor, Malaysia

Received 19 June 2013; Accepted 20 July 2013

Academic Editors: S. Blanes and L. You

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.

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 Runge-Kutta (SDIRK), and singly implicit Runge-Kutta 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 A-stable 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 2-point 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.

Figure 1: Stability region of the 2-point VSSBBDF when .

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.

Figure 2: Stability region of the 2-point VSSBBDF when .

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.

Figure 3: Stability region of the 2-point VSSBBDF when .

A-stability 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.

Table 1: Numerical results for problem (1) (the circuit problem).

TOL
Method
IST
FS
TS
MAXE

NBDF
23
0
23

VSBBDF
63
0
63

NVSBBDF
43
0
43

VSSBBDF
44
0
44

NBDF
121
0
121

VSBBDF
77
0
77

NVSBBDF
62
0
62

VSSBBDF
70
0
70

NBDF
365
0
365

VSBBDF
107
0
107

NVSBBDF
101
0
101

VSSBBDF
142
0
142

Table 2: Numerical results for problem (2) (the cosine problem).

TOL
Method
IST
FS
TS
MAXE

NBDF
126
38
164

VSBBDF
72
4
76

NVSBBDF
65
1
66

VSSBBDF
72
0
72

NBDF
249
70
319

VSBBDF
159
7
166

NVSBBDF
116
0
116

VSSBBDF
149
0
149

NBDF
474
100
574

VSBBDF
457
19
476

NVSBBDF
398
4
402

VSSBBDF
472
0
472

Table 3: Numerical results for the nonlinear problem (3).

TOL
Method
IST
FS
TS
MAXE

NBDF
65
4
69

VSBBDF
45
0
45

NVSBBDF
48
1
49

VSSBBDF
49
0
49

NBDF
107
7
114

VSBBDF
63
0
63

NVSBBDF
70
4
74

VSSBBDF
80
0
80

NBDF
134
10
144

VSBBDF
99
0
99

NVSBBDF
110
8
118

VSSBBDF
169
0
169

Table 4: Numerical results for problem (4).

TOL
Method
IST
FS
TS
MAXE

NBDF
48
8
56

VSBBDF
46
0
46

NVSBBDF
48
0
48

VSSBBDF
50
0
50

NBDF
100
15
115

VSBBDF
66
0
66

NVSBBDF
73
0
73

VSSBBDF
96
0
96

NBDF
167
19
186

VSBBDF
117
1
118

NVSBBDF
142
0
142

VSSBBDF
224
0
224

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.

Figure 4: Graph of versus TOL for problem (1).

Figure 5: Graph of versus TOL for problem (2).

Figure 6: Graph of versus TOL for problem (3).

Figure 7: Graph of versus TOL for problem (4).

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 A-stable 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.

J. Álvarez and J. Rojo, “An improved class of generalized Runge-Kutta methods for stiff problems. I. The scalar case,” Applied Mathematics and Computation, vol. 130, no. 2-3, pp. 537–560, 2002.

K. Dekker and J. G. Verwer, Stability of Runge-Kutta Methods for Stiff Nonlinear Differential Equations, vol. 984, North-Holland, Amsterdam, The Netherlands, 1984.

E. Hairer and G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic 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.

B. Gustafsson and W. Kress, “Deferred correction methods for initial value problems,” BIT Numerical Mathematics, vol. 41, no. 5, pp. 986–995, 2001.

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.

G. G. Dahlquist, “A special stability problem for linear multistep methods,” BIT Numerical Mathematics, vol. 3, pp. 27–43, 1963.

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.

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. 1-2, pp. 117–130, 2000.

H. Musa, M. B. Suleiman, and F. Ismail, “A-stable 2-point block extended backward differentiation formula for solving stiff ordinary differential equations,” AIP Conference Proceedings, vol. 1450, pp. 254–258, 2012.

Z. B. Ibrahim, K. I. Othman, and M. Suleiman, “Implicit -point block backward differentiation formula for solving first-order stiff ODEs,” Applied Mathematics and Computation, vol. 186, no. 1, pp. 558–565, 2007.

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.

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.

J. Williams and F. de Hoog, “A class of -stable advanced multistep methods,” Mathematics of Computation, vol. 28, pp. 163–177, 1974.

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.

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.

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.

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.

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