Formulation of a New Implicit Method for Group Implicit BBDF in Solving Related Stiff Ordinary Differential Equations

This paper proposed a new alternative approach of the implicit diagonal block backward differentiation formula (BBDF) to solve linear and nonlinear first-order stiff ordinary differential equations (ODEs). We generate the solver by manipulating the numbers of back values to achieve a higher-order possible using the interpolation procedure. The algorithm is developed and implemented in C ++ medium. The numerical integrator approximates few solution points concurrently with off-step points in a block scheme over a non-overlapping solution interval at a single iteration. The lower triangular matrix form of the implicit diagonal causes fewer differentiation coefficients and ultimately reduces the execution time during running codes. We choose two intermediate points as off-step points appropriately, which are proven to guarantee the method's zero stability. The off-step points help to increase the accuracy by optimizing the local truncation error. The proposed solver satisfied theoretical consistency and zero-stable requirements, leading to a convergent multistep method with third algebraic order. We used the well-known and standard linear and nonlinear stiff IVP problems used in literature for validation to measure the algorithm's accuracy and processor time efficiency. The performance metrics are validated by comparing them with a proven solver, and the output shows that the alternative method is better than the existing one.


Introduction
In applied science and engineering problems, the inclusion of ordinary differential equations (ODEs) in a dynamic model is vital in demonstrating the state of the system. The higher-order and complex phenomena of ODEs barely have analytical solutions, and it advocates using numerical methods to approximate the solutions for the problems. Frequent dynamic ODE models are governed by the unique behavior identified as stiffness. The stiffness caused the solutions to become unbearably slow, and only a few ODEs integration methods can be applied in all spatial locations. On top of that, a wide range of stiffness may exist in the standard practice set of ODEs. Thus, selecting low-cost and economical numerical methods is urging based on local conditions. Although stiff ODEs are the matured area in practices, the need to produce quality and improved methods in line with the systems' demand and difficulty cannot cease [1].
Our primary interest in this paper is solving an initial value problem (IVP) of ordinary differential equations (ODEs) of order 2 (two), which can be expressed:  [2].
According to [1], stiffness features are classified through the eigenvalues of the Jacobian matrix that exist in the left half-plane of the complex plane, alternatively stated that the i  must have a maximum negative magnitude. The second condition for stiffness is outlined by the absolute value of the stiffness ratio between the maximum and minimum of eigenvalues that must contain maximum negative values.
Stiff ODEs produced a part of solutions called the transient phase. The solution rapidly decays zero as t increases and requires the stability equation to take a smaller step size. It became impracticable for forwarding step methods [3].
As reported by [1], the implicit methods are significantly fitted in solving stiff ODEs because of their higher-order accuracy as compared to explicit ones. However, [4] set up the limitation for stability and proved that stable multistep-method did not exceed order over two when solving stiff ODEs. The alternative approaches in [5] and [6] had proposed to overcome the limitation. These approaches involved modifications of the multistep method procedure for achieving a higher-order and stable method. In 1964, [7] introduced a hybrid multistep method for the IVP to overcome the Dalquist's barrier in [4]. The hybrid terms classified as n h f    where 1 / 2 or 1/3.   in the method and provided stable and convergent solutions. Besides, [8] proposed a strategy to improve the accuracy and order of the methods by adding the off-step points and future points. Relevant research on the hybrid multistep methods can be found in [9]- [11].
In 1952, [12] had overcome stiffness in their kinetic studies and stated that backward differentiation formula (BDF) was the best method for solving stiff ODEs. [13] improvised the BDF method with the new approach of calculating a few solution points simultaneously in a block at each iteration of the algorithm. It was identified as the block backward differentiation formula (BBDF) and produced the solutions with minimal time cost and lessened the number of integration steps. Since then, the study on BBDF captivated the attention and grew in numbers as seen in [15]- [19]. In recent studies by [17][18], a new approach of BBDF method proved that diagonally implicit block has effectively solved stiff ODEs.
The alternative formula introduced in this paper is a combination of diagonally implicit block BDF and off-step points. We simply call the formula DI2OBBDF. We present the formulation of the coefficients for DI2OBBD in the next section.

Derivation of the Method
This section presented the derivation of the coefficients for DI2OBBDF by using Lagrange Interpolation Polynomial.
According to [17] and [18], the derivation for diagonal implicit methods by using (2)  x y (x n , y n ).
For second point, For the third point of the method The last point 2 n y  of formula is obtained by interpolates the points comprise The next step is to introduce the variable substitution (5) and (6) 146 Formulation of a New Implicit Method for Group Implicit BBDF in Solving Related Stiff Ordinary Differential Equations respectively. Then, we differentiate once for (3), (4), (5) and (6) with respect to s, following with substitution each formula with 0, 1 /2, 1, 3 /2 s  respectively, and the corrector formula is as follows: We stored the coefficients of DI2OBBDF in the algorithm to avoid the former tedious computation of the divided difference. We implemented the algorithm in PECE mode where P stands for predictor, E for evaluation of (1) and C for corrector in (7). The same technique in (2) is applied for the derivation of predictor formula. The interpolation points involve

Order of the Method
This section addressed the calculation of an order for DI2OBBDF in (7) by referring to several definitions in [1]. Beforehand, the DI2OBBDF can be written in the following matrix form: The corrector formula in (7) associated with linear operator can be written as follows: Taylor's series are used for expanding the function Substituting (13) into (12) gives The equations (14)  . Hence, the corrector method (7) is an order three method.

Convergent of the Method
The following theorem stated in [1] defines the convergent criteria for DI2OBBDF. (10) to be convergent are that it must be consistent and zero stable.

Theorem 1. The necessary and sufficient conditions for the linear multistep method (LMM) of
The terms r  and r  where 1, 2,..., 11. r  are used to classify the elements of the matrix in (7) as follows:   (7) is consistent.

Theorem 2.
The necessary and sufficient condition for the method given by (10) to be zero stable is that it satisfies the root condition. [1] We associated the DI2OBBDF with the scalar test equation, ( ) ( ) y x y x    with notation H h  for verification of root condition. We solved the formula by the determinant of the matrix and obtain the sets of roots that satisfied 1 i   Hence, DI2OBBDF satisfied the criteria of a convergent method.

Numerical Result
This section demonstrated the efficiency of the derived method by comparing it with a proven solver, implicit 2-point block BDF (BBDF) in [13]. We validate the derived method's performance via test problems that consist of linear and nonlinear stiff problems. The problems are listed as follows and are taken from [17]: We validate the efficiency of the proposed method with respect to maximum error in (16), and execution time with fixed step-size, h. We define the notations used in the table as follows: h : Step size. MAXE : Maximum error. TIME : CPU times in microseconds.

DI2OBBDF
: Diagonally Implicit Two-point BBDF with two off-step points.. BBDF : Implicit BBDF [13] 148 Formulation of a New Implicit Method for Group Implicit BBDF in Solving Related Stiff Ordinary Differential Equations Computation of maximum error is defined as Numerical result is as follows: The accuracy and execution time of the methods are presented in Tables 1-4. In conjunction, we plotted 10 Log ( ) MAXE log 10 (MAXE) versus 10 Log ( ) TIME graphs for each problem in Figures 1-4 to support the result obtained.        In Figure 1-4, we plotted the accuracy curves for the DI2OBBDF and BBDF. The figure depicts the scale of Log (MAXE) versus Log (TOL) for each solver, and the small magnitude ordinates indicate the minimum error and the faster the time for running code. Based on the result, DI2OBBDF converged fastest for all test problems than other proven solvers since the growth of error curves in figures was smaller in values and shorter in the distance.

Conclusion
We constructed an implicit diagonal BBDF with off-step points to accurately approximate linear and nonlinear stiff IVP solutions. The analysis through the order, consistency, and zero stability resulted in a convergent method. We found that the significant contribution of DI2OBBDF is the ability to produce multiple solutions in one iteration. Meanwhile, the implicit diagonal form of DI2OBBDF led to Jacobian Newton's coefficients matrix in a lower triangular form, thus shortening the time required to run the codes. The presence of off-step points optimized the local error by increasing the accuracy. The numerical outputs of four test problems in Section 5 demonstrate that DI2OBBDF is a powerful tool for solving linear and nonlinear ODEs. Therefore, DI2OBBDF is significantly useful for efficient approximations in a broad scope of ODEs in engineering and science.