Derivation of Diagonally Implicit Block Backward Differentiation Formulas for Solving Stiff Initial Value Problems

The diagonally implicit 2-point block backward differentiation formulas (DI2BBDF) of order two, order three, and order four are derived for solving stiff initial value problems (IVPs). The stability properties of the derived methods are investigated. The implementation of the method using Newton iteration is also discussed. The performance of the proposed methods in terms of maximum error and computational time is compared with the fully implicit block backward differentiation formulas (FIBBDF) and fully implicit block extended backward differentiation formulas (FIBEBDF). The numerical results show that the proposed method outperformed both existing methods.


Introduction
Many scientific and engineering problems which arise in real-life applications are in the form of ordinary differential equations (ODEs), where the analytic solution is unknown.The general form of first order ODEs is given in the following form: =  (, ) ,  () =  0 ,  ≤  ≤ . (1) In the early 1950s, Curtiss and Hirschfelder [1] realized that there is an important class of ODEs which is known as stiff initial value problems (IVPs).There are various definitions of stiffness given in the literature.Generally, stiff problems are problems where certain implicit methods perform better than explicit ones.For simplicity, we choose the definition of stiff problem given by Lambert [2].Definition 1.The system of ( 1) is said to be stiff if (1) Re(  ) < 0,  = 1, 2, . . ., ; (2) max  | Re(  )| ≫ min  | Re(  )|, where   are the eigenvalues of the Jacobian matrix,  = (/).
Much research has been done by the scientific community on developing numerical methods which permit an approximate solution to (1).The most commonly used numerical method is block method.This classical method is introduced by Milne [3] to compute previous -blocks and calculate the current block where each block contains -point.The -block -point method is given by a matrix finite difference equation of the form: where   and   are properly chosen × matrix coefficients.This method is extended by Shampine and Watts [4] with convergence and stability properties of one step block implicit method.Then Fatunla [5] proposed block -point method followed by Majid et al. [6,7] with the 2-point block methods.However, the most widely used multistep method for solving stiff ODEs is block backward differentiation formulas (BBDF).This method has been claimed by Ibrahim [8] to be one of the suitable numerical methods for solving stiff IVPs.Furthermore, Ibrahim et al. [9,10] proposed the fully implicit -point block backward differentiation formulas (FIBBDF).The following equations represent the formulas of fully implicit 2-point block backward differentiation formulas of order three (FI2BBDF(3)) and fully implicit 3point block backward differentiation formula of order three (FI3BBDF(3)).
FI2BBDF (3).One has In a related study, Ibrahim et al. [11] plotted the stability region of (3).Since all region in the left half plane is in the stability region, the FI2BBDF(3) is A-stable and suitable for solving stiff problems.Nasir et al. [12] extended the order of formula (3) which is called fifth order 2-point BBDF for solving first order stiff ODEs.Recently, Musa et al. [13,14] modified formulas (3) and (4) to compute more than one solution value per step using extra future point.This method is called fully implicit block extended backward differentiation formulas (FIBEBDF).The formulas of fully implicit 2-point block extended backward differentiation formula of order three (FI2BEBDF(3)) and fully implicit 3point block extended backward differentiation formula of order three (FI3BEBDF(3)) are given in the following forms.

FI2BEBDF(3). One has
Numerical results in the literature showed that the FIBEBDF performed better as compared to FIBBDF in terms of accuracy.Unfortunately, the execution time of the FIBEBDF is slower than FIBBDF.This is due to the fact that FIBEBDF has an extra future point which requires more computation time.In order to gain an efficient numerical approximation in terms of accuracy and computational time, the diagonally implicit method must be considered.The study of diagonally implicit for multistep method attracted several researchers such as Alexander [15], Ababneh et al. [16], and Ismail et al. [17].However, Lambert [2] stated that there is some confusion over nomenclature to identify the diagonally implicit.Some authors use the term of diagonally implicit to describe any semi-implicit method.Therefore, the definition of diagonally implicit block method is given by Majid and Suleiman [18] as follows.
Definition 2. Method (2) is defined to be diagonally implicit if the coefficients of the upper-diagonal entries are zero.
From this motivation, we established Definition 2 by introducing the definition of diagonally implicit 2-point BBDF method as follows.Definition 3. We consider that  11 ,  12 ,  21 , and  22 are coefficients of  +1 and  +2 in the matrix form Method ( 2) is defined to be diagonally implicit if  12 is zero, whereas  11 and  22 are equal.
Therefore, the main purpose of this paper is to develop a new diagonally implicit multistep method for solving stiff ODEs.This paper is organized as follows: in Section 2, the diagonally implicit two-point block backward differentiation formulas (DI2BBDF) will be derived.Next, the stability properties of the derived methods are analyzed in Section 3. Section 4 discusses the implementation of the methods using Newton iteration.Standard test problems are selected in Section 5, whereas the performance of the proposed method is shown in Section 6. Finally in Section 7 some conclusions are given.

Formulation of the Method
In this section, we will derive the DI2BBDF of order two, order three, and order four with constant step size to compute the approximated solutions at  +1 and  +2 concurrently.Contrary to the fully implicit method that has been proposed by Ibrahim et al. [11], the first point of diagonally implicit formula has one less interpolating point.The derivation using polynomial   () of degree  in terms of Lagrange polynomial is defined as follows: where for each  = 0, 1, . . ., . (2).DI2BBDF(2) will compute two approximated solutions  +1 and  +2 simultaneously in each block using two back values.This formula is derived using interpolating points  −1 , . . .,  +1 to obtain the first formula  +1 of DI2BBDF(2):

Derivation of DI2BBDF
Replacing  = ℎ +  +1 into (10) yields Equation ( 11) is differentiated once with respect to  at the point  =  +1 .Evaluating  = 0 gives The same technique is applied for the second point  +2 of DI2BBDF(2).This formula is derived using  −1 , . . .,  +2 as the interpolating points and produces Substituting  = ℎ +  +2 into (13) yields Differentiating ( 14) with respect to  at the point  =  +2 gives Considering ℎ +1,+2 =   ( +1,+2 ), the corrector formula of DI2BBDF( 2) is given as follows: The order of the method is distinguished by the number of back values contained in the formulas.Adopting a similar approach as the derivation of DI2BBDF of order two, we will construct the DI2BBDF of orders three and four with different number of interpolating points.
Mathematical Problems in Engineering

Stability Analysis
In this section, we will plot the graph of stability for DI2BBDF(2), DI2BBDF(3), and DI2BBDF(4) using Mathematica software.Based on Dahlquist [19], the linear multistep method (LMM) is able to solve stiff problems if it satisfies the following definition.
Definition 4. The LMM is A-stable if its region of absolute stability contains the whole of the left-hand half-plane Re(ℎ) < 0.
In Figures 1, 2, and 3, we observed that the intervals of unstable region for DI2BBDF(2), DI2BBDF(3), and  DI2BBDF( 4) are [0, 5.33333], [0, 8.66667], and [0, 13.86667], respectively.All the graphs of stability regions are combined as shown in Figure 4.The region outside the green line is the stable region of DI2BBDF(2); the region outside the red line is the stable region of DI2BBDF(3); and the region outside the blue line is the stable region of DI2BBDF(4).
Clearly, the unstable region becomes larger when the order of the method increases.From Definition 4, DI2BBDF(2) is A-stable, while DI2BBDF(3) and DI2BBDF(4) are almost A-stable since the stability region covers the entire negative half plane.Therefore, we can conclude that the proposed method is suitable for solving stiff problems.

Implementation of the Method
In this section, the proposed methods will be implemented using Newton iteration.We begin by converting ( 16), (23), and (30) in general form as follows: where  1 and  2 are the back values.Equation ( 41) is transformed into matrix form as follows: where Implementing Newton's method to (42) produces Iteration for ( 44) is given by where  (+1) +1 denotes the ( + 1)th iteration.Equation ( 45) can be rearranged into the following equation: Replacing ( 44) into ( 46) yields where (/)( () +1,+2 ) is the Jacobian matrix of  with respect to . (2).Formula ( 16) is written in the form of (41):

Newton Iteration of DI2BBDF(4).
We rewrite formula (30) in the form of (41) and obtain Substituting ( 52) into (47) will produce with  (+1) +1,+2 =  (+1) +1,+2 −  () +1,+2 being the increment.All the formulas derived are implemented in predictorcorrector computation which is symbolized as PECE mode.P and C indicate one application of the predictor and the corrector, respectively, and E indicates one evaluation of the function , given  and .The approximation calculations for  +1 and  +2 in PECE are as follows.
To approximate the solution to  () +1,+2 , we apply the two-stage Newton type iteration.The iteration process for DI2BBDF(2) is done as follows.

Test Problems
In this section, linear and nonlinear stiff problems are tested using C programming to examine the efficiency and reliability of the proposed method.The following problems are commonly found in engineering and physical sciences, particularly in the studies of vibrations, electrical circuits, and chemical reaction.

Numerical Results
The performance of the derived methods is compared with the existing methods in terms of maximum error and execution time.We consider 10 −2 , 10 −4 , and 10 −6 as the step size, ℎ.Tables 1 and 2

Discussion
This section is divided into discussion of maximum error and computational time.

Maximum Error.
From the numerical results in Tables 1-6, we observe that DI2BBDF(2), DI2BBDF(3), and DI2BBDF(4) outperformed the existing methods in terms of maximum error.This is expected because the diagonally

Table 1 :
The accuracy for problem 1.

Table 2 :
The accuracy for problem 2.

Table 3 :
The accuracy for problem 3.

Table 4 :
The accuracy for problem 4.

Table 5 :
The accuracy for problem 5.

Table 6 :
The accuracy for problem 6.