L-stable and A-stable numerical method of order two for stiff differential equation

This article presents a finite difference method of order two for stiff differential equation that will overcome the effects of stiffness since it is A-stable. And, reflect the asymptotic behavior of the solution of the stiff problem since it is L-stable. The classical explicit Runge–Kutta methods of order two are not suitable for stiff problems since the numerical solution will not overcome the effects of stiffness and will not reflect the asymptotic behavior of the solution of the stiff problem (for example, Dahlquist problem). And, they are not both A-stable and L-stable. On applying Euler’s implicit method, the numerical solution will overcome the effects of stiffness and will reflect the asymptotic behavior of the solution of the stiff problem. And hence, it is both A-stable and L-stable. But it is of order one method. A new numerical method which is both A-stable and L-stable and of order two for the numerical solution of a stiff differential equation is presented in this article. It will overcome the effects of stiffness and reflects the asymptotic behavior of the solution. The new method works well for large values of step size and works well in large domain. The new method is a modified form of the mid-point rule for the stiff problems. The rate of order of convergence is proved to be two both theoretically and numerically. Experimental results show the performance of the method based on the metrics such as stability function, stability region, order star fingers, the rate of order of theoretical and numerical convergence, absolute and relative errors, percentage of numerical solution accuracy and local and global truncation errors both numerically and graphically.


Introduction
The problem considered in this article is the stiff differential equation It is a question of calculating the unknown curve which starts at t = a and satisfies the given stiff differential eq. (1). Aiken (1985) explained the notion of stiffness and computational aspects of stiffness. It is a collection of information given by many Scientists and Engineers on stiffness. Selvakumar k_selvakumar10@yahoo.com Numerical methods are tools allowing to approach solutions to problems which can have complex developments or which cannot be solved analytically. In engineering, engineers must tackle problems related to other disciplines such as the mechanics of structures or rocks, biology, chemistry or physics. Before solving these problems, it is important to define and adopt a rational framework.
In 1769-1770, Leonhard Euler in his book, Institutionum Calculi Integralis, derived a numerical method for the numerical solution of the problem (1) and it is of the form Atkinson Kendall (1980) where {a = t 0 , t 1 , t 2 , . . . , t N = b} is a sequence of points in [a, b], using the step size h, t i = t 0 + ih and h = t i+1 − t i , i = 0, 1, 2, . . . , N − 1. Method (2) is consistent with stiff problem (1) when the step size h approaches zero (h → 0).
In the literature, for stiff problems, it is proved that method (2) is of order one.
To design a numerical method, the designer must have a deeper knowledge of the stability of the differential equation. For a designer, some tips are given by Mazzia and Trrigiante (1994). Gear (1981) examined the development of numerical solving techniques from problem identification to the ever final preparation of automatic codes for the solution of classes of similar problems. It is a study of recent works which have advanced the state of the art or which offer the promise of the future. It addresses several issues that are only part of the development path.

Main results of this article
A new numerical method which is both A-stable and L-stable and of order two for the numerical solution of a stiff differential equation is presented in this article. It will overcome the effects of stiffness and reflects the asymptotic behavior of the solution. And hence, it is both A-stable and L-stable. The new method is a modified form of the mid-point rule for the stiff problems. The rate of order of convergence is proved to be two both theoretically and numerically. Experimental results show the performance of the method based on the metrics such as stability function, stability region, order star fingers, theoretical and numerical rate of the order of convergence, absolute and relative errors, percentage of numerical solution accuracy and local and global truncation errors both numerically and graphically. The new method works well for large values of step size and works well in large domain.
In brief, this article presents a new numerical method for the numerical solution of the first-order stiff differential eq. (1). In Sect. 2, a new numerical method is proposed for the numerical solution of the problem (1). In sect. 3, stability analysis using the amplification factor is given and in sect. 4, stability regions, stability function, order star fingers, order star finger region, relative order star region and relative stability regions for the stiff problem are presented. In sect. 5, the theoretical rate of the order of convergence is proved. In sect. 6, the numerical rate of the order of convergence is proved. Finally, in sect. 7, experimental results show the performance of the method both numerically and graphically.
Throughout this article, C is a constant independent of the nodal point i, and the step size h. The term classical used by Doolan et al., in Doolan et al. (1980) refers to the numerical methods derived from Taylor's series expansion which are not exponentially fitted. Explicit in the sense, the unknown term y i+1 does not occur on the right hand side of the difference equation. And so, no iteration will be involved to evaluate the solution y i+1 . Implicit in the sense, the unknown term y i+1 occurs on the right hand side of the difference equation. And so, iteration will be involved to evaluate the solution y i+1 for the convergence of the solution. The classical explicit Runge-Kutta methods of order two refer to Jain (1984) From this representation, we list three explicit Runge-Kutta methods of order two, namely 1. If b 1 = 1 2 , b 2 = 1 2 , c 2 = 1, and a 21 = 1, then one will get the explicit Runge-Kutta method of order two. It is also named Euler-Cauchy method or Improved Euler method. 2. If b 1 = 0, b 2 = 1, c 2 = 1 2 , and a 21 = 1 2 , then one will get the explicit Runge-Kutta method of order two. It is also named as mid-point rule. 3. If b 1 = 1 4 , b 2 = 3 4 , c 2 = 2 3 , and a 21 = 2 3 , then one will get the explicit Runge-Kutta method of order two. It is also named as Heuns method.

Numerical method
A new numerical method proposed for the numerical solution of the stiff problem (1) is where {a = t 0 , t 1 , t 2 , . . . , t N = b} is a sequence of points in [a, b], using the step size h, t i = t 0 + ih and h = t i+1 − t i , i = 0, 1, 2, . . . , N − 1. Method (3) is consistent with the problem (1) when the step size h approaches zero (h → 0). The Butcher tableau for the method (3) is This article focuses on the method which is fully stable on the entire right half of the complex plane since the exact solution of the Dahlquist problem is stable on the entire right half of the complex plane. The new method must be unconditionally stable. The new method should work well for large values of step size and work well in large domain. The method must be both A-stable and L-stable.

Stability result
The stability analysis of the stiff problem is presented using the amplification factor. On applying the method (3) to the Dahlquist problem y = −λy, t ∈ [0, 1], λ > 0, y(0) = 1 ( 4 ) whose exact solution is y(t) = ex p(−λt). The corresponding difference equation is The corresponding amplification factor Q(λh) is

Stability
A numerical method is stable if Using (7), the maximum step size for which the numerical method is stable can be determined. Equation (7) holds for all values of h in [0, 1] and no restriction on h for stability.
On applying the explicit method of Runge-Kutta of order two (Improved Euler method) where k 1 = h f (t n , y n ), k 2 = h f (t n + h, y n + k 1 ) to the Dahlquist problem (4), we have the amplification factor Q(λh) as The expression (7) will not hold for all values of h, but only for 0 ≤ λh ≤ 2. The numerical solution is stable only in a small region on the right half of the complex plane and the numerical solution will not reflect the asymptotic behavior of the solution. But, the exact solution of the Dahlquist problem is stable on the entire right half of the complex plane. So, the explicit method of Runge-Kutta is not suitable for the stiff problem (1). Again, on applying the explicit Runge-Kutta method of order two (the mid-point rule), 2 h, y n + 1 2 k 1 ) to the Dahlquist problem (4), we have the amplification factor Q(λh) as The expression (7) will not hold for all values of h, but only for 0 ≤ λh ≤ 2.
Again, on applying the explicit Runge-Kutta method of order two ( Heuns method ), where k 1 = h f (t n , y n ), k 2 = h f (t n + 2 3 h, y n + 2 3 k 1 ) to the Dahlquist problem (4), we have the amplification factor Q(λh) as The expression (7) will not hold for all values of h, but only for 0 ≤ λh ≤ 2.
It is observed that the explicit Runge-Kutta methods of order two (Improved Euler method, the mid-point rule and Heuns method) have the same amplification factor Q(λh) as It is observed that the numerical solution is stable only in a small region on the right half of the complex plane and the numerical solution will not overcome the stiffness and will not reflect the asymptotic behavior of the solution. Finally, on applying the implicit Trapezoidal method of order two (4), we have the amplification factor Q(λh) as The expression (7) will hold for all values of 0 ≤ λh ≤ 2. The numerical solution is not stable fully on the right half of the complex plane and the numerical solution will not reflect the asymptotic behavior of the solution. So order two implicit Trapezoidal method is not suitable for the stiff problem (1) in the sense that the numerical solution is not stable fully on the right half of the complex plane.
In this article, we have focused on the method which is fully stable on the entire right half of the complex plane since the exact solution of the Dahlquist problem is stable on the entire right half of the complex plane. And, the method should be both A-stable and L-stable. The above observations motivate to design the numerical method (3) to the stiff problem (1). The exact solution of the Dahlquist problem is stable on the entire right half of the complex plane is the main key point. And, hence the method (3) is designed so that it is unconditionally stable to the Dahlquist problem (4), on the full right half of the complex plane, that is, it will overcome the stiffness and will reflect the asymptotic behavior of the solution.

Stability region
The stability region, order star finger region, relative order star region and relative stability regions of the solution of the stiff problem are presented using stability function to show boundedness, convergence and order of numerical solution. Also, stability function is applied to show the numerical method (3) is A-stable and L-stable.

Stability region of the exact solution
From the exact solution of the Dahlquist problem (4), the stability function of the exact solution is given by where z = hλ, z is a complex number and the stability region is the full right half of the complex plane shown in Fig. 1. It represents the bounded region {z : z ∈ C, Re(z) ≥ 0}.

Stability region of the numerical solution
On applying a new method (3) to the problem (4), the stability function corresponding to the numerical solution is given by where z = hλ, z is a complex number and the corresponding stability region is a region outside the contour R1(z) = 1 + z+ z 2 2 in the complex plane which is presented in Fig. 2. It represents the bounded region {z : z ∈ C, abs(R1(z)) ≥ 1} = . It is observed that the stability region includes the entire right half plane and also, in the left half plane the region outside the region {z : z ∈ C, abs(R1(z)) ≤ 1}. R(z) is the (0,2) Pade approximation to exp(-z).

Order star finger
The number of fingers in the absolute stability region will be considered for the order of the numerical method. In Fig. 3, there are two fingers in the absolute stability region. The numerical method (3) is of order two since the number of fingers in the absolute stability region is two.

Order star finger region of the numerical method
The bounded region given in Fig. 4 is the order star finger region of the numerical method (3).

The relative stability region
The bounded region given in Fig. 5 is the relative stability region or the order star of the first kind for the numerical method.

Absolute relative stability region of the numerical method
For the method (4), absolute stability region is obtained relative comparison with one is given in Fig. 6. It represents the bounded region {z : z ∈ C, abs(R1(z)) ≤ 1}

Order graph
The degree of the stability function R(z) is the order of the method.
Result If R(z) = P(z) Q(z) and deg (P(z)) = k and deg (Q(z)) = j where k ≤ j then deg(P(z)) ≤ k + j. Using this result, the order of the method (3) is two since deg (R(z)) = 2.

A-stable
If the stability region of the solution of the numerical method is the full right half-plane then the numerical method is Astable. The stability region of the numerical solution of the numerical method (3) is the union of the full right half-plane and some region in the left half-plane from Fig. 2. And so, the numerical method is A-stable. The new method will overcome the effects of stiffness since it is A-stable.  The stability regions of explicit Runge-Kutta methods and implicit Trapezoidal method of order two are given in Figs. 7 and 8, respectively. From the figures, it is clear that both the methods are not covering the entire right half of the complex plane. In the case of explicit methods covers only a small region of stability in the neighborhood of the point (2, 0). In case of implicit Trapezoidal method, there is a gap in the neighborhood of the point (2, 0). Hence, the explicit Runge-Kutta methods (Improved Euler method, the mid-point rule and Heuns method) and implicit Trapezoidal method of order two are not A-stable since they cannot overcome the effects of stiffness.

L-stable
A-stable method is said to be L-stable if the stability function R(z) → 0 as z → ∞. On applying this condition to the numerical method (3), it is observed that the stability function The stability function R(z) for explicit Runge-Kutta methods (Improved Euler method, the mid-point rule and Heuns method) of order two is As z → ∞, R(z) → ∞ and so explicit Runge-Kutta methods of order two are not L-stable. The stability function R(z) for implicit Trapezoidal method of order two is As z → ∞, R(z) → −1 and so the implicit Trapezoidal method of order two is not L-stable. Hence, the explicit Runge-Kutta methods (Improved Euler method, the midpoint rule and Heuns method) of order two and implicit Trapezoidal method of order two will not reflect the asymptotic behavior of the solution of the stiff problem.

Theoretical rate of order of convergence
In this section, the theoretical rate of order of convergence of the numerical method (3) is proved to be two. The main result of this section is shown in the following theorem.
Theorem 1 If y(t) is the solution of the given differential eq.
(1) and y i is the numerical solution of the numerical method (3), then, for i=0(1)N, we have an error estimate of the form, where C is independent of i and h.

Proof
The solution y(t) of (1) at the nodal point t = t i+1 can be expressed in terms of solution y(t) at the nodal point at t = t i , using Taylor's series as follows: where h = t i+1 − t i , i = 0, 1, . . . N − 1. The numerical method (3) can be expressed as where h = t i+1 −t i and y i = y(t i ), i = 0, 1, . . . N . Equation (11) can be rewritten as From eqs. (11) and (13), we have where e i = y(t i ) − y i , i = 0, 1, . . . N . The error per step, local truncation error (LTE) is given by The error at a given time t is termed as global truncation error (GTE) can be obtained from (15) by rewritting (15) as follows: From eq. (14), we have to find the error estimate for the absolute error. For i = 0, 1, 2,...,N-1, we have continuing like this finally, From eqs. (17-20), we have, Let k = max y (t j ) for j = 0(1)N − 1, then, Now for i = 0 and i = N, we have the estimate. And for i = i, 2,...,N-1, we shall find the error estimate by taking equation (19) we have, The right hand side of eq. (22) is same as in eq. (24), and so, for i = 1(1)N-1. Finally, from (17), (24) and (25) the desired estimate (10) follows for i=0 (1) Hence, the proof.
Theoretical rate of order of convergence If y(t) is the solution of the given differential eq. (1) and y i is the numerical solution of the numerical method. And, if we have an error estimate of the form,  (27) and (28), the theoretical rate of the order of convergence of the numerical method (3) is of two (n = 2). When the exact solution is not known how to find the numerical rate of the order of convergence is the next arising question. To answer this question in the next section, an alternate numerical rate of the order of convergence is derived.

Numerical rate of order of convergence
In this section, the numerical rate of the order of convergence of the numerical method (3) is proved as two. The main result of this section is stated in the following theorem.
where C is independent of i and h. (1)N. For i = 0, w 0 = 0, and for i = 1(1)N we have, from (3), where Using the procedure in Theorem. 1, eq. (30) gets the form From (31), the error per step is given by Equation (31) can be rewritten as, and hence the error at a given time t is given by Now, eq. (33) is for i = 1(1)N. Adopting the procedure followed in Theorem. 1, for i = 1(1)N, we have from (31), Now, for i = 1(1)N, we have the relation Following the procedure in Theorem. 1, Taking absolute value on both sides, Let k = max y (t m ) , for m = 0, 1, 2,...N-1, then, for i = 1(1)N, From (30) and (38), we have the required result, for all i = 0(1)N Hence the proof.
Numerical rate of order of convergence Let y(t) be the solution of the given differential eq. (1) and y h i and y h 2 2i are the numerical solution of the numerical method using step sizes h and h 2 , respectively. Then, we have an error estimate of the form, Here, n is the numerical rate of order of convergence of the method and it can be obtained from the estimate (40 (40) and (41), the numerical rate of order of convergence of the numerical method (3) is of two (n=2).

Experimental results
In this section, both numerical and graphical results for four stiff problems are given. Test problems considered are Test problem 1. y = −λy, y(0) = 1. whose solution is y(t) = exp(−λt).
Using Dahlquist problem (Test Problem 1, λ = 1), experimental results are provided concerning theoretical and numerical rate of order of convergence, absolute and relative errors, and their percentage. Using test problems 1 and 3, comparison is made with other order one and order two numerical methods. Graphical results are given using the test problems 1, 2, and 3.

Theoretical rate of order of convergence
If the exact solution is known then the theoretical rate of order of convergence of a numerical method P N is defined as . N refers to the number of nodal points for a particular step size h. The average rate can be defined as P = 1 8 P N . It is observed from Table 1 for the Dahlquist problem (Test Problem 1, λ = 1), the theoretical rate of the order of convergence is two.

Numerical rate of order of convergence
If the exact solution is not known then the numerical rate of order of convergence of a numerical method P N is defined as . N refers to the number of nodal points for a particular step size h. The average rate can be defined as P = 1 7 P N It is observed from Table 2 for the Dahlquist problem (Test Problem 1, λ = 1), the numerical rate of the order of convergence is two.
From Tables 1 and 2, it is observed that for the Dahlquist problem (Test Problem 1, λ = 1), we get the rate of the order of convergence as two.

Order graph plot
The order graph for the solutions of the Dahlquist problem (Test Problem 1, λ = 1) and the numerical method (3) is plotted in Fig. 9.
It is a plot of h and Q(h). From the exact solution of Dahlquist problem (Test Problem 1), Q (h)=exp (-h) and from the numerical method (3), we have Q(h) = Table 3 Exact and numerical solutions, absolute and relative errors-test problem 1 in [0, 1], λ = 1

Experimental results concerning the error
From the derivation of Theorem 1, it is observed that the local and global truncation errors get increased from time steps t = t 0 = a to t = t 1 and up to t = t N = b since the numerical method (3) is a recurrence relation between two consecutive time steps t = t i and t = t i+1 the error in time step t = t i get added with the error in the step t = t i+1 .
Similarly. absolute and relative errors get increased from time steps t = t 0 = a to t = t 1 and up to t = t N = b. This is illustrated in Table 3 using Dahlquist problem (Test Problem 1, λ = 1). It is observed from Table 3; the value of the absolute error and cumulative absolute error increases from the nodal point i = 0 to i = 16 and reaches maximum values 2.283E-04 and 2.7383E-03, respectively. Similarly, the value of the relative error and cumulative relative error increases from the nodal point i = 0 to i = 16 and reaches maximum values 6.2148E-04 and 5.2821E-03, respectively. Accumulation of error at each nodal point causes this effect.

Experimental results concerning percentage
The ratio of numerical and exact solution s i = y i y(t i ) and its percentage is calculated. The expected ratio value is 1 at all nodal points and its percentage is 100 at all nodal points. And, the percentage of relative error is zero at all nodal points. The sum of the ratio percentage and the relative error percentage at each nodal point is 100.
It is observed from Table 4 using Dahlquist problem (Test Problem 1, λ = 1); at all nodal points, the ratio is 1 from i = 1 to i = 16. The percentage of the ratio is 100 at all nodal points from i = 1 to i = 16. The maximum percentage of relative error is 0.0062148.

Comparative results concerning absolute error
Using Dahlquist problem (Test Problem 1, λ = 16 in [0, 10]), Table 5 presents errors due to all schemes concerning maximum absolute errors. The numerical method (3) proposed in this article is compared with explicit Euler Forward method, implicit Euler backward method, the explicit Runge-Kutta methods (Improved Euler method, the midpoint rule and Heuns method) of order two, and implicit Trapezoidal method of order two. In Table 5, R-K methods refers to explicit Runge-Kutta methods. From Table 5, it is observed that explicit Runge-Kutta methods, and implicit Trapezoidal method are of order two. Euler's forward and backward methods are of order one and the new method (3) is of order two. It is observed from the Table 5, the method (3) is better than other order one and two methods for large values of the step size h.
In Table 6, using Test Problem 3 in [0, 25], the numerical method (3) is compared with other methods. Explicit Euler method of order one and explicit Runge-Kutta methods of order two are not suitable for large step size for the Problem 3, and so implicit Euler method and Trapezoidal methods can be compared with method (3). It is observed from Table 6,   10, 11, and 12, respectively, using the step size h = 1 16 . The numerical solution is compared with the exact solution and the absolute error is also plotted. It is observed that the numerical solution converges with the exact solution and reflects the asymptotic behavior of the solution of the stiff problem.
Finally, the numerical method (3) will overcome the effects of stiffness and will reflect the asymptotic behavior  (3) is both A-stable and L-stable.

Conclusions
This article presents a finite difference method of order two for stiff problems that will overcome the effects of stiffness since it is A-stable. And, reflect the asymptotic behavior of the solution of the stiff problem since it is L-stable. The new method works well for large values of step size and works well in large domain. The classical explicit Runge-Kutta methods of order two are not suitable for stiff problems since the numerical solution will not overcome the effects of stiffness and will not reflect the asymptotic behavior of the solution of the stiff problem (for example, Dahlquist problem). And, they are not both Astable and L-stable. On applying Euler's implicit Backward method, the numerical solution will overcome the effects of stiffness and will reflect the asymptotic behavior of the solution of the stiff problem. And hence, it is both A-stable and L-stable. But it is of order one method.
A new numerical method which is both A-stable and L-stable and of order two for the numerical solution of a stiff differential equation is presented in this article. It will overcome the effects of stiffness and reflects the asymptotic behavior of the solution. The mew method is a modified form of the mid-point rule for the stiff problems. The rate of order of convergence is proved to be two both theoretically and numerically. Experimental results show the performance of the method based on the metrics such as stability function, stability region, order star fingers, theoretical and numerical rate of the order of convergence, absolute and relative errors, percentage of numerical solution accuracy and local and global truncation errors both numerically and graphically.
Calculating the shape of the unknown curve which starts at a given point and satisfies the given stiff differential equation is our problem. The shape of the unknown curve is got by the new A-stable and L-stable numerical method with order two convergence.
In brief, for stiff problems, it is proved that the numerical method (3) is of order two. This means, 1. The error per step is proportional to the cube of the step size, 2. The error at a given time is proportional to the square of the step size.
Both theoretical and numerical rate of the order of convergence is proved to be two. The new method will overcome the effects of stiffness since it is A-stable. And, also, the new method has stronger property than the A-stability. Stability regions, order star, order star finger region, relative order star region, relative absolute region and order graph are presented for the stiff problem. We hope that Scientists and Engineers will get in deep of the development of numerical methods in the path of development of stability region in the complex plane, percentage of numerical solution accuracy and theoretical and numerical rate and the average rate of the order of convergence. We suggest to use the order two method presented in this article instead of the classical explicit Runge-Kutta methods of order two since it is unconditionally stable to the Dahlquist problem (4) on the full right half of the complex plane. And, it will overcome the effects of stiffness and reflects the asymptotic behavior of the solution. The future research work which will be this method will be applied to system of stiff initial value problems and also for problems of higher dimension.