A NON-LINEAR FIFTH ORDER A-STABLE EXPLICIT ONE-STEP METHOD FOR STIFF SYSTEMS ARISING IN CHEMICAL REACTIONS

In this paper, a non-linear explicit one step method is presented for solving stiff differential systems which characterize several kinds of linear reactions and diffusion from biochemistry, physiology, etc. The accuracy and stability properties of the method are investigated and shown to yield at least fifth-order and A-stable. The method is consistent and convergent. Some differential systems arising in chemical reactions are solved to illustrate the performance and accuracy of the method. Received: November 21, 2013 c © 2014 Academic Publications, Ltd. url: www.acadpubl.eu Correspondence author 342 E.R. El-Zahar, Y.S. Hamed, H.M. Habib AMS Subject Classification: 65L05


Introduction
Stiff problems are characterized by the fact that the numerical solutions of slow movements are considerably perturbed by nearby rapidly changing solutions, as it usually occurs in problems arising from chemical reactions.In fact, the first appearance of the term "stiff" was in the paper by Curtiss and Hirschfelder [3]on problems in chemical kinetics.Many phenomena of interest in physiology and biochemistry are characterized by reactions among several chemical species and diffusion in various mediums; see ( [14], [10], [12]).In a closed system, both reactions and diffusion are governed by the following system of differential equations.y ′ (t) = M y(t), which guarantees conservation of the total amount of y(t) for any t > 0.
The popular methods for differential reaction systems are simple explicit schemes such as Euler's method, Runge-Kutta methods, etc.However, it is wellknown that conditional stability, the typical weak point of explicit methods, is very fatal for stiff problems.In the past few decades, many studies on numerical methods for stiff ODEs have been done in various aspects ; see ( [15], [16], [13], [1], [12], [5], [4], [6]).
High order A-stable methods allowing large step sizes with stiff systems.A fourth-order A-stable method for numerically solving stiff differential systems arising in chemical reactions was presented in [4], and the method approximates the solution very well.In fact, extension of this method to fifth order method suffers some stability restrictions.Recently,a new non-linear explicit one-step method for solving stiff initial-value problems was presented in [6].The accuracy of the method depends on some parameter φ inserted in Taylor expansion and determined from the error analysis.The method is A-stable and at least third order.The reason of inserting a parameter φ in Taylor expansion is to overcome the stability restrictions in higher order integration schemes.
The aim of this paper is to propose a non-linear fifth-order A-stable explicit one step method for numerically solving stiff ODE systems which characterize several kinds of linear reactions and diffusion from biochemistry, physiology, etc.The method is based on the algorithm presented in [6].We prove absolute stability and convergence of the method.Some relations between the proposed method and that presented in [4] are derived.Three ODE systems arising in chemical reactions are solved to illustrate the performance and accuracy of the method

A Fifth Order A-Stable Method
Consider the following system of ODEs where , and it is assumed that f satisfies all the requirements of the uniqueness theorem in order for (1) to have a unique solution.The interval [0, b] is divided into a number of subintervals [ t j , t j+1 ] with t 0 = 0 and t j = jh, j = 1, 2, .. such that h is the step size.Suppose that we have solved numerically the problem in (1.b) up to a point t j and have obtained a value y j as an approximation of y ( t j ) .Assuming the localization hypothesis [11],y j = y ( t j ) , we are interested in obtaining an approximate value, y j+1 , for the true value y(t j+1 ).For that purpose, following the ideas in [6], we consider the following Taylor's expansions of y(t) about t j where φ is some parameter which will be determined later from the error analysis.From ( 2) and (3) we have and from ( 3) and ( 4) we have Thus we get where y j+1 ≃ y(x j+1 ), and it is assumed that 120y ′ j −60hy ′′ j +20h 2 y ′′′ j −5h 3 y (4) j + φh 4 y (5) j = 0. Now the parameter φ will be determined from the following 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), where collecting terms in hleads to It is clear that the relation ( 6) has at least fourth order of accuracy.
From equating to zero the coefficient of h 5 in ( 7) we get and a local truncation error given by From ( 6) and ( 8) the numerical scheme is readily obtained, which may be written in the form , y ′ j = 0, (10) where j = M 5 y j .After applying this scheme we will take as approximation for the true solution of (1) at t j+1 the value 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).

Convergence of the Method
We establish the convergence of the method by showing that the method is consistent and stable.
Subtracting y j from both sides of (10) and dividing the result by h, leads to Taking limit as h tends to zero, on both sides of (11), we have 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, Eq.( 10) can be rewritten as Setting z = λh in the above equation, the amplification factor is therefore which has modulus less than one on the left-half complex plane, and thus the method is A-stable [7].Using MATLAB we plot the absolute and relative stability regions for the method in figures 1 and 2. As shown in Figure 3 fingers do not overlap the imaginary axis and there are no poles in the left half-plane which confirms that the method is A-stable [8].This implies that the method is consistent and stable and hence the method is convergent.Remark 1.We observe that if we set φ = 0 in Eq.( 6), we obtain as a particular case of the present method, the scheme given by that is, the fourth-order A-stable scheme in [4] Remark 2. Applying the algorithm presented in [4], using fifth order Taylor series expansions results in a fourth order conditionally stable method, where at φ = 1, Eq. ( 6) becomes with local truncation error and stability function The absolute and relative stability regions for the scheme ( 15) are shown in figures 3 and 4.
And thus the extension of the method presented in [4] to fifth order method suffers some stability restrictions.

Application on Chemical Reaction Systems
In this section, we apply the present method on some chemical reaction systems.Figure 5 describes a circular reaction with 3 substancesy 1 , y 2 and y 3 , with initial values y 1 (0) = 1,y 2 (0) = 2 and y 3 (0) = 3.The system of ODEs is as follows For the stiff system (18), the stiffness ratio is 1.7809 E+17 and the drawback of explicit methods is more severe, where for (18) Euler's method is stable only if h ≤ 1.978E − 03, third order Runge-Kutta method (RKM) is stable only if h ≤ 2.845E − 03, classical fourth order RKM method is stable only if h ≤ 2.754E − 03 and the four-stage fifth-order RKM as in [2] is stable only if h ≤ 3.1819E − 03.The problem has been integrated on the interval [0, 3] and the results are presented in Table 1 for different step size, h.The errors have been defined as the maximum of the absolute errors on the nodal points in the integration interval E max = y(t j ) − y j 1 .

3691e-006
In Table 1, the maximum absolute errors generated by classical fifth order Runge-Kutta method, RK5, show us that the numerical solutions do not approximate the true solutions of the stiff system correctly using 1/256integration step size.We insert dashes (-) to indicate this phenomenon.This is because fifth order RKM is not A-stable method, which causes it to suffer some stability restrictions.Obviously, 1/256 integration step size is 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 while the new method performs better even for large step size.
The problem has been integrated on the interval [0, 10] and the results are presented in Table 2. Finally, we consider the chemical reaction defining the dehydration of talc [9]as shown in figure 8.The rates of the reactions are well described by first-order kinetics (rate proportional to amount) so the differential equations are: where y 1 (t) refers to talc, y 2 (t) refers to anthophyllite, y 3 (t)refers to enstatite, and y 4 (t) refers to quartz.Expressing the constituents in terms of volume and for a temperature of 830 • C and a pressure of 1000 bars, the coefficients are: a 11 = −2.592day −1 ,a 21 = 1.584 day −1 ,a 22 = −8.496e− 1day −1 , a 31 = 5.616e − 1 day −1 , a 32 = 6.912e − 1 day −1 , a 41 = 3.024e − 1 day −1 , a 42 = 7.2e − 2 day −1 .The stable solutions and the absolute errors result from present method are shown in figures 9 and 10.

Conclusions
In this article, we introduced a new non-linear fifth order explicit one-step method for numerically solving stiff ODE systems which characterize several kinds of linear reactions and diffusion from biochemistry, physiology, etc.Unlike most of explicit methods, the present method is absolutely stable in spite of its   explicitness.In addition to the method is of high order allowing large step sizes with stiff systems.Some important relations between the present method and that presented in [4] are presented in this paper.We have applied the method on three systems arising in chemical reactions and presented the numerical results in tables and figures.It can be observed that the numerical results verify the theoretical ones and the present method approximates the solution very well.

Figure 4 :
Figure 4: Relative stability and order stars for the method Eq. (15)

Figure 6 :
Figure 6: Stable solutions obtained using the present method (a) and the unstable solutions using fifth order RKM (b) at h = 1/256.

Figure 10 :
Figure 10: Absolute errors result from the present method at h = 1/8.

Table 2 .
Numerical results for system (19) using the present method