A NON-LINEAR ABSOLUTELY-STABLE EXPLICIT NUMERICAL INTEGRATION ALGORITHM FOR STIFF INITIAL-VALUE PROBLEMS

The time-step in integration process has two restri ctions. The first one is the time step restriction due to accuracy requirement τac and the second one is the time-step restriction du e to stability requirement τst. The most of explicit methods have small stability regio ns and consequently small τst. It obliges us to solve stiff problems with small step size τst << τac. The implicit methods work well with stiff problem s but these methods require more work per step than the explici t methods. In this study, a non-linear absolutly st able explicit one step numerical integration algorithm i s proposed for solving non linear stiff initial-val ue problems in ordinary differential equations. The al gorithm is based on deriving a non-linear relation between the dependent variable and its derivatives from the well known Taylor expansion. The accuracy of the method depends on some unknown parameter inserted in Taylo r expansion and determined from the error analysis. The accuracy and stability properties of the method are investigated and shown to yield at least third order and A-stable. The results obtained in the numerical exp eriments show the efficiency of the present method in solving stiff initial value problems.


INTRODUCTION
The initial value problems with stiff ordinary differential equations occur in many fields of applied sciences, particularly in the studies of electrical circuits, vibrations, chemical reactions, biological economic systems. The problem of stiffness has been known for some time and has attracted the attention of many numerical analysts leading to surveys of methods for stiff problems. Explicit numerical integration methods are intrinsically faster than commonly used implicit methods. Implicit integration methods solve a system of equations for each solution step often requiring iterative solution methods (Burden and Faires, 2010;Lambert, 1991) to satisfy nonlinear algebraic equations. Solution of the equations is the main computational cost in the analysis of large ODE systems. Implicit methods are commonly used in solving stiff problems in ODEs because of their stability (Fatunla, 1982;Gupta et al., 1985;Butcher, 2000). Lack of stability causes the normally efficient explicit methods to be unsuitable for stiff problems but recently many authors introduced and developed explicit methods to solve stiff problems (Ahmad et al., 2004;Ahmad and Yaacob, 2005;Hairer et al., 1993;Lambert, 1973;Novati, 2003;Otunta and Ikhile, 1999;Egbako and Adeboye, 2012;Wu and Xia, 2001;2007;Niekerk, 1987;1988;Wu, 1998). Traditional linear multistep numerical integration methods, both explicit and implicit, are based on polynomial approximations in the time domain, while rational methods are based on rational function approximations. In this study, we introduce a new non-

AJAS
linear absolutely stable explicit one-step rational method that can be used to solve stiff problems effectively. The method is based on deriving a non-linear relation between the dependent variable and its derivatives from the well known Taylor expansion. The accuracy of the method depends on some unknown parameter inserted in Taylor expansion and determined from the error analysis. The accuracy and stability properties of the method are investigated and shown to yield at least third-order and A-stable. Some numerical examples are presented to illustrate the performance of the method.

MATERIALS AND METHODS
Consider the initial-value problem given by Equation (1): where it is assumed that fsatisfies all the requirements of the uniqueness theorem in order for (1) to have a unique solution. The interval [a, b] is divided into a number of subintervals [x j ,x j+1 ] with x 0 = a and x j = a + jh, such that h is the step size. Suppose that we have solved numerically the problem in (1) up to a point x j and have obtained a value y j as an approximation of y(x j ). Assuming the localization hypothesis (Lambert, 1991) where, φ is some parameter which determined later from the error analysis. From (2) and ( (5) Thus we get Equation (6): where, Now the parameter φ is determined from the local truncation error of the approximation (6).

Local Truncation Error
The local truncation error T j+1 is readily obtained from subtracting (6) from Taylor series expansion in (2) and collecting terms in h Equation (7): 2 j j j j 1 j (4) j j j j 3 4 5 2 j 2( 1)y y 3(y ) T , 12y 2( 1)y y y y h h O(h ). 24(y ) It is clear that the relation (6) has at least second order of accuracy.
From equating to zero the coefficient of h 3 in (7) we get Equation (8): and a local truncation error given by Equation (9): From (6) and (8) the numerical scheme is readily obtained, which may be written in the form Equation (10): 12h(y ) y y , y 0. 12(y ) 6hy y h (3(y ) 2y y ) After applying this scheme we will take as approximation for the true solution of (1) at x j+1 the value Science Publications AJAS y j+1 given by (10). Repeating the procedure along the nodes on the integration interval we will obtain a discrete solution for the problem in (1).

Consistency and Stability
Subtracting y j from both sides of (10) and dividing the result by h, leads to Equation (11): , y 0. h 12(y ) 6hy y h (3(y ) 2y y ) Taking limit as h tends to zero, on both sides of (11), we have Equation (12): Suggesting that the scheme defined by (10) is consistent.
In order to examine the present method for the stability, let us consider the differential equation: where, λ is a complex constant and Re(λ)<0. For this equation, Equation (10) can be rewritten as Equation (13): Setting z = λh in the above equation, the amplification factor is therefore Equation (14): 2 2 12 6z z R(z) . 12 6z z which is the (2, 2)-Pade' approximation to the exponential e z and thus the method is A-stable (Hairer et al., 1993). The stability region of the method consists of the left-half complex plane. Using MATLAB we plot the stability region for the method and the region as given in Fig. 1.

Remark 1
We observe that when the differential equation in (1) is autonomous, y f (y) ′ = and verifies the Equation (15): then 0 φ = and we obtain as a particular case of the present method, the scheme given by: that is, the first-order scheme in (Fatunla, 1982), which results to be exact when the differential equation in (1) verifies (15), i.e., the solution of the differential equation is of the form:

Remark 2
We observe that when the solution of the differential equation in (1) verifies the Equation (16): then φ = 0 and we obtain as another particular case of the present method, the scheme given by: f y ,f y and f 0, that is, the second-order scheme in (Lambert, 1973;Niekerk, 1987;1988), which results to be exact when the differential equation in (1) verifies (16), i.e., the solution of the differential equation is of the form: , p, q, r, s . r sx

Remark 3
We observe that the present method has at least fourth-order of accuracy when the solution of the differential equation in (1) is one of the following forms: where p,q,r,s,m are real constants. For these solution forms the coefficient of h 4 in the error Equation (9) is vanished. This occurs for example when the differential equation in (1) is a Reccati equation with constant coefficients, i.e., y′(x) = p+qy(x)+r(y(x)) 2 . In fact it may be verified easily that the present method is exactly the same with the fourth order method in (Niekerk, 1988) when the solution of the differential equation in (1) verifies the Equation (17):

RESULTS AND DISCUSSION
To test the proposed method, we have integrated four initial value problems. All the computations were done in MATLAB environment with double precision arithmetic on a personal computer with a Pentium based processor (2019 MHz).

A Stiff Equation
The present method was tested on the stiff problem taken from Equation (18)  The problem has been integrated on the interval [0,0.5] and the results are presented in Table 1 for different numbers of steps, NI. The errors have been defined as the maximum of the absolute errors on the nodal points in the integration interval: For comparison, the numerical solutions errors obtained using the third-order rational methods in (Lambert, 1973;Niekerk, 1987;1988) and the present method are shown in Table 2. Table 2 shows us that the numerical solutions obtained by the present method and the method in (Niekerk, 1988) approximate the true solution to great extend more than that obtained by the two methods in (Lambert, 1973;Niekerk, 1987).

A Stiff System
The above method may be also applied to a system of equations. If we have y,f(x,y)∈R m in (1), we have just to consider the formula in (10) for every component. Let be the stiff system taken from Equation Science Publications The exact solution is: The results for every component of the solution on the interval [0, 1] appear in Table 3. For comparison, the numerical solutions errors obtained using the third-order methods in (Lambert, 1973;Niekerk, 1987;1988) and the present method are shown in Table 4 and 5.

AJAS
In Table 4 and 5, the maximum absolute error generated by the third-order methods in (Lambert, 1973;Niekerk, 1987) show us that the numerical solutions do not approximate the true solutions of the stiff system correctly using 128 integration steps. We insert dashes (-) to indicate this phenomenon. This is because the thirdorder rational methods in (Lambert, 1973;Niekerk, 1987) are not A-stable methods. Obviously, 128 integration steps are not enough to meet the stability restrictions. Therefore, there is a need to increase the number of integration steps in order to obtain the accurate numerical solutions. Table 4 and 5 also show us that the third order method in (Niekerk, 1988), L-stable method, performs better compare to the present method for 256 integration steps, but not for 512 integration steps.

Singularly-Perturbed Problems
Singularly-perturbed IVPs usually contain a small parameter that multiplies the first-order derivative. These problems are characterized by the presence of thin layers where the solution varies very fast, whereas away from the layers the solution behaves regularly and varies slowly (Roos et al., 1996). These problems are stiff and their numerical treatment presents some major computational difficulties (Ramos, 2005;Kumar et al., 2009;El-Zahar, 2012;2013;El-Zahar and EL-Kabeir, 2013). The present method was tested on two singularly perturbed problems taken from (Ramos, 2005).
Consider the first non-linear singularly perturbed IVP given by Equation (20): which has exact solution: ( sin x / ) y(x) 1 . 1 e − ε = + This problem exhibits an initial layer near x = 0. Table 6 shows the numerical results obtained with the present method when the integration is performed on the interval [0, 1] for a small value of the perturbation parameter, ε = 0.01. Table 7 also shows us that the present method performs better compare to the methods in (Lambert, 1973;Niekerk, 1987;1988).
The second non-linear singularly perturbed IVP given by Equation (21): whose solution is: This problem exhibits an initial layer near x = 0. Numerical results are presented in Table 8 and 9 and show that the present method performs better in solving stiff singularly perturbed IVPs. Table 3. Error for problem 4.2 at different numbers of steps, NI.

CONCLUSION
In this article, we introduced a new non-linear explicit one-step numerical integration algorithm, which is simple to use and easy to implement, for solving stiff initial-value problems. The method is based on deriving a non-linear relation between the dependent variable and its derivatives from the well known Taylor expansion. The accuracy of the method depends on the parameter φ which inserted in Taylor expansion and determined from the error analysis. The result of the analysis shows that the method is A-stable and at least third order. These characteristics suggested their use for solving stiff initialvalue problems. The method can be classified as a rational one-step method. The results obtained in the numerical experiments show the efficiency of the present method in solving stiff and singularly perturbed IVPs at low cost as it is shown by the CPU values in the tables.

ACKNOWLEDGEMENT
The researcher would like to gratefully acknowledge the King Abdulaziz City for Science and Technology, Kingdom of Saudi Arabia, for granting and supporting for the project No. PC-11-0111.