An accurate block solver for stiff initial value problems

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.


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][3][4].Also well known are methods based on Richardson extrapolation and deferred correction schemes [5][6][7][8].Others are multistep methods like the backward differentiation formula (BDF) (see, e.g., [4]).Initially, the BDF methods are known to be A-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 A-stable methods.Methods of this type are known as extended backward differentiation formula [10][11][12].Some block methods for stiff problems are also known to be A-stable and have order greater than two (see, e.g., [13][14][15][16][17][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  = 3/8 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.
Define the linear difference operator   associated with the method (3) by where  =  = 1 stands for the first point and  =  = 2 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,  3,1 is normalized to 1 and for the second point,  4,2 is normalized to 1; the following 2-point implicit block method is obtained as For optimal performance, we choose different parameters from those used in [15].We choose  = 3/8 and  = 1, 2 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  = 0 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.

Stability of the Method
This section investigates the stability regions of the methods ( 6), (7), and (8).When  = 1, 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 ( 2 −  − ) to obtain To show that the method is zero stable, we let ℎ = 0 in (15) and the following equation for the first characteristics polynomial is obtained as ISRN Applied Mathematics 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  = 1 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 A-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  = 2) as When ℎ = 0 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 A-stable.
Similarly, the stability polynomial for the method (8) (corresponding to  = 10/17) is given by and when ℎ = 0 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 A-stable.
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.

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  (+1) +1 denote the ( + 1)th iterate and Then (33) can be written as which is equivalent to where LU decomposition method is used to solve the system (35).

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  +1 +2 refers to the formula of order  + 1 and   +2 refers to the formula of order .
If the LTE is less than or equal to a specified error tolerance (TOL) (i.e., LTE ≤ TOL), the step size (from the previous block) is maintained as constant (equivalent to the formula when  = 1) 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 ( = 0.3) different from that in [15].If new ℎ > (1.7 × old ℎ), then the step size ℎ becomes ℎ = 1.7 × old ℎ (equivalent to the formula when  = 10/17).Otherwise if LTE > TOL, the step size is halved (equivalent to the formula when  = 2).

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  → 0. We choose  = 10 −3 .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  1 = −1 and  2 = −( −1 + 2).This allows the degree of stiffness to be regulated.We choose  = 10 −5 .(4) This physical problem is taken from [21] as follows: The exact solution is as follows: The eigenvalues are −2, −40 + 40, and −40 − 40.

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, 10 −2 , 10 −4 , and 10 −6 .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 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.

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  = 3/8, A-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 A-stable, but almost A-stable when the step ratio is 1.In addition, choosing  = 0 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.
, 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 3 :
Numerical results for the nonlinear problem (3).