A New Generalized Taylor-Like Explicit Method for Stiff Ordinary Differential Equations

: A new generalised Taylor-like explicit method for stiff ordinary differential equations (ODEs) is proposed. The algorithm is presented in its component and vector forms. The error and stability analysis of the method are developed showing that it has an arbitrary high order of convergence and the L-stability property. Moreover, it is veriﬁed that several integration schemes are special cases of the new general form. The method is applied on stiff problems and the numerical solutions are compared with those of the classical Taylor-like integration schemes. The results show that the proposed method is accurate and overcomes the shortcoming of the classical Taylor-like schemes in their component and vector forms.


Introduction
Most of the mathematical models describing real-world phenomena are often stiff ordinary differential equations (SODEs). Such problems involve a wide range of temporal scales and solving these SODEs requires careful treatment. Efficient numerical and semi-analytical numerical methods for solving stiff problems must have good accuracy, wide region of stability and low computational effort. It is well known that explicit linear multistep methods are not absolutely stable. In addition, most of the implicit methods are absolutely stable and work adequately with SODEs, but they involve a higher computational load per step than the explicit methods  Wu and Xia [9][10][11] presented a vector form for A-stable explicit one step Taylor-like method for solving stiff ODE systems. Wu and Xia [9] proposed a vector form for the two low accuracy methods [5] to be applied on stiff ODE systems. Wu and Xia [10] proposed a vector form for the sixth-order Taylor-like explicit method [7] to be applied on stiff ODE systems. Wu and Xia [9,10] showed the superiority of the Taylor-like methods when formulated in vector form compared to their component form. More recently, Wu and Xia [11] derived a general form of the Taylor-like explicit method and derived its corresponding vector form. One of the main advantages of using Taylor-like methods is that the approximate solution is given as an arbitrary order piecewise analytical function defended on the sub-intervals of the whole integration interval. This property offers different facility for adaptive error control [19,20]. Moreover, the Taylor-like method [11] is an arbitrary high order A-stable method that avoids extremely small stepsizes during the integration procedure. To avoid the analytical computation of the successive derivatives involved in the Taylor-like methods, numerical differentiation [21,22], automatic differentiation [23], differential transformation [24][25][26] and Infinity Computer with a new numeral system [27][28][29][30][31][32] can be used. In fact, Taylor-like explicit methods [5,7,[9][10][11] have computational drawbacks with zero-component derivative or zero-vector norm in their component or vector forms, respectively. Moreover, the computational errors due to the round-off, particularly in the case of long time intervals or large time-step sizes, may lead to very small values of the derivatives or very large values of the derivatives' ratio at some grid points. Consequently, poor results are obtained. To overcome these limitations a new generalised Taylor-like explicit method for SODEs and stiff ODE systems is present. The method is developed both in its component and vector forms. The error and stability analysis of the method are presented. It is shown that the new method has an arbitrary order of convergence and the L-stability property. Indeed, many other integration schemes are essentially special cases of the proposed general form method. The method is applied on stiff test problems and the numerical results are compared with those in the literature. The results show that the proposed method is accurate and avoids the shortcoming of the classical Taylor-like explicit methods both in their component and vector forms.

Generalised Taylor-Like Method
Consider the initial-value problem given by: where Assume that the solution of (1) can be approximated byȳ where a, b and c n are unknowns and determined later from the local truncation error and stability analysis of the method. By considering Taylor's expansions of y(t) andȳ(t) about t j , we have: where y The approximation in (2) is constructed so that y j+1 andȳ j+1 agree and their derivatives up to order m + 1, that is Equating the coefficients of the same powers of h L , L = 0, 1, . . . , m + 1, in (5) by zero, results in the following system of nonlinear equations Solving the first and last equations in (6) results in Consequently, solving the remaining equations in (6) results in From (2), (7) and (8) we obtain the following integration algorithm:

Local Truncation Error
The local truncation error T j+1 is readily obtained from subtracting (9) from the Taylor series expansion in (3) and collecting terms in h It is clear that relation (9) has at least m + 1 order of accuracy.

Stability Analysis
In order to examine the integration scheme (9) for the stability, let us consider the differential equation, where λ is a complex constant and Re(λ) < 0. For this equation, Equation (9) results in Setting z = λh in the above equation, the amplification factor is given by R(z) = e z and thus we have obtained the following A-stable and L-stable method with at least (m + 1) order of accuracy (12) where k = 1, 2, . . . , m + 1.
The local truncation error of (12) is given by

Consistency
Subtracting y j from both sides of (12) and dividing the result by h leads to Taking the limit as h tends to zero, on both sides of (14), yields showing that the method defined in (12) is consistent. Thus we have the following theorem: The generalised Taylor-like explicit method (12) is L-stable and convergent with at least (m + 1) order of accuracy.
For a fixed value of m and different values of k (k = 1, 2, . . . , m + 1) the method formulated in (12) results in (m + 1) L-stable different integration schemes with the same order of accuracy O(m + 1).

Remark 1.
For m = 0 and k = 1 the Equation (12) is reduced to which is the first order L-stable method in [5,9].

Remark 2.
For m = 0 and k = 2 the Equation (12) is reduced to which is the second order L-stable methods in [5,9].

Remark 3.
For m = 5 and k = 6 the Equation (12) is reduced to which is the sixth-order L-stable method in [7].

Remark 4.
For m = 4 and k = 6 the Equation (12) is reduced to which is the sixth-order L-stable method in [9].
Remark 5. For m = 5 and k = 7 the Equation (12) is reduced to which is the seventh-order L-stable method in [8].

Remark 6.
When k is set to m + 2, Equation (12) is reduced to and the local truncation error (13) is reduced to that is the classical general form of Taylor-like m + 2 order L-stable method in [11].

Extension to Vector Form
The method (12) can also be extended directly to systems of ODEs by using the new definitions of vector product and quotient defined in [9,12] where a, b ∈ C m and ρ = 1 The vector version of the proposed method is of the form The local truncation error of (25) is given by where ε is a vector and ε 2 = O(h m+3 ). The error in (26) is derived in a similar way as in [9,10].
The case k = m + 2 can be found in [11]. It is clear that the new method can overcome the restrictions in the classical Taylor-like explicit method in its component applicable form (i.e., y (m+2) j y (m+1) j = 0) and in its vector applicable form (i.e., . With the proposed method, the value of k can be adapted automatically to use any arbitrary pairs y . . , m + 1, while maintaining the order of convergence and the L-stability property.

Numerical Results
In this section, we provide six numerical experiments to illustrate the theoretical results obtained in Sections 2 and 3. All numerical experiments are carried out using Matlab 8.3. The test problems were collected from the literature and the results are compared in the follow-up. Problem 1. Consider the following SODE [8]: The theoretical solution of (27) is given by Problem 1 is solved numerically using the generalised Taylor-like method (GTL) in its component form. The results are shown in Tables 1-3. Table 1 compares the solution error presented in [8] using the seventh order Sin and Explicit Taylor-like methods (STL7 and ETL7) with one given by the new method of seventh order (GTL7) having m = 6 and k = 7. The results show that the GTL7 leads to a more accurate solution than the STL7 and ETL7. Table 2 lists the maximum solution error E max = max 0≤t i ≤0.5 |y(t i ) − y i | using the time-step size h = 0.05 at different values of m and k, k ≤ m + 1.
The results show that as the order m + 1, or the value of k, increases, the solution error decreases. Moreover, increasing k is more effective than increasing m for improving the solution accuracy, and setting k = m + 1 leads to more accurate results. Table 3 where ε is a parameter controlling the stiffness. The theoretical solution is given by From the theoretical point of view, we have y (k) (0) = 0 , k = 3, 5, 7, . . . , which results in the computational overflow when adopting the high order classical ETL method. Problem 2 is solved numerically at ε = 1/200 using the GTL component form with forth, fifth and sixth orders (i.e., the GTL4, GTL5, GTL6). The solution error, maximum solution error and the computed order of convergence are listed in Table 4. The results show that the GTL is not only stable and accurate but also overcomes the overflow in computation at t = 0 by reducing the value of k from 6 to 2, without losing the convergence order and the L-stability property. Figure 3 shows the exact and numerical solutions at h = 0.05 . Figure 4     In fact, since the GTL overcomes the overflow in computations by reducing the value of k, the accuracy may be little decreased while the order of convergence remains constant. Being that the effect occurs more in the component than in the vector form, we conclude that the component form GTL is more suitable for scalar stiff problems than for stiff ODE systems where the GTL vector form can be applied. Problem 3. Consider the following stiff ODE system [ [7,10,16]]: The theoretical solution is given by y 1 = 0.5 e −0.5t + e −20t (sin(20t) + cos(20t)) y 2 = 0.5 e −0.5t + e −20t (sin(20t) − cos(20t)) y 3 = 0.5 e −0.5t − e −20t (sin(20t) + cos(20t)) . (32) Due to the round-off error in computations, the ratio between y 6 j and y 5 j , j = 2, 3 may be very large at some grid points and the classical component form of the ETL6 [7,10] results in low accuracy, or even in an overflow, as shown in Table 5. On the other hand, the present component form of the GTL6 can overcome this drawback easily by reducing the value of k from 6 to 5 for h = {0.01, 0.1} at t = {0.55, 0.3}, respectively. Table 6 reveals that the vector forms of the ETL6 and GTL6 can overcome the overflow in computations with the component form of the ETL6 at h = 0.1 and both of them have accurate results since there are no vectors with zero norms. Figure 5 depicts the exact and numerical solutions at h = 0.1. Figures 6 and 7 show the solution error of different-order vector form GTL at h = 0.1 and h = 0.01 respectively.       Problem 4. Consider the following stiff ODE system [2,13,15].
The theoretical solution is given by Problem 4 is solved using both the component and vector forms of the ETL6 and GTL6. The numerical results are listed in Tables 7 and 8. Table 7 shows that the component form GTL6 results in more accurate results than the ETL6 for h = 0.005. Furthermore, the GTL6 can overcome the overflow in computations for h = 0.001 by reducing the value of k from 6 to 5 at t = 0.283. Table 8 shows that the vector form GTL6 leads to the same results as ETL6 for h = 0.01 due to the absence of zero-vector norms. For h = 0.005, ETL6 exhibits an overflow while GTL6 overcomes this overflow by reducing the value of k from 6 to 5 at t = 76.93. Figure 8 shows the exact and numerical solutions at h = 0.1. Figures 9-11 illustrate the solution error of different-order vector form GTL at h = 0.001, 0.005 and 0.1, respectively.
Tables 9 and 10 compare the results achieved by ETL6 [7,10] with those obtained by GTL6.
From the theoretical point of view, y 3 = 0 at t = 0.01 and so, as shown in Table 9, the component form ETL6 results in computational overflow, while THE GTL6 overcomes the overflow in computations at zero derivative component by changing the value of k from 6 to 5 at t = 0.01. Table 10 shows that both the ETL6 and GTL6 have the same accurate results in their vector forms due to the absence of zero-vector norms. Figure 12 shows the exact and numerical solutions at h = 0.01. Figure 13 depicts the solution error of different-order vector form GTL at h = 0.01.
This problem has no known analytical solution, and therefore, a numerical approximate solution using the built in BVP4C MATLAB solver [33][34][35] is taken as the reference solution for comparison. The problem is solved using different-order vector form GTL at h = 0.1 and the solution error is represented in Figure 14 at h = 0.1. The same results are obtained using the ETL due to the absence of zero derivatives or zero vectors. Moreover, the exact and numerical solutions are illustrated in Figure 15.

Conclusions
This paper presented a new generalised Taylor-like explicit method for SODEs and stiff ODE systems. The new method was formulated in its general component and vector applicable forms. The error and stability analysis of the method were discussed showing that it has at least (m + 1) order of convergence and the L-stability property. It was shown that several other integration schemes are essentially special cases of the proposed method.The algorithm is applied to six stiff test problems and the results are analyzed in several figures and tables and compared with those obtained by means of the classical Taylor-like methods. The results show that the new method is not only stable and accurate but also overcomes the shortcomings of the classical Taylor-like explicit methods in their component and vector forms by adapting the value of k during the computations without losing the convergence order and the L-stability property.

Conflicts of Interest:
The author declares no conflict of interest.