A Comparison of Explicit Semi-Analytical Numerical Integration Methods for Solving Stiff ODE Systems

: In this study, a comparison among three semi-analytical numerical integration algorithms for solving stiff ODE systems is presented. The algorithms are based on Differential Transform Method (DTM) which are Multiple-Step DTM (MsDTM), Enhanced MsDTM (E-MsDTM) and MsDTM with Padé approximants (MsDTM-P). These methods can be classified as explicit one step semi-analytical numerical integration methods. The error and stability analysis of each method is presented. New important relationships among the methods are introduced. To demonstrate our results, a comparison of the accuracy, stability and computational efficiency of the methods is presented through solving some linear and nonlinear problems arising in applied science and engineering.


Introduction
Many mathematical problems arising from the real world cannot be solved completely by analytical means. One of the most important mathematical problems arising in applied science and engineering is stiff ODE systems. Stiff ODE systems, which have significantly different timescales of change, occur in many fields of engineering particularly in the studies of electrical circuits, chemical and biochemical reactions, optimal control, vibrations and fluid mechanics; see for example (Schiesser, 1993;El-Zahar, 2015). 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. It is well known that explicit numerical methods are intrinsically faster than commonly used implicit methods while implicit methods are commonly used for solving stiff problems because of their stability (Lambert, 1991;Hairer and Wanner, 1996;Butcher, 2008). In fact, the explicit numerical methods which recently are presented in (Wu, 1998;Wu and Xia, 2001;Novati, 2003;Ahmad et al., 2004;Ahmad and Yaacob, 2005;Ebady et al., 2012;Egbako and Adeboye, 2012;El-Zahar, 2013a;El-Zahar et al., 2014a) can also satisfy the stiff problems.
Although the treatment of HPM, ADM and VIM as multiple-step methods speeds up the convergence rate and improves the accuracy, these methods have some limitations which hamper their applications.
An analytical integration is required for each step, each solution series term in HPM and ADM and each iteration in VIM and so some difficulties may arise as follows: • The nonlinear and non-homogeneous parts of the ODE may be ill-conditioned such that the integration became very complicated • By increasing the number of iterations or the number of solution series terms, the number of terms of approximate solution may increase, so rapidly that the integration becomes both complicated and time consuming • "In VIM for some differential equations, we may have more than one Lagrange multiplier. Let L 1 and L 2 be different Lagrange multipliers for an ODE. We assume that L 1 is more accurate than L 2 and then by using L 1 the more fast the approximation to the exact solution can be obtained, but the computation will be complicated. The mentioned notation for L 2 is totally opposite" (Heydari et al., 2013;2014) • In ADM, for nonlinear problems, we need to evaluate the Adomian polynomials that mostly require tedious algebraic calculations • In HPM, we should suitably choose an initial guess, or infinite iterations are required Due to the above mentioned difficulties, the Multiple-step DTM (MsDTM), which does not require analytical integration, evaluation of the Lagrangian multiplier, or rather difficult computation for finding the Adomian polynomials, is one of the most effective multiple-step semi-analytical numerical methods. EL-Zahar, (2015) presented an adaptive step size MsDTM for solving singularly perturbed ODE systems, also known as stiff ODE systems. The results showed that the proposed method is an accurate and efficient method compared to classical MsDTM and RK4 method. However, for large values of the admissible local error or when a fixed step size is used, the problem of stability arises. A studying of this obstacle and how to overcome it is introduced in this study through presenting a comparison of the stability and accuracy of three semianalytical numerical differential transform based methods. These methods are MsDTM, Enhanced MsDTM (E-MsDTM) and MsDTM with Padé approximants (MsDTM-P). These methods can be classified as explicit one step semi-analytical numerical methods. The stability and error analysis of each method is presented. New important relationships among these methods are introduced. To demonstrate our results, a comparison of the accuracy, stability and computational efficiency of the methods is presented through solving some linear and nonlinear problems of practical importance in applied science and engineering.

Multiple-Step DTM
The basic definition and the fundamental theorems of the MsDTM are given in (Odibat et al., 2010;Keimanesh et al., 2011;Gokdogan et al., 2012b;Yildirim et al., 2012;Erturk et al., 2012;El-Zahar, 2012b;Khader and Megahed, 2014). For convenience of the reader, we present a review of the method as follows.
Consider the following non-linear initial-value problem: Let [t 0 , T] be the interval over which we want to find the solution of the initial value problem (1). In actual applications of the DTM, the Nth-order approximate solution of the initial value problem (1) can be expressed by the finite series (Zhou, 1986;Nik and Soleymani, 2013;Rashidi et al., 2013;2014;Al-Amr, 2014;Benhammouda et al., 2014): where: Equation (2) and (3) imply that the concept of differential transformation is derived from the Taylor series expansion. The following theorems can be deduced from (2) and (3) Theorem 6 Theorem 7 Using some fundamental operations of DTM, we have the following recurrence relation: where, F(k, Y(k)) is the differential transform of the function f(t, y(t)). The differential transform Y(k) of the unknown function y(k) can be obtained by solving the iterating algebraic system (4). In order to speed up the convergence rate and to improve the accuracy of resulting solutions, the entire domain [t 0 , T] is usually split into sub-intervals and the algorithm of MsDTM is applied as follows: Assume that the interval [t 0 , T] is divided into M subintervals [t m , t m+1 ], m = 0,1,…….M−1 of equal step size h = (T-t 0 )/M. The main ideas of the MsDTM are as follows. First, we apply the DTM to Equation (1) over the interval [t 0 , t 1 ] to obtain the following approximate solution: Using the initial conditions y 0 (t 0 ) = C 0 . For m≥1 and at each subinterval [t m , t m+1 ] we will use the initial conditions y m (t m ) = S m-1 (t m ) and apply the DTM to Equation (1) over the interval [t m , t m+1 ], where t 0 in Equation (5) is replaced by t m . The process is repeated and generates a sequence of approximate solutions y m (t), m = 0,1,…….M−1, for the solution y(t): Finally, the MsDTM assumes the following solution:

Error Analysis of MsDTM
The local truncation error of the approximate solution S m (t) is estimated by the inequality (Do and Jang, 2012;Bervillier, 2012;El-Zahar, 2013b): It is well known that stiff ODE systems are associated with the existence of eigenvalues with very large magnitudes which consequentially results in solutions usually having exponential behavior.
Therefore, if we let From the local truncation error (8), it is clear that, for a fixed order N, as the step size h = (t-t m ) tends to zero the approximate solution S m (t) tends to the exact solution y m (t) and thus the method is consistent, while numerically, from (9), for a fixed step size, the convergence of the MsDTM only begins when N is of order O(Kh).

Stability Analysis of MsDTM
In order to examine the MsDTM for the stability, let us consider the differential equation: where, λ is a complex constant and Re λ < . For this equation, the MsDTM solution, Equation (6), yields: Setting z = λh in Equation (12) results in the amplification factor: Which is the N th order Taylor expansion of e z and thus the MsDTM is not A-stable method.

Enhanced MsDTM
Since the approximate solution obtained by MsDTM may yield an increasing error due to the nonlinearity (Do and Jang, 2012), a modification of MsDTM, which is called Enhanced MsDTM (Do and Jang, 2012) is presented as a multi-stage integration method as follows: Suppose that we have solved the problem (1) up to a point t m-1 , m = 1,2,..,M-1 and have obtained a value S m-1 (t m ) as an approximation of y m-1 (t m ). Assuming the localization hypothesis (Lambert, 1991), y m (t m ) = S m-1 (t m ), we are interested in obtaining an approximate value for the true one y m (t m+1 ). For that purpose, the following method is developed.
The solution y m (t) verifies the following IVP: Equation (14) can be written as: We can replace y m (τ) in Equation (15) by the approximate solution S m (τ) from Equation (6) to obtain an improved approximation 1 ( ) m S t to y m (t): τ and thus we can update the value S m (t m+1 ) by: Then, (6). For each subinterval, repeat this process, Equation (17), to update the initial condition of the followed subinterval, Y m+1 (0) and consequently to obtain enhanced approximation to y m+1 (t).
The above process, Equation (17), is called Enhanced MsDTM with one iteration (E 1 MsDTM). To apply Enhanced MsDTM with two iterations, we solve the integral equation, Equation (16) and hence we get an enhanced initial condition In fact the above process, Equation (18), is a multiple-step successive approximations method and it is equivalent to multiple-step Picard iteration method (Djang, 1948;Lal and Moffatt, 1982;Wazwaz, 2011) with DTM solution as the zeroth approximation of y m (t).

Error Analysis of E-MsDTM
It is well known that (Innocentini, 1999;Youssef and El-Arabawy, 2007) when the function ( , ( )) f t y t satisfies Lipschitz condition with respect to y(t) in a rectangle R , irrespective of the choice of the initial approximation, the successive approximations to the solution of the problem (1). Also, the error of the approximate solution ( ) I m S t can estimated as follows: It is clear that the choice of the initial approximation, 0 i.e., N = 0, results in the following bounded error: which is the classical Picard iteration method error (Youssef and El-Arabawy, 2007). From the local truncation error (19), it is clear that, for a fixed number of iterations I, as the step size h = (tt m ) tends to zero the approximate solution ( ) S t y t → and thus the method is consistent, while numerically, for a fixed step size, the convergence of E I -MsDTM only begins when I is of order O(Lh).

Stability Analysis of E-MsDTM
In order to examine the stability of E-MsDTM with N th order and I iterations, we consider again Equation (10).
Using E 1 -MsDTM, the approximate solution of Equation (10) is given by: where, 0 ( ) m S t is the zeroth approximation of ( ) m y t and given by: From Equation (21) and (22) we have: for Equation (10), then Equation (23) becomes: and can be written as: which is the (N+1) th order approximate solution of (10) obtained by DTM. Replacing the solution in Equation (22) by that obtained in Equation (25) and repeating the above algorithm I times results in: which is the (N+I) th order approximate solution of Equation (10) obtained by DTM and thus the amplification factor of E I -MsDTM with N th order and I iterations is given by: which is the (N+I) th order Taylor expansion of e z and thus the E I -MsDTM is not A-stable method.
From the stability analysis of MsDTM and E-MsDTM, we get the following theorem.

Theorem 8
Enhanced-MsDTM with order N and iterations I and the MsDTM with order (N+I) have the same stability region.
Also, we observe that applying the above procedure from Equation (21) to (26) on Equation (10) in a vector form results in the following theorem.

Theorem 9
Enhanced-MsDTM with order N and iterations I is equivalent to the MsDTM with order (N+I) for the homogenous autonomous ODE system:

MsDTM-Padé Approximants
The goal of Padé approximants is to make the maximum error of MsDTM, Equation (8), as small as possible by maintaining the accuracy far outside the radius of convergence of the series solution (Ehle, 1973;Rashidi and Mohimanian Pour, 2010a;2010b;Lu, 2012;Kanth and Aruna, 2013;Domairry and Hatami, 2014) and to improve the stability property of the solution. Let the Padé approximation of the MsDTM series solution S m (t) on [t m , t m+1 ] is the quotient of two polynomials p m (t) and q m (t) of order p and q respectively and p+q = N (Ehle, 1973). We use the notation S m,pa (t) to denote the quotient: Let the polynomials used in (28) are: The polynomials in (28) are constructed so that S m (t) and S m,pa (t) agree at t = t m and their derivatives up to N. Since S m (t) has a Taylor expansion defined in Equation (6), the formal power series: Determines the coefficients of p m (t) and q m (t). Multiplying (30) by q m (t) gives: When the left side of (31) multiplied out and the coefficients of the powers of (t-t m ) k are set equal to zero for k = 0,1,….,N, the result is a system of N+1 linear equations in the N+1 unknown coefficients of p m (t) and q m (t). Solving this linear system we get S m,pa (t) rational approximation of S m (t). From the local truncation error (30), it is clear that, for a fixed order N, as (t-t m )→0, S m,pa (t)→S m (t) and thus the method is consistent.

Stability Analysis of MsDTM-P
For the stability of MsDTM with Padé approximants, it is well known that the Padé approximants of order (p, q) are A-stable if and only if 2 ≥ p-q ≥ 0 (Ehle, 1973;Hairer and Wanner, 1996, Theorem 4.12), that is, the main diagonal and two sub diagonals in Padé approximants are A-stable. The absolute stability of MsDTM with diagonal and two sub diagonals Padé approximants is shown in Fig. 2.

Numerical Experiments
To demonstrate our results, a comparison of the accuracy, stability regions and computational efficiency of the MsDTM, E-MsDTM and MsDTM-P is presented through solving some linear and nonlinear practical problems with known exact solutions arising in electrical and chemical engineering. All calculations are carried out by MAPLE 14 software in a PC with a Pentium 2 GHz and 512 MB of RAM. Figure 3 is a schematic diagram of a separately excited D.C motor where R the armature resistance; L the armature inductance; I the armature current; ν the terminal voltage of the motor; e the back electromotive force; T the load torque; J the torque of inertia and θ = ω ɺ is the motor speed. Taking

Example 1
as dimensionless variables, where K b is the constant of the torque, Ω is a nominal operating point, System (32) with parameters listed in Table 1 has been integrated using fourth order MsDTM Figure 4-6 show the stable and unstable solutions of (32) obtained by using 4MsDTM, E 1 -4MsDTM and 4MsDTM-P at different values of time-step size h. In each figure, we choose the time step and the plotting range so that the differences between the obtained solutions are clear. Figure 4 shows that the 4MsDTM solutions obtained at h = 0.0025, 0.0027 (h<h s4 ), converge to the true solutions, while the solution obtained at h = 0.0028 (h>h s4 ), is unstable and diverges far from the true solution. To overcome the stability restriction on using 4MsDTM with h = 0.0028 for (32), we should use smaller time-step size or higher order MsDTM, where the 5MsDTM has a larger stability region (h s5 = 3.2202e-003) than 4MsDTM. As shown in Fig. 5, although the E 1 -4MsDTM solutions, at h 0.0025, 0.003 (h<h s5 ), exhibiting spurious oscillations at the layer region (Hairer and Wanner, 1996;Shyy et al., 1992;Yost and Rao, 2000), they remain bounded and decrease with time and thus the solutions are stable. At h = 0.0035 (h>h s5 ), the E 1 -4MsDTM solutions oscillate wildly and quickly exit the range of the graph and therefore, there is a need to increase the order of the method or the number of iterations I, where E 1 -5MsDTM or E 2 -4MsDTM has a larger stability region (h s6 = 3.557e-003) than E 1 -4MsDTM. Figure 6 shows that the 4MsDTM-P solutions are stables and converge to the true solutions even for larger time-step size.
Results in Table 2 respectively, do not approximate the true solutions of (32) correctly. We insert dashes (-) to indicate this phenomenon. Obviously, these step sizes are not enough to meet the stability restrictions where 4MsDTM and E 1 -4MsDTM are not A-stable methods. and therefore, there is a need to decrease the time-step to obtain accurate solutions, while the A-stable 4MsDTM-P performs better even for large time step.
Results in Table 2 also show us that the E 1 -4MsDTM solutions are accurate better than the 4MsDTM solutions, while the 4MsDTM-P solutions are more accurate than E 1 -4MsDTM solutions Table 3 and 4 present a comparison of the maximum absolute error I E ω and the CPU used time, respectively, for solving (32) using the NMsDTM, E I -4MsDTM and NMsDTM-P at different orders, N, number of iterations, I and time-step size, h.
As shown from Table 3, the results of (I+4) MsDTM and E I -4MsDTM have the same accuracy and stability results, which confirms theoretical results in (Theorems 8,9), while the results of (I+4) MsDTM-P are more accurate. In addition, for a fixed order in NMsDTM or a fixed number of iterations in E I -4MsDTM, as the step size decreases the numerical solution converges to the exact one. While, for a fixed step size h, the convergence of (NMsDTM, E I -4MsDTM) only begins when (N, I) are (O(Kh), O(Lh)), respectively, where for (32), K = L = O(ε −1 ) ( [Vulanovic, 1989, Theorem 1]; (El-Zahar and EL-Kabeir, 2013, Lemma 3.2); Zhao and Xiao, 2010;Zhao et al., 2014). Table 4, show us that, for the same step size h, the (I+4) MsDTM has the lowest CPU usage and the E I -4MsDTM has the highest CPU usage in solving (32). Moreover, increasing the order N of NMsDTM is computationally cheaper than increasing the order of NMsDTM-P, while increasing the number of iterations I in E I -4MsDTM is computationally more expensive, where for I>1 successive analytical integrations are required (I-1)-times and therefore the MsDTM is the lowest computational cost while the E-MsDTM is the highest computational cost.

Example 2
Figure 7 describes a circular reaction with 3 substances A, B and C, with initial values A(0) = 1, B (0) = 2 and C(0) = 3. The system of ODEs describing the circular reaction in Fig. 7 is as follows: For system (33), the eigenvalues are 1.0110 03, 1.5964 01, 0.00 e e λ =− − and therefore the stiffness ratio is 63.33 and we have ( 4 = 2.7550e-003 s h , 5 = 3.1820e-003 s h ). The system has been integrated on the interval [0, 1] and the numerical results are presented in Fig. 8-10 and Table  5 and 6 using MsDTM, E-MsDTM and MsDTM-P at different orders, N, number of iterations, I and timestep size h. The present numerical results for system (33) provide support for the earlier analysis of the results presented for system (32), which confirm our theoretical results.
The exact solution of (34) is: where:    The problem has been integrated on the interval [0, 1] and the numerical results are presented in Table 7-9 using MsDTM, E-MsDTM and MsDTM-P at different orders, N, number of iterations, I and time-step size h.
Results in Table 7 show that, the E 1 -4MsDTM solutions are accurate better than the 4MsDTM solutions, while the 4MsDTM-P solutions are more accurate and have a good stability property.
Results in Table 8 show that, the E I -4MsDTM results are little better than the (I+4) MsDTM results, while the (I+4) MsDTM-P results are more accurate than the (I+4) MsDTM and E I -4MsDTM results.
Results in Table 4, 6, 9 show us that, the CPU usage in solving (34) is extremely higher than CPU usage in (32) and (33) and thus is due to the nonlinearity of (34). But we still have MsDTM with the lowest computational cost and E-MsDTM with the highest computational cost.

Conclusion
In this study, we have presented an effective comparison among three explicit one step semianalytical numerical methods for solving stiff ODE systems. The methods are the MsDTM, E-MsDTM and MsDTM-P. The error and stability analysis of each method is presented. The error analysis shows that the methods are consistent and so they closely match the given IVP (1) when the step size h is sufficiently small. The stability analysis shows that the E I -NMsDTM have stability regions larger than those with NMsDTM and both of them are conditionally stable methods, while NMsDTM-P is A-stable method and larger time steps can be used hence suitable for the stiff ODE systems. In addition, the (N+I) MsDTM and E I -NMsDTM have the same stability regions (Theorem 8) and are equivalent for linear homogenous autonomous initial-value ODE system (Theorem 9). Also, the present analysis of the E I -NMsDTM shows that it is equivalent to the multiplestep Picard iteration method with NDTM solution as the zeroth approximation of y m (t). We have applied the methods on three linear and nonlinear practical problems with known exact solutions and the numerical results are presented in figures and tables for comparison. The numerical results confirm our theoretical ones and additionally show that: • As the step size decreases or the order increases or the number of iterations for E-MsDTM increases, the accuracy of the methods increases • For the same order and step size, the E-MsDTM solutions are accurate better than the MsDTM From the view of computations, increasing the order of the MsDTM or MsDTM-P is computationally cheaper than increasing the number of iterations in E-MsDTM.
We conclude that the MsDTM-P is more suitable than classical MsDTM and E-MsDTM for solving stiff ODE systems because of its stability, higher accuracy and a relatively low computational cost.