A Class of Hybrid Multistep Block Methods with A — Stability for the Numerical Solution of Stiff Ordinary Differential Equations

: Recently, block backward differentiation formulas (BBDFs) are used successfully for solving stiff differential equations. In this article, a class of hybrid block backward differentiation formulas (HBBDFs) methods that possessed 𝐴 (cid:3398) stability are constructed by reformulating the BBDFs for the numerical solution of stiff ordinary differential equations (ODEs). The stability and convergence of the new method are investigated. The methods are found to be zero ‐ stable and consistent, hence the method is convergent. Comparisons between the proposed method with exact solutions and existing methods of similar type show that the new extension of the BBDFs improved the stability with acceptable degree of accuracy.


Introduction
The aim of this paper is to develop a class hybrid multistep block method based on block backward differentiation formulas (BBDFs) for the numerical solution of first order ordinary differential equations (ODEs) of the form: where is continuous and second variable Lipschitzian over its existence domain, so that the existence of Solution (1) is guaranteed. Equation (1) frequently arise in many applications in science and engineering. Many of these applications often results to stiff ODEs. Stiff systems are found in description of atmospheric phenomena, chemical reactions occurring in living species, chemical kinetics (e.g., explosions), engineering control systems, electronic circuits, lasers, mechanics, and molecular dynamics [1][2][3].
Several other examples are given in [4]. There is no unique definition of stiffness. Various definition describing stiffness are given in the literature [4][5][6]. The first identification of stiff equations as a special class of problems is due to a chemist [5]. He states that stiff equations are equations where certain implicit methods, in particular BDF, perform better, usually tremendously better, than explicit ones [6] and systems containing very fast components as well as very slow components while [4] states that "stiff problems are problems for which explicit methods don't' work".
In many practical applications, the difficulty associated with stiff ODEs is in choosing an efficient integration scheme since the choice of numerical methods is imposed by numerical stability. Besides numerical stability, the efficiency of numerical methods also depends on whether the method is explicit or implicit. Explicit schemes require the least amount of computation per time step but the allowable time step is severely restricted by the stability requirements, see [7]. Furthermore, the implementation of implicit method can be very costly in terms of computer time since solvers for stiff equations require the evaluation of the Jacobian, ⁄ .
The most widely used method to numerically solve (1) which is stiff are the backward differentiation formulas (BDFs) which is popularly known as the Gear method, see [8,9]. Recently, a great deal of effort has been proposed in the literature in the development of block methods for solving (1). The block methods are extensively discussed by [10][11][12][13][14][15] and are proven to be efficient in improving the accuracy and reduced the number of function evaluation. This approach is extended to various type by increasing the order, varies the step size and orders with extended stability region. All the methods, mentioned above, produced solutions at multiple points simultaneously at single iteration. Independent solution at selected grid points was generated using block method without overlapping. Despite the success of the block approach, there has been no guidance or recommendation on the choice of the most suitable numbers of interpolation points or the suitable order in improving the accuracy with lesser computational time when solving a very stiff problem. In deriving a good numerical method, the challenge is to increase the stability region while maintaining or even improving efficiency. Therefore, our objective is to construct a new class of hybrid methods by extending the existing block methods by [10] to address the setbacks associated with solving stiff ODEs. The improvement is obtained by considering [10] with a suitable hybrid points to find numerical solutions of (1).
The paper is organized as follows: in Section 2, we present the formulation of the hybrid method by presenting a review of classical block methods and hybrid block methods from the literature; Section 3 includes the analysis of the new method in terms of order, zero stability, and region of stability. Implementation of the method is discussed in Section 4. In Section 5 numerical tests are performed on first-order differential equation problems possessing stiffness and comparisons with MATLAB solvers are applied in providing numerical solutions for initial value problems of ordinary differential equations. Finally, conclusions are presented in Section 6.

Derivation of the Method
Various attempts have been proposed for solving (1) using block methods ranging from constant step size to variable step size, fully implicit and then to diagonally implicit methods, see [10,12,13]. We will illustrate the idea of computation in block by reviewing the theory of block methods available in the literature as follows.

Review of Block Methods
The idea of using block methods for solving ODEs was initially developed by [14,15] followed by [16]. Later, [17] proposed block methods in the form of predictor-corrector for solving special second-order ODEs of type non-stiff while [18] proposed parallel solution of solving (1). Adopting the definition by [17], the general form of a block method can be expressed as follows: The method derived in the form of (2) has the advantage of producing point solutions simultaneously at each application. Lagrange Interpolation was adopted by some authors to generate continuous linear multistep method, among them are [11][12][13].

Hybrid Block Method
Recently many authors have extended existing block method by adding hybrid points which is also referred as off-step points to find numerical solutions for (1). Some of the works on hybrid methods when solving (1) has been developed by see [19][20][21][22]. From [21], the focus is on the construction of explicit hybrid methods with nonnegative coefficients, which are a class of multistep methods incorporating a function evaluation at an off-step point. In [22], researchers derived a two-step hybrid methods for ′′ , , with oscillatory or periodic solutions, through the usage of the exponential fitting technique. Recently, [23] proposed a method by considering two intermediate points and the approximation of the true solution by an adequate polynomial and imposing collocation conditions. However, there are some drawbacks ranging from choosing the suitable order of the method and the number of off-step points to be included as hybrid methods.

Formulation of Hybrid Block
The method derived in this section is a hybrid block multistep method. The formulation is based on the existing BBDFs by adding off-step points in the formula that possessed stability property can be generated. The formulas are commonly derived by using an integrated form of the ODEs and replacing the integrals with interpolating polynomial. We consider the grid points, , , , and , with ℎ and three off-step-points / , / , / , / , and / , / as illustrated in Figure 1. The basic idea to formulate hybrid block backward differentiation formulas (HBBDF) is to replace , in (1) by an interpolating polynomial of degree . Assume the solution to (1) is approximated by the polynomial : The polynomial associated with interpolating points / , / , , , / , / , , , / , / and , is given by: Substituting ℎ in (4), we obtain: Generally, if we differentiate (3), yields the following: In our case, we differentiate (5) Upon substituting 1/2, 0,1/2, 1 respectively into (7) we obtained the coefficients for the HBBDF method as follows:

Analysis of the Proposed Method
In this section our objective is directed towards the discussion of the properties of the HBBDF method which includes order of the new method, zero stability and consistency which are the necessary conditions for convergence of the method. Region of absolute stability of the HBBDF will also be studied to determine its suitability of solving stiff ODEs.

Order
We start by defining the order and error constant as given in [24].

Definition 2. (Order)
The linear multistep method associated with the linear difference operator are said to be of order if ⋯ 0 and 0.

Definition 3. (Error constant)
The term is called the error constant and it implies that the local truncation error is given by: In order to check the consistency of the method, first, we will determine the order of the method. By rearranging the equation in (8), we get the following equation: Based on the matrix form of Equation (2), Equation (9) can be written as follows: We defined the linear difference operator associated with the hybrid BBDFs method in Equation (2) which is given by [20] as: where and are constants. Expanding ℎ 2 ⁄ and ′ ℎ 2 ⁄ in Equation (10) using Taylor series at , and collecting the coefficients of ℎ, we obtained: By letting the collected terms of , , ′′ in (11)    We obtained 0, hence the HBBDF is of order 5. Since the HBBDF method is of order 5 1, it is consistent by definition. The error constant is just above.

Convergence of the Derived Method
In this section we investigate the convergence of the HBBDFs. The convergence of a numerical method is generally determined by analysing consistency and zero stability given in the following theorem.

Theorem 1. (Convergence):
The necessary and sufficient conditions for the linear multistep method to be convergent are that it be both consistent and zero-stable.
In order to check the consistency of the method, we use the following definition given by [24]. □

Definition 4. (Consistent)
The linear multistep method is said to be consistent if it has the order of p 1. Hence, HBBDF in (8) is consistent. Next, we discuss the zero stability associated with HBBDFs.

Definition 5. (Zero stable)
The block method is zero stable provided the roots R , j 1 1 k of the first characteristics polynomial ρ R specified as: This satisfies R 1, and for those roots with R 1, the multiplicity must not exceed two.
We find the stability polynomial of HBBDF by applying linear test problem into Equation (8) Next, we transformed Equation (14) into matrix form as below: which is equivalent to .
Denote ℎ ℎ , the associated stability polynomial, This yields the roots 0,0, 1, and 0.00999. All of the roots have modulus less than one which satisfied the zero stable conditions given in definition 5. Thus, we conclude that HBBDF is zero stable. We have shown that the HBBDF method is consistent and zero stable. From theorem 1, we conclude the HBBDF converged.

Stability of Mathematical Components
Before we construct the region of absolute stability of the HBBDF method derived in Section 2, we give some definitions adopt from [25].

Definition 6. (A stable)
A numerical method is said to be A stable if the left-hand half-plane with Re h 0, h hλ, is contained by the absolute stability regions of the method used. This means that the stability region covers the entire negative left half plane.

Definition 7.
A α stable) A method is said to be A α stable if it is stable for all λ such that π α arg λ π α.
The boundary of the stability region is given by the set of points determined by R e , 0 θ 2π for which |R| 1. Note that the set of values of λ for which the method is stable is called the stability domain. Below we present the stability region which corresponds to the HBBDF method drawn in the hλ plane. The region of absolute stability is shown in Figure 2. The shaded region represents the set of points for which | | 1. Evidently, from Figure 2, this region of stability includes all the negative real axis which corresponds to an stable method. Next, we make comparison of stability region of HBBDF with the existing BBDF of the same order. The region for both methods are presented in Figure 3. Looking at the stability regions in Figure 3, we see that the HBBDF is stable while the BBDF(5) is stable. The absolute instability interval of both methods is tabulated in Table 1.

Implementation of HBBDF Method
In this section, we discussed the application of Newton's method in the implementation of the HBBDF method for integrating stiff ODEs. Let denote the iteration, denote the 1 iterative value of and denote as the differences between and 1 iteration values of as follows: The code developed will be using the PECE (predict-evaluate-correct-evaluate) mode. denote the predicted value of .

Numerical Results and Discussion
In this section, we present some numerical results in order to illustrate the performance of the new method in solving stiff ODEs. Three sets of stiff problems are tested under various step sizes 10 , 10 , and 10 . The numerical results will be compared with stiff solver ode15s in MATLAB and with the existing BBDF by [26] of the same order. For error calculations, the error formula is given by: (20) In (20), is the exact solution for the problem considered and is the computed solution.
Problem 1. Source [9] 100 1, 0 1, 0 10, 100    Tables 3, 5 and 7 gives the comparison of numerical solution at ten points equally spaced in the interval of integration to exact solution when ℎ 10 .         Here we observe that the graph of the approximated solution and the graph of the exact solution coincide with each other

Conclusions
In this paper we have derived a class of hybrid BBDF for the numerical integration of stiff ordinary differential equations. By restricting the method to three off-step points, the results presented in show that the stability region for the HBBDF with three off-step points are larger than the stability of the BBDF method without the off-step points even though both methods are of the same order. We also establish the convergence of the new method by investigate the zero-stability and consistency which are the necessary conditions for convergence. Our numerical experiments support the conclusion that HBBDF solve the problem with acceptable efficiency when compared to MATLAB solver, ode15s, especially when we are integrating at smaller step sizes. It is recommended to extend the HBBDF to a variable step size variable order, VSVO code.

Acknowledgments:
The authors are thankful to the referees for their valuable comments.

Conflicts of Interest:
The authors declare no conflict of interests regarding the publication of this paper.

HBBDF
Hybrid Block Backward Differentiation Formula of order 5 BBDF(5) Fifth Order Block Backward Differentiation Formula by [26] ode15s Variable order solver based on numerical differentiation formulas, NDFs. h Step size