3-Point block backward differentiation formula with an off-step point for the solutions of stiff chemical reaction problems

A major challenge in simulating chemical reaction processes is integrating the stiff systems of Ordinary Differential Equations (ODEs) describing the chemical reactions due to stiffness. Thus, it would be of interest to search systematically for stiff solvers that are close to optimal for such problems. This paper presents an implicit 3-Point Block Backward Differentiation Formula with one off-step point (3POBBDF) for the solutions of first-order stiff chemical reaction problems. In deriving the method, the Lagrange polynomial was adopted as the basis function. The paper further analyses the basic properties of the 3POBBDF which include order of accuracy, consistence, zero-stability, and convergence. The stability region as well as the interval of instability of the method was also computed. To demonstrate the accuracy of the proposed approach, some famous stiff chemical reaction problems such as Robertson problem and Chemical AKZO were solved, and the results obtained were compared with those of some existing methods. The results obtained clearly show that the 3POBBDF performs better than the existing methods with which we compared our results.


Introduction
Most chemical reaction problems have been known to be stiff in nature, little wonder they have attracted the attention of mathematicians recently. According to [1,2], chemical reactions in alive species, atmospheric phenomena, mechanics and molecular dynamics, and chemical kinetics, extensively describe the stiff systems. If the time constants fluctuate significantly in magnitude, then the system is known as stiff. The time constants here refer to the rates of decay of perturbations [3]. If a system diverges considerably in terms of time constants, we must use smaller time steps to deal with expeditiously decaying terms.
In this paper, an implicit 3-Point Block Backward Differentiation Formula with one off-step point (3POBBDF) is derived for the solutions of stiff chemical reaction problems of the form, where f is assumed to be continuous and satisfies the existence and uniqueness theorem.
For Ordinary Differential Equations (ODEs) such as non-linear or linear, stiff or non-stiff ODEs, different numerical methods have been developed [4][5][6]. It must be noted that employing the wrong method for a model can result in a slow or incorrect solution [7]. Problems with the form (1.1) are generally divided into two categories. The first category is a non-stiff ODEs, in which explicit methods are employed together with certain error control. Another category is stiff ODEs. [8] was the first to use the term "stiff" in a scenario requiring chemical kinetics. The only way to find answers to stiff problems is to employ implicit methods because explicit methods are slow and sometimes fail to deliver an accurate solution [9]. The Backward Differentiation Formula (BDF) is one of the well-known classes of implicit methods for finding the solution of stiff ODEs. Developing a mathematical model for stiff ODEs considers a few factors such as selection of step size, stability, accuracy, and computational cost [10,11].
Stiffness is defined in different ways in literature as it does not have a specific definition. Equation (1.1) is said to be stiff according to Lambert [12] if the Jacobian matrix f y for the eigenvalues t (x) satisfies these conditions: where t are the eigenvalues and the ratio max t=1,2,…,m| Re( t )| min t=1,2,…,m| Re( t )| is called stiffness ratio.
Many studies have been directed into the improvement of multistep block methods for the solution of (1.1). Several studies have recently focused on improving the performance of different forms of multistep block methods encompassing estimated solution accuracy and computational time; see [10,[13][14][15][16][17][18][19]. While several block methods (based on Adams formulas) were proposed previously emphasizing proposed method is that the solutions are approximated with off-step points concurrently and provide better accuracy. The organization of the paper is as follow; Section 2 briefly describes the formulation of the method. Analysis of convergence properties is examined in Section 3. Section 4 elaborates the implementation of the derived method. List of tested problems with results are presented in section 5. Section 6 includes discussions of the results. Graphical representation of the results is displayed in section 7. consists of Numerical problems with their results respectively. Furthermore, Sections 8 describes the conclusions.

Derivation of the 3POBBDF
This section consists of a detailed explanation of how the proposed 3POBBDF for solving (1.1) is formulated. The two starting values, x n−1 and x n of identical step size which is fixed as h = x n+1 − x n and one off-step point at x n+ 5 2 is considered, as shown in Fig. 1.
The interval [a, b] in (1.1) is split into a sequence of blocks in a 3-point block method (see Fig. 1). For finding out the solution of (1.1), y n−1 and y n are used as the values of previous block to evaluate the future points at y n+1, y n+2 , y n+ 5 2 and y n+3 with constant step size simultaneously. The off-step point y n+ 5 2 is added to all three formulas to increase the accuracy. The 3POBBDF we shall derive is a k-step Linear Multistep Method (LMM) defined according to [12] as, h represents the step size, j and j are the unknown constants that needs to be decided and k indicates the number of steps of the method. Presuming that j ≠ 0 , and both 0 and 0 are not equal to zero. The LMM (2.1) is extended by adding an off-step point to generate a new form of LMM given by, The function f n+T in Eq. (2.2) is evaluated by using Lagrange interpolation polynomial (P(x)) with interpolation points such as ( x n−1 , y n−1 ), ( x n , y n ), ( x n+1 , y n+1 ), ( x n+2 , y n+2 ), ( x n+ 5 2 , y n+ 5 2 ) and ( x n+3 , y n+3 ). The process is given as, whereas, Expanding the polynomial L j by substituting s = x−x n+3 h implying x = sh + x n+3 . The interpolating polynomial is further differentiated with respect to s provides the coefficient values of and as presented in Table 1.
Furthermore, (2.2) is shaped up after substituting all values taken from Table 1. Hence, the corresponding formulas for the points y n+1 , y n+2 , y n+5∕2 , and y n+3 take the following form: Equation (2.4) is the 3POBBDF for the solution of stiff chemical reaction problems of the form (1.1).

Analysis of basic properties of the 3POBBDF
The following section comprises the analysis of basic properties of the proposed method. These properties include order, consistency, zero-stability and convergence. The region of absolute stability of the method shall also be determined.

Definition 3.1 (Order and Error Constant)
The LMM (2.1) associated with the linear difference operator is said to be of order p if C 1 = C 2 = … = C p = 0 and C p+1 ≠ 0 . The term C p+1 ≠ 0 is called the error constant of the method According to [12] and [50], the local truncation error related with (2.1), for determining the order and error constant of the proposed method is well-defined by the linear difference operator L(y(x)) as, The function y(x + Th) with its derivative y � (x + Th) are expanded and collecting terms in (3.1) will be obtained such as, On the application of Eq. (3.2) on 3POBBDF derived in Eq. (2.4) using the coefficients presented in Table 1, we obtain As a result, the error constant of the proposed method is obtained as 3) According to Definition 3.1, c 0 = c 1 = … = c 5 = [0000] T . The term c 6 represents the error constant. Therefore, this can be demonstrated that the proposed method is of order five.

Consistency of the method Definition 3.2 (Consistency)
The LMM (2.1) is said to be consistent if it is of order p ≥ 1.
It is therefore clear that the 3POBBDF is consistent since it has order p = 5 which satisfies Definition 3.2

Zero-stability of the method Definition 3.3 (Zero-stable)
If no root of the characteristic polynomial has a modulus greater than one and every root with modulus one is simple, then LMM (2.2) is said to be zero-stable To be a convergent method, zero-stability is one of the main criteria which will be discussed here. [13] proposed the scalar test to determine the stability of the method (2.4) as where λ represents the complex constant with Re(λ) < 0. Equation (3.4) is substituted in Eq. (2.4), therefore it precedes as, The matrix is then created by inscribing Eq. (3.5) into the matrix form as follows:  .7), the zero-stability is obtained by replacing H = h = 0 into (3.7), this yield, Solving R(t, 0) = 0 , gives the roots for t , as t = 0, 0, −0.000186375, 0.999767 . Since all the roots lie within |t| ≤ 1 described in Definition 4.1; hence, it is concluded that the method 3POBBDF is zero-stable.

Convergence of the method Definition 3.4 (Convergence)
Consistency and zero-stability are required conditions for the LMM (2.1) to be convergent Therefore, we conclude that the newly derived 3POBBDF is convergent since it is consistent and zero-stable.

Stability region of the method Definition 3.5 (A-stable)
The term "A-stable" refers to a method that is stable for all λ in the left-half plane (Re( ) ≤ 0) . The stability region for A-stable method covers up the complete negative left-half plane In this section, the stability region of the 3POBBDF(5) presented in Eq. (2.4) is plotted by using the stability polynomial given in (3.7). The boundary of the stability region is described by the set of points defined by t = e i ,0 ≤ ≤ 2 . The root condition (|t| ≤ 1) of the stability polynomial must be tested at multiple grid points in order to determine the boundary of the stability region. Figure 2 shows the region of absolute stability for the system in (2.4). It is observed from this figure that the stability region lies outside of the bounded region. Since most of the region is in the left half-plane therefore, according to Definition 3.5, the method is A-stable. Thus, it can be implied that the constructed 3POBBDF method is appropriate for solving stiff problems.

Comparison of stability regions
The stability region of the proposed 3POBBDF (5) is depicted in the Fig. 3 below, as opposed to the existing stability region of the Block Backward Differentiation Formula of Order 6 (BBDFO (6)) developed by [42]. On the complex h plane, Fig. 3 illustrates the region including both methods, namely the 3POBBDF(5) and the BBDFO (6).
The interval of instability of the current 3POBBDF(5) and the BBDFO(6) [42] are then compared. The derived 3POBBDF(5) has a larger stability area, as shown in Table 2. This implies that the proposed 3POBBDF(5) is more computationally reliable than BBDFO(6) developed by [42].

Implementation and Algorithm of the 3POBBDF
This section discusses the implementation and algorithms of the proposed method.

Implementation of the 3POBBDF
By using Newton iteration, implementation of the 3POBBDF(5) method having one off-step is briefly mentioned in the generalized form. An increment notation is introduced for stipulating the iteration. In which y (i+1) n+1 will represent the (i + 1) th iterative value of y n+1 [51] .
The main limitation for the implementation of the 3POBBDF(5) predominantly is that it is not self-starting. The back values are therefore calculated using the Euler method in this case.

Algorithm of the 3POBBDF in Code
This section will elaborate on how the developed code works with the proposed method. Initially, the code will begin with three points having one off-step block method. It is employed in PE(CE) r , the type where P and C stand for predictor and corrector correspondingly. Whereas E indicates the evaluation of function f .The power r = 2 shows the number of iterations as presented by [51] that is required to converge three points having one off-step block method corrector formulae. The algorithm of the proposed method in the code is given below, Step 1: P = Predictor values y n+m for m = 1, 2, 5 2 , 3, Step 2: E = Evaluate f n+m for m = 1, 2, 5 2 , 3, Step 3: C = Corrector value y n+m for m = 1, 2, 5 2 , 3, Step 4: E = Evaluate f n+m for m = 1, 2, 5 2 , 3.

Test problems
To test the reliability of the proposed 3POBBDF(5) method, some numerical results are obtained by applying the method on some well-known chemical reaction problems. Comparison is done with other multistep methods. The 3POBBDF(5) method has been coded in C + + .

Problem 5.1 (Stiff Chemical Reaction Problem)
A chemistry problem with a stiff system is considered below, with initial valuey 1 = 0,y 2 = 1 and This problem is solved using h = 10 −5 for 3POBBDF (5). Table 3 shows the integration results for this problem at x = 2 . Compared with the formulas of y � 1 = −0.013y 2 − 1000y 1 y 2 − 2500y 1 y 3 ,  [48], and Ismail's method [49], the 3POBBDF(5) method provides more accurate results in terms of Absolute error.

Problem 5.2 (Chemical AKZO Nobel Problem)
The chemical AKZO noble problem is a chemical process given by a stiff system with six non-linear differential equations taken from [52] is considered. Mathematically the problem is described as, The function F(y) is defined by where r i and F in are auxiliary variables given by The initial vector y 0 = (0.437, 0.00123, 0, 0, 0, 0.367) T . Background of the Chemical AKZO Nobel Problem: This problem comes from AKZO Nobel Central Research in Arnhem, the Netherlands. It defines a chemical process that mixes MBT and CHA with a continuous supply of oxygen. The important consequent species is CBS. AKZO Nobel problem contains following reaction equations,   Table 4.

Problem 5.3 (Robertson Problem)
A stiff system of three non-linear ODEs is involved in the Robertson problem. In 1966, Robertson proposed it [53]. [9] came up with the name ROBER. The file [CHA] � , [CHA], [CHA] 2 , rober.f contained the software part of the problem which can be found in [54]. The problem is described mathematically in the form.
with y ∈ ℝ 3 , x ∈ [0, X] and the function F is specified by with the initial value y 1 = 1,y 2 = 0 and y 3 = 0.
Origin of the Problem: Robertson [53] defined the ROBER problem as the kinetics of an autocatalytic reaction. The composition of the reactions is as follows: A , B, and C are the chemical species and k 1 , k 2 and k 3 are the rate constants. The following mathematical model, consisting of a set of three ODEs, can be put up under some ideal conditions [55] and the assumption that the mass action law is implemented to the rate functions.  with (y 1 (0), y 2 (0), y 3 (0)) T = (y 01 , y 02 , y 03 ) T , where y 1 , y 2 and y 3 refers to the concentrations of A, B, and C respectively, and y 01 , y 02 and y 03 are concentrations at time t = 0 . In the test problem, k 1 = 0.04 , k 2 = 3 × 10 7 and k 3 = 10 4 are the numerical values of rate constants with the initial concentrations y 01 = 1 , y 02 = 0 and y 03 = 0.
The cause for stiffness is the substantial variation in reaction rate constants. This system features a modest, very rapid initial transient phase for problems arising from chemical kinetics. This is followed by a highly smooth fluctuation of the components, which would be acceptable for a numerical technique with a large stepsize. For x ∈ [0,4000], the ODE system is integrated. The numerical results by proposed method at the points.

Problem 5.4 (Stiff Chemical Reaction Problem)
Consider the stiff system IVPs.
with the initial conditions as y 1 (0) = 1 and y 2 (0) = 1 having exact solutions as, This problem is solved at x = 50 by the new method and compared the results with forth order Two step hybrid method with one off-step point [48]. Here the stepsize h = 0.05 is used for compared and proposed methods. The smaller stepsize can also be used to get more accurate results. Table 6 shows the numerical results.

Discussion
The results generated by the developed 3POBBDF (5) method are displayed in Tables 3,  4, 5, 6. In Table 3, comparison on the basis of absolute error was made using the 3POBBDF(5) of order 5 at the point x = 2 with the methods of Khalsaraei and Molayi [48] and Ismail's method [49] at the stepsizeh = 10 −5 . In Table 4, 3POBBDF (5) method is integrated at h = 10 −5 and the direct comparison has been made with the two step hybrid method with one off-step point [48] by taking the solution points atx = 180 . Table 3 shows that the 3POBBDF (5) is more accurate and performs relatively better than the method [48] and significantly better than the method [49]. However, Table 4 shows the superiority of 3POBBDF(5) over the method [48]. Table 5 shows the numerical integration results for the problem 7.3 at x = 0.4, 40 and 4000 using the proposed 3POBBDF(5) approach in comparison to that of [48] at h = 0.001 . The comparison of the absolute errors in Table 5 for different points of x demonstrates the dominancy of 3POBBDF(5) over the method given in [48] in terms of accuracy ath = 0.05 . According to Table 6, it is observed that the 3POBBDF(5) method of order 5 is more effective than the two-step hybrid method of forth order with one off-step point [48]. In Table 6, it can be seen that the 3POBBDF(5) method of order 5 performs better than the fourth order Two step hybrid method with one off-step point [48]. As shown by tables (3-6), 3POBBDF (5) is well suited to nonlinear chemical reaction problems of a stiff nature.

Graphical representations of numerical results
The tested chemical reaction problems (5.1)-(5.4) are also compared with ode15s. Figures 4, 5, 6, and 7 are the graphs demonstrating the accuracy of the 3POBBDF(5) compared with the stiff solver ode15s. It can be examined from the graphs that the results of the 3POBBDF(5) method and the results of ode15s are nearly identical at specific intervals. These figures demonstrate that the proposed method is well suited for stiff chemical reaction problems, since the solutions produced by the 3POBBDF(5) coincide with the well-known ode15s code of MATLAB.

Conclusion
In the present paper, an implicit 3POBBDF(5) method has been derived for the numerical solution of stiff systems of first-order IVPs arising from chemical reactions such as Chemical AKZO, Robertson problem, and stiff chemical problems. The approach is based on the block BDF, in which at each step of the integration, three approximate solutions are generated instantaneously. Based on stability analysis of the 3POBBDF(5), the method is consistent and zero-stable; thus, the 3POBBF(5) is convergent. Since it has an A-stability property, the method is declared suitable for solving stiff ODEs. Comparison between the stability regions was also made which shows that the 3POBBDF(5) has a smaller interval of instability than the BBDFO(6) suggested by [42]. The numerical results obtained through the 3POBBDF(5) and compared to [48] and [49] evidenced that by applying the 3POBBDF(5) method, the accuracy of the numerical solutions in terms of absolute error at specific points is improved. Hence, the proposed method can be successfully applied on stiff systems generated from chemical reactions because of their high order accuracy and wider stability region. Therefore, it can be concluded that the 3POBBDF(5) could be an appropriate substitute solver for stiff ODEs. To improve efficiency, further research can be implemented by using variable step sizes for solving ODEs would be beneficial.