On the Integration of Stiff ODEs Using Block Backward Differentiation Formulas of Order Six

In this research, a six-order, fully implicit Block Backward Differentiation Formula with two off-step points (BBDFO(6)), for the integration of first-order ordinary differential equations (ODEs) that exhibit stiffness, is proposed. The order, consistency and stability properties of the method are discussed, and the method is found to be zero stable and consistent. Hence, the method is convergent. The numerical comparisons with the existing methods of a similar type are given to demonstrate the accuracy of the derived method.

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 a d-dimensional vector assumed to be continuously differentiable, and if it satisfies the Lipschitz conditions as given in [1], which guarantees the existence and uniqueness of the solution of Equation (1). The development of numerical methods for the solution of Equation (1) is influenced by the types of ODEs, such as linear or non-linear ODEs, or stiff or non-stiff ODEs. Note that using the wrong type of method for a model can produce slow and/or inaccurate results. Generally, the problems in the form of Equation (1) can be classified into two types. The first type is non-stiff ODEs, for which all of the components evolve simultaneously and on comparable time-scales. Non-stiff problems are often solved by using explicit methods with some error control. The second type is stiff ODEs. The first appearance of the term "stiff" is in the paper by [2] on the problems in chemical kinetics. There are various definitions of stiffness given in the literature since there is no universally accepted definition. The problem in Equation (1) is said to be describing the stiff problem if: i.
It contains widely varying time scales, i.e., some components of the solution decay much more rapidly than others [3]; ii.
The step size is dictated by the stability requirements rather than the accuracy requirements [3]; iii.
Explicit methods do not work, or work only extremely slowly [4]; Symmetry 2020, 12 iv.
All of its eigenvalues have negative real parts, and the stiffness ratio (the ratio of the magnitudes of the real parts of the largest and smallest eigenvalues) is large [3]; v.
No solution component is unstable (no eigenvalue of the Jacobian matrix has a real part which is at all large and positive) and at least some components are very stable (at least one eigenvalue has a real part which is large and negative) [5].
For this article, we consider the definition of stiffness to be that given by [6] which stated that the linear system of Equation (1)  is called the stiffness ratio as a measure of stiffness.
The consideration of stability properties is particularly important when developing methods for solving stiff ODEs, and implicit methods which possess an A-stable region are preferred. Several numerical methods of various types have been proposed for the integration of both types of problems in Equation (1). Among the most popular and widely used methods for the solution of the stiff systems in Equation (1) are the Backward Differentiation Formulas (BDFs) proposed by [7]. In MATLAB, Gear's methods are implemented in the stiff solver ode15s.
For many years, the development of more advanced and efficient methods for the solution of ordinary differential equations (ODEs) has been an active research area. Most of the developments have been based on block methods. These block methods generate multiple solutions simultaneously and in single iterations. As noted by [8][9][10][11][12][13][14][15][16][17][18], these methods are effective for solving the problems in Equation (1). The block methods not only reduce the computational time but also require a smaller number of function evaluations per step. These block methods include the variable order Adams methods, the Runge-Kutta methods and Backward Differentiation Formulas. Inspired by Gear's method, [11][12][13][14][15][16][17][18] developed Block Backward Differentiation Formulas (BBDFs) which are proven to be efficient in solving stiff ODEs. For a detailed study on the theory and development of BBDFs, one can refer to [19]. Recently, the method that incorporated a BBDF proposed by [11] has been modified and extended, and it is now an improved method due to its enhancements in reducing the total number of integration steps and computational time, and in providing better accuracy of the approximation-see [12,20]. Therefore, in this paper, we formulate a new method called the Block Backward Differentiation Formula with two off-step points (BBDFO(6)) which is extended from the BBDF method developed by [14] by adding off-step points. These block methods with off-steps were studied by [21]. A closely related method with off-step points was described by [22][23][24]. A similar approach is adapted from [25,26] in adjusting the selection of the most suitable off-step points. We extend the theory to [11] and establish new coefficients of stability which is an improvement on the work of [21]. We shall restrict our discussion to solving stiff ODEs. This paper is organized as follows: in Section 2, we present the derivation of the coefficients of the new method; Section 3 includes the order, zero stability, and region of stability of the proposed method. In Section 4, numerical tests are performed on the first-order differential equation problems possessing stiffness and are compared to related works. Finally, the conclusions are presented in Section 5.

Derivation of BBDF with Off-Step Points, BBDFO
In this section, a detailed discussion of the derivation of the proposed method for the solution of Equation (1) will be presented. We considered the three points y n−2 , y n−1, y n of equal step size h given by h = x n+1 − x n as the starting values and the two off-step points as y n+1/2 and y n+3/2 with half the step size, as illustrated in Figure 1.
Hence, the BBDFO formulas are given by In the following sections, the order, zero stability, region of stability and convergence of the new BBDFO method will be discussed.

Order, Convergence and Stability of the Method
This section is concerned with the order, convergence and stability of the proposed method. The following definitions, as stated in [6], give the meanings of order, consistency, zero stable and A-stability by relating them to the definition of convergence.

Definition 3.1. (Convergence)
The necessary sufficient conditions for a linear multistep method to be convergent are that it must be consistent and zero-stable.

Definition 3.2. (Order)
Linear multistep methods associated with the linear difference operator are said to be order p if C 1 = C 2 = . . . = C p = 0 and C p+1 0.

Definition 3.3. (Consistent)
The linear multistep method Equation (2) is said to be consistent if it has order p ≥ 1.

Definition 3.4. (Zero stable)
The linear multistep method Equation (2) is said to be zero stable if no root of the first characteristics of the polynomial has a modulus greater than one, and if every root with modulus one is simple.

Definition 3.5. (A-stable)
If a method is stable for all λ in the left-half plane (Re(λ) ≤ 0), then the method is said to be A-stable. This means the stability region covers the entire negative left-half plane.

Order of the Method
The formula in Equation (12)  We expressed Equation (13) by the following matrix form By letting the collected terms of y(x), y (x), y (x) in Equation (16) as C p where p = 0, 1, 2 . . . , we can determine the order of the method and the error constant of the method.
CThe associated linear difference operator, the order and error constants for the method in Equation (12), can be defined by extending the corresponding definitions for a linear multistep method.
Expanding on Equation (14) with Taylor's series of expansions of y(x + jh/2) and y (x + jh/2) about x, the terms can be collected to obtain    We have This yields the error constant which shows that the order of the method is six. Next, we examine the consistency of the method. In accordance with Definition 3.3., the BBDFO(6) method is consistent since the order is greater than one.

Zero Stability of the Method
Next, we will discuss the zero stability of the method which is one of the criteria for the method to be convergent. The stability of the method in Equation (3) is determined using the linear test equation proposed by [27] as y = λy where λ is a complex constant with Re(λ) < 0. The substitution of Equation (17) y n−7/2 y n−3 y n−5/2 y n−2 By inputting H = hλ, the stability polynomial R(t, H) associated with the method in Equation (3) is determined by the location of the roots by solving the characteristic equation det At 2 − Bt − C = 0.
Since all the roots lie within |t| ≤ 1, we can conclude that the method for BBDFO(6) is zero stable. In the next subsection, we determine the absolute stability region for the BBDFO(6) method.

Stability Region of the Method
In this subsection, we plot the stability region of the method from the stability polynomial given in Equation (19). The boundary of the stability region is given by the set of points determined by t = e iθ , 0 ≤ θ ≤ 2π. It is necessary to test the root condition (|t| ≤ 1) of the stability polynomial at numerous grid points in the stability space to obtain the boundary of the stability region. The region of absolute stability for the method in Equation (12) is as shown in Figure 2.

Stability Region of the Method
In this subsection, we plot the stability region of the method from the stability polynomial given in Equation (19). The boundary of the stability region is given by the set of points determined by = , 0 ≤ ≤ 2 . It is necessary to test the root condition (| | ≤ 1) of the stability polynomial at numerous grid points in the stability space to obtain the boundary of the stability region. The region of absolute stability for the method in Equation (12) is as shown in Figure 2.  Figure 2 shows the stability region which corresponds to the BBDFO(6) method. The stable region lies outside the closed region and it covers the entire left-half plane. From Definition 3.5, we can conclude that the BBDFO(6) method is A-stable.  Figure 2 shows the stability region which corresponds to the BBDFO(6) method. The stable region lies outside the closed region and it covers the entire left-half plane. From Definition 3.5, we can conclude that the BBDFO(6) method is A-stable.

Stability Comparison
The following figure shows the stability region of the new BBDFO(6) as compared to the existing stability region of the order six Block Backward Differentiation Formula without the off-step points, denoted as BBDF(6). For further discussion of the BBDF(6) method, see [28]. Figure 3 shows the region of both methods, i.e., the BBDFO(6) and the BBDF(6) in the complex hλ plane. The following figure shows the stability region of the new BBDFO(6) as compared to the existing stability region of the order six Block Backward Differentiation Formula without the off-step points, denoted as BBDF(6). For further discussion of the BBDF(6) method, see [28]. Figure 3 shows the region of both methods, i.e., the BBDFO(6) and the BBDF(6) in the complex ℎ plane. Next, we compare the interval of instability of the new BBDFO(6) with the BBDF(6) by [28]. It is clear, from Table 1, that the derived BBDFO(6) has a larger stability region.

Numerical Results
In this section, the three problems that exhibit stiffness, given by , will be tested by using the proposed method with the step sizes ℎ = 10 −3 , 10 −4 , 10 −5 and 10 −6 . The efficiency of the method is compared to the exact solutions to illustrate the accuracy of the method. Problem 1. (See [29]).  Next, we compare the interval of instability of the new BBDFO(6) with the BBDF(6) by [28]. It is clear, from Table 1, that the derived BBDFO(6) has a larger stability region.

Numerical Results
In this section, the three problems that exhibit stiffness, given by λ, will be tested by using the proposed method with the step sizes h = 10 −3 , 10 −4 , 10 −5 and 10 −6 . The efficiency of the method is compared to the exact solutions to illustrate the accuracy of the method.  Exact solution: y(x) = e −1000x + 1.
Exact solution: For error calculations, the error formula is given by where y(x n ) is the exact solution and y n is the computed solution.
The MATLAB solver ode15 s was chosen to compare the methods since it is the most popular code for solving stiff problems. It is a variable order solver which is based on the numerical differentiation formulas (NDFs) and optionally comes with the backward differentiation formulas (BDFs), also known as Gears' method. The numerical results are presented in Tables 2-4.             As seen in Figures 4-6, the global errors produced by the BBDFO(6) are smaller global errors compared with the BBDF(6) and the ode15 s as the value of h becomes smaller.

Conclusions
In conclusion, a new method with off-step points, the BBDFO of order six, is introduced as the numerical solution of stiff initial-value problems. The stability analysis of the BBDFO indicates that the method is consistent and zero stable; therefore, the BBDFO(6) is convergent. The developed method is suitable to solve stiff ODEs since it possesses an A-stability property. The interval of instability for the BBDFO(6) is smaller than that of the BBDF(6) proposed by [28]. The numerical results obtained using the BBDFO(6) were compared with existing methods of a similar type and performed with the MATLAB solver ode15s. The BBDFO(6) is found to be more efficient at certain step sizes. Hence, we can conclude that the BBDFO(6) could be a compatible alternative solver for stiff ODEs.