Stability Analysis of Singly Diagonally Implicit Block Backward Differentiation Formulas for Stiff Ordinary Differential Equations

: In this research, a singly diagonally implicit block backward differentiation formulas (SDIBBDF) for solving stiff ordinary differential equations (ODEs) is proposed. The formula reduced a fully implicit method to lower triangular matrix with equal diagonal elements which will results in only one evaluation of the Jacobian and one LU decomposition for each time step. For the SDIBBDF method to have practical signiﬁcance in solving stiff problems, its stability region must at least cover almost the whole of the negative half plane. Step size restriction of the proposed method have to be considered in order to ensure stability of the method in computing numerical results. Efﬁciency of the SDIBBDF method in solving stiff ODEs is justiﬁed when it managed to outperform the existing methods for both accuracy and computational time


Introduction
Many problems in engineering, physical and social sciences are reduced to quantifiable form through mathematical modelling involving systems of ordinary differential equations (ODEs).These problems sometimes exhibit a phenomenon known as stiffness.It associates with components that are decaying at widely differing rates.For this article, our main concern is the linear system of first order stiff ODEs of the form where y T = (y 1 (x), y 2 (x), . . ., y d (x)), f T = ( f 1 (x), f 2 (x), . . ., f d (x)) and µ T = (µ 1 (x), µ 2 (x), . . ., µ d (x)).
Equation ( 1) is said to be linear if f (x, y) = A(x)y + Φ(x), where A(x) is a constant d × d matrix and Φ(x) is an d-dimensional vector.
There are various definitions of stiffness that exist in the literature.However, we consider the one given by [1] which stated that the linear system (1) is said to be stiff if 1. Re(λ i ) < 0, i = 1, 2, . . ., d and Most realistic stiff systems do not have analytical solutions, hence a numerical approach have to be used [2].An exact solution of (1) can be approximated by using either the one-step method or the linear multistep method (LMM).Euler's is the simplest form of one-step method whereas the Runge-Kutta (RK) method is the most famous family under this class.These methods use the solution of current point, y n , as initial value to compute solution at the next point, y n+1 .

max
On the contrary, LMM uses information from the previous points to calculate the next values.Adams and BDF method are examples of families under the LMM.A classical approach for numerical methods in finding an approximation of y n+1 at the point x n+1 is computed one step at a time.
In order to accelerate the computational process, block method is introduced with the idea of producing simultaneously k-blocks where each block contains r-point approximation, y n+1 , y n+2 , . . ., y n+r , at each iteration of the algorithm.For the multistep block method, points from the previous block calculated are used in generating approximation solutions for the new block.The block method is first proposed by [3] which is later extended by several researchers such as [4][5][6].A modified block by Adams method for higher order ODEs had been discussed by [7] followed by [8] for solving nonstiff higher order ODEs.
As for stiff ODEs, fixed coefficients BBDF introduced by [9] had proved to optimize performance in both accuracy and computational time when solving stiff problems.Not only that, but the method also managed to outperform the non-block variable step variable order BDF (VSVOBDF) method by [10] (as cited in [9]).Numbers of improvement had been made to the original BBDF method, one of them is the diagonally implicit 2-point BBDF (DI2BBDF) method by [11] that gives better accuracy for high dimension problem of ODEs than the BBDF method.The method is derived based on the motivation by [12] which defined the diagonally implicit as the method with its coefficients of the upper-diagonal (read: lower triangular) entries are zero.
Over the years, conventional RK methods evolutionized to form the best solver for ODEs.The commonly improved RK methods involved an implementation of singly diagonally implicit properties which had proved to produce better computational time when compared with the existing methods.Among the initial works on the singly diagonally implicit RK (SDIRK) method is [13], which defined a singly diagonally implicit as a method with equal diagonal elements, α ii = γ, as shown in the Butcher tableau below.
The proposed problem of stiff ODEs is also considered by [14] through Newton-type iterations that solve for linear systems at each stage with coefficient matrix of the form I − hα ii ∂ f ∂y .Based on [15], singly implicit RK is a transformed method whose RK matrix has just one real s-fold eigenvalue.If we have α ii = γ, then the class of diagonally implicit RK (DIRK) is called SDIRK methods [16].
The definition of SDIRK by Norsett is further discussed by [17][18][19] in their research.As summarized by [20], when having the lower triangular matrix with equal diagonal elements, the stored LU-factorization of a single such matrix can be used repeatedly.Hence, only one evaluation of the Jacobian and one LU decomposition will be needed for each time step.By having these properties, degree of implicitness can also be reduced as it involved less computational process which will results in less execution time [21].
Therefore, developing the SDIBBDF method involves the hybrid process of implementing qualities from SDIRK method to block multistep method.As we know, both methods are of different families hence, our main concern is the compatibility of the derived method to solve for stiff ODEs.Derivation of the SDIBBDF method is shown in the next section.

Derivation of 2-Point SDIBBDF Method
Our main objective is to develop a method that capable of computing two solutions, y n+1 and y n+2 , simultaneously in less expensive environment with accurate approximation to the exact solution of stiff ODEs.The idea is illustrated as shown in Figure 1.
We consider the linear difference operator for block multistep method (2) as Equation ( 3) is expanded by using the Taylor series, and the terms of derivative y are collected to produce L s (y(x The constant C q is defined as Since the SDIBBDF method proposed is of order 2, k = 3 and q = 2 are substituted into Equation (4) to get By using MAPLE, Equation ( 6) is solved simultaneously, and the coefficients obtained are substituted into Equation (2).The formula derived is rewritten in matrices form as shown below: Therefore, the general formula of 2-point SDIBBDF method of order 2 is as follows: In the next subsection, we will discuss on stability of the derived SDIBBDF method.

Stability Analysis
One of the practical characteristics for a method to be useful is that it must have a region of absolute stability.1) is said to be absolutely stable in a region R for a given H if and only if for that H, all the roots, r s = r s (H) of the stability polynomial of the linear k-step method, π(r, H) = ρ(r) − Hφ(r), satisfy |r s | < 1, s = 1, 2, . . ., k where H = hλ and ρ(r) and φ(r) are the first and second characteristic polynomials respectively.Otherwise the method is said to be absolutely unstable.

Definition 1. The LMM in Equation (
In order to analyse the stability properties of the proposed method, a stability graph of the SDIBBDF method has to be constructed.First, characteristic polynomial of the method is determined by referring to the following statement.= λy (9) to Equation ( 7), we get where A and B are properly chosen r × r matrix coefficients and m represents the block number.
We consider ρ(t) = det(At − B) as the first characteristic polynomial of the SDIBBDF method.By substituting Equation (10) Based on Definition 2, when we solve for ρ(t) = 0, stability polynomial, R(H), of Equation ( 8) is obtained.
In order to determine zero stability of the SDIBBDF method, we refer to the following definition.
Definition 3. The LMM in ( 1) is said to be zero-stable if no root of the first characteristic polynomial, ρ(ζ), has modulus greater than one, and if every root with modulus one is simple.
Solving R(H) = 0 yields the following roots: Therefore, Equation ( 13) proved that the SDIBBDF method is zero stable.
To plot R(H), we let for each value of H, R is a complex number.Boundary of the stability region is the set of all H such that R(H) is on the unit circle of for some θ ∈ [0, 2π].Equation ( 14) is expanded for various value of θ in steps of 2π n from 0 to 2π.Next, the points are plotted by using MAPLE as shown in Figure 2 while Figure 3 represents the output which is the stability graph of 2-point SDIBBDF method.By referring to [22] a method is suitable for solving stiff problems because of its A-stability property as reviewed in the following statement.Definition 4. A numerical method is said to be A-stable if its region of absolute stability contains the whole left-hand half-plane, Re(hλ) < 0.
As observed in Figure 3, the area inside the closed region is unstable whereas the stable part lies outside the region.Coincide with Definition 4, we can conclude that the 2-point SDIBBDF method of order 2 is A-stable.
Suppose that the following definition is applied.
Definition 5.The LMM in Equation ( 1) is said to have region of absolute stability, A , where A is a region of the complex ĥ-plane, if it is absolutely stable for all ĥ ∈ A .The intersection with real axis is called the interval of absolute stability.
For SDIBDDF method, its interval of absolute stability is 4 < H < 0. Comparison in terms of stability region is made between the proposed method, fully implicit BBDF by [9] and DIBBDF method by [11] as shown in Figure 4. Analysis on Figure 4 concludes that the unstable region of the proposed method is smaller when compared with the other methods.Interval of the unstable area for each method is presented in Table 1.From the table shown, it is noted that the unstable area of DIBBDF method is wider followed by the BBDF and SDIBBDF method.This proved that the SDIBBDF method has wider stability area than the comparing methods.

Step Size Restriction
As stated in Definition 1, an LMM is said to be absolutely stable when |R(H)| < 1, otherwise it is unstable.There are two parameters involved, h and λ, but it is only their product, H, that needs to be taken into accounts.
By solving |R(H)| < 1, we found that the SDIBBDF method is stable everywhere except when H ∈ [0, 4].However, h must lies within a certain range in order for the SDIBBDF method to be stable.By substituting endpoint of the interval into characteristic polynomial in Equation ( 12), we obtain Next, Equation ( 15) is presented in the form of |ε| < 1 which equal to Solving Equation ( 16) yields |H| < 0.6240609349.
Please note that H = hλ, h < 0.6240609349 λ . ( Thus, Equation ( 17) is the step size restriction of SDIBBDF method.By taking an example of a stiff ODEs with eigenvalue, λ = −100, h < 0.6240609349 −100 , the suitable step size for solving the problem must be h < 0.006240609349.
If we take h > 0.006240609349, the method will be unstable therefore, numerical solutions with large maximum error is produced.

Numerical Results
In the previous section, we claimed that the proposed method is A-stable hence, suitable for solving stiff problems.To validate this finding, we tested the SDIBBDF method to solve for single and system of stiff ODEs.Efficiency of the method is measured in terms of its accuracy and execution time.
Accuracy is analysed based on the maximum error, where T is the total number of steps and N is the number of equations.By considering the error as with y i and y(x i ) are the approximate and exact solutions of Equation ( 8) respectively.
The following notations are used SDIBBDF : Singly diagonally implicit BBDF method DIBBDF : Diagonally implicit BBDF method by [23] BBDF : BBDF method by [24] ode15s : VSVO solver based on the numerical differentiation formulas (NDFs) ode23s : modified Rosenbrock formula of order 2 h : step size MAXE : maximum error TIME : computational time To verify performance of the methods, the following test problems of stiff ODEs are solved.
Test Problem 5 Exact solution: Tables 2-6 represent the numerical approximations of SDIBBDF, DIBBDF, BBDF, ode15s and ode23s for stiff ODEs with various step size.The results are illustrated in the form of efficiency curves as shown in Figures 5-9.Figures 5a-9a indicate performance of the methods based on step sizes tested.Meanwhile, Figures 5b-9b give visual representations on efficiency of the methods in terms of computational time to approximate solutions.Matlab solvers, ode15s and ode23s, are chosen as comparing methods since both solvers are designed for solving stiff problems.Ode15s is a variable order solver based on the numerical differentiation formulas.For certain occasions, it uses the BDFs method.On the other hand, ode23s is a one-step solver based on a modified Rosenbrock formula of order 2. Since results for SDIBBDF, DIBBDF and FIBBDF methods are computed by using C++ programming languages, the TIME of those methods cannot be compared with the one from Matlab.
From the results, one can observed that Matlab solvers, ode15s and ode23s, managed to obtain smaller MAXE than the SDIBBDF method when computing solutions for h = 10 −2 for every test problems.Numerical results in Table 3 show that SDIBBDF method produced slightly bigger MAXE than the fully implicit method when h = 10 −2 .This is due to the step size restriction as discussed earlier in Section 2.3.Since test problem 2 has λ = −100 hence, it requires h < 0.006240609349.
Not only that, but step size restriction is also experienced when h = 10 −2 is used to solve test problem 3 with λ = −100 and −1.Notice that the SDIBBDF, DIBBDF and BBDF method produced bigger MAXE for test problem 4 and 5 than the other test problems for stepsize h = 0.01.The errors obtained by both test problems are affected by the nature of the system where test problem 4 involved a highly stiff ODEs with λ = −1000 and −1 while test problem 5 consists of complex eigenvalues, λ = −40 − 40i, −40 + 40i and −2.
On the contrary, efficiency of the SDIBBDF method is proven when smaller stepsizes are used.It is capable of outperforming the other comparing methods for h = 10 −4 , h = 10 −6 and h = 10 −8 for each problems of stiff ODEs tested.
As for TIME, the proposed method is capable of executing solutions faster than the fully and diagonally implicit methods for each test problem.This is motivated by the implementation of singly diagonally implicit properties that had proved to improve efficiency of the method.Figures 5b-9b show that smaller step sizes require longer computational time due to larger number of steps involved.

Conclusions
The 2-point SDIBBDF method of order 2 for solving stiff ODEs is successfully derived by implementing the SDIRK properties of one-step method to the multistep block method.Stability analysis of the proposed method shows that it is zero-stable and with the A-stable characteristics which makes it fit for solving stiff problems.
Based on Section 2.3, we analysed the relation between step size and eigenvalue of stiff ODEs which gives effect on the MAXE produced.Hence, we will be able to set the best step size to used for solving stiff ODEs based on its eigenvalues in order to ensure good performance of SDIBBDF method to approximate solutions of the problems.
From numerical results, it is proved that the SDIBDDF method produced better accuracy for smaller stepsizes than the other comparing methods, and capable to execute solutions faster than the DIBBDF and fully implicit BBDF method.Effects of the step size restriction are also shown on the numerical results and analysed further.
Therefore, we can conclude that the SDIBBDF method can be used as an alternative solver for solving stiff ODEs.

i
|Re(λ i )| >> min i |Re(λ i )|, where λ i are the eigenvalues of A and the ratio S = max i |Re(λ i )| min i |Re(λ i )| is called the stiffness ratio.

Figure 1 .
Figure 1.2-point SDIBBDF method of constant step size.Points x n−1 and x n are the back values used to evaluate solutions of future points, x n+1 and x n+2 , with constant step size.Hence, we proposed the 2-point SDIBBDF method of order 2 with k+s−2 ∑ j=0

Figure 2 .
Figure 2. MAPLE codes to construct stability graph of 2-point SDIBBDF method.

Table 1 .
Interval of unstable area on real and imaginary axis.

Table 2 .
Numerical results for test problem 1.

Table 3 .
Numerical results for test problem 2.

Table 4 .
Numerical results for test problem 3.

Table 5 .
Numerical results for test problem 4.

Table 6 .
Numerical results for test problem 5.
(a) Graph of log h vs. log MAXE (b) Graph of log TIME vs. log MAXE