#### Abstract

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.

#### 1. 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:

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), ;(2), 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 3-point block backward differentiation formula of order three (FI3BBDF(3)).

*FI2BBDF(3).* One has*FI3BBDF(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 3-point block extended backward differentiation formula of order three (FI3BEBDF(3)) are given in the following forms.

*FI2BEBDF(3).* One has*FI3BEBDF(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 , and are coefficients of and in the matrix formMethod (2) is defined to be diagonally implicit if is zero, whereas and 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.

#### 2. 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 and 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:wherefor each .

##### 2.1. Derivation of DI2BBDF(2)

DI2BBDF(2) will compute two approximated solutions and simultaneously in each block using two back values. This formula is derived using interpolating points to obtain the first formula of DI2BBDF(2): Replacing into (10) yieldsEquation (11) is differentiated once with respect to at the point . Evaluating givesThe same technique is applied for the second point of DI2BBDF(2). This formula is derived using as the interpolating points and producesSubstituting into (13) yieldsDifferentiating (14) with respect to at the point givesConsidering , 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.

##### 2.2. Derivation of DI2BBDF(3)

The first point of DI2BBDF(3) is derived using interpolating points , and we haveSubstituting into (17) givesEquation (18) is differentiated once with respect to at the point . Substituting will obtainThe derivation process continues for second point of DI2BBDF(3) using as the interpolating points. We haveReplacing into (20) producesThe resulting polynomial above is differentiated once with respect to at the point . Substituting will giveConsidering , the corrector formula of DI2BBDF(3) is given by

##### 2.3. Derivation of DI2BBDF(4)

The interpolating points are used to obtain the first point of DI2BBDF(4):We define and produceDifferentiating (25) once with respect to at the point and evaluating will produceThe derivation continues for the second point of formula by using as the interpolating points:Substituting into (27) will produceEquation (28) is differentiated once with respect to followed by substituting ; we haveTherefore, the corrector formula of DI2BBDF(4) is obtained as follows:

#### 3. 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 .

We consider the simplest test equationwhere the eigenvalues , , satisfy to analyze the stiff stability of the method. Substituting (31) into (16) will produce the following matrices form:whereLet and evaluating the determinant of from (32), the stability polynomial of DI2BBDF(2) is obtained as follows:The similar approach is applied to obtain the stability polynomial of DI2BBDF(3) and DI2BBDF(4). Substituting (31) into (23), we havewhereSubstituting (31) into (30) will givewhere

We compute the determinant of from (35) and (37) to obtain the stability polynomials of DI2BBDF(3) and DI2BBDF(4), respectively. Consider

Next, the boundary of the stability region will be determined by substituting into (34), (39), and (40). The graphs of stability region for all formulas are given in Figures 1, 2, and 3.

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.

#### 4. 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 and are the back values. Equation (41) is transformed into matrix form as follows:where Implementing Newton’s method to (42) producesIteration for (44) is given bywhere denotes the th iteration. Equation (45) can be rearranged into the following equation:Replacing (44) into (46) yieldswhere is the Jacobian matrix of with respect to .

##### 4.1. Newton Iteration of DI2BBDF(2)

Formula (16) is written in the form of (41):Applying (48) into (47) will yield matrix form as follows:with being the increment.

##### 4.2. Newton Iteration of DI2BBDF(3)

Rearranging formula (23) in the form of (41), we haveReplacing (50) into (47), this will yield matrix form as follows:with being the increment.

##### 4.3. Newton Iteration of DI2BBDF(4)

We rewrite formula (30) in the form of (41) and obtainSubstituting (52) into (47) will producewith being the increment.

All the formulas derived are implemented in predictor-corrector 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 and in PECE are as follows.(1)P (Predict): .(2)E (Evaluate): .(3)C (Correct): .(4)E (Evaluate): .To approximate the solution to , we apply the two-stage Newton type iteration. The iteration process for DI2BBDF is done as follows.(1)Compute the values for , where(2)Calculate the corrected value for with the value from Step .(3)Solve for second stage iteration.(4)The final values for are obtained from the second stage iteration of .The similar iteration process is applied for DI2BBDF(3) and DI2BBDF(4).

#### 5. 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.

*Problem 1 (linear).* One hasExact solution isInitial condition is(source: Ibrahim [8]).

*Problem 2 (linear).* One hasExact solutions areInitial conditions are(source: Musa et al. [13]).

*Problem 3 (linear).* One hasExact solution isInitial condition is(source: Ibrahim et al. [10]).

*Problem 4 (linear).* One hasExact solutions areInitial conditions are(source: Musa et al. [14]).

*Problem 5 (nonlinear).* One hasExact solution isInitial condition is(source: Musa et al. [14]).

*Problem 6 (nonlinear).* One hasExact solution isInitial condition is(source: Musa et al. [14]).

#### 6. Numerical Results

The performance of the derived methods is compared with the existing methods in terms of maximum error and execution time. We consider , , and as the step size, . Tables 1 and 2 present the performance comparison of DI2BBDF(2), DI2BBDF(3), and DI2BBDF(4) with FI2BBDF(3) and FI2BEBDF(3), whereas Tables 3–6 exhibit the comparison of proposed methods with FI3BBDF(3) and FI3BEBDF(3). In addition, the graphs of log(MAXE) against log() are illustrated as shown in Figures 5–10. The following notations are used in Tables 1–6: : step size; MAXE: maximum error; TIME: time execution using high performance computer (HPC); DI2BBDF(2): diagonally implicit 2-point block backward differentiation formulas of order 2; DI2BBDF(3): diagonally implicit 2-point block backward differentiation formulas of order 3; DI2BBDF(4): diagonally implicit 2-point block backward differentiation formulas of order 4; FI2BBDF(3): fully implicit 2-point block backward differentiation formulas of order 3; FI2BEBDF(3): fully implicit 2-point block extended backward differentiation formulas of order 3; FI3BBDF(3): fully implicit 3-point block backward differentiation formulas of order 3; FI3BEBDF(3): fully implicit 3-point block extended backward differentiation formulas of order 3.

#### 7. Discussion

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

##### 7.1. 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 implicit method has less differentiation coefficients in order to prevent the cumulative error. The graphs in Figures 5–10 depict the scales of maximum error versus step size, , for the proposed methods as compared to existing methods. There is a slight drop in accuracy for the proposed methods as the order increases. This is due to the increase in cumulative errors during the computation as more interpolating points are used. Among the methods of order three, DI2BBDF(3) is more accurate compared with FI2BBDF(3), FI3BBDF(3), FI2BEBDF(3), and FI3BEBDF(3). For all test problems, we conclude that the accuracy of the proposed method increases as the step size becomes smaller.

##### 7.2. Computational Time

In Tables 1-2, it can be seen that the execution time taken by DI2BBDF(2), DI2BBDF(3), and DI2BBDF(4) is comparable with that of FI2BBDF(3) and FI2BEBDF(3). The computational time increases as the step size becomes smaller. However, Tables 3–6 show that the proposed methods compute faster than FI3BBDF(3) and FI3BEBDF(3). This could be justified by the fact that the 3-point method has to perform extra computation on the Jacobian matrix since it involves a 3 by 3 matrix, while the 2-point method has to compute Jacobian matrix of dimension 2 by 2 only. Among the proposed methods, DI2BBDF(2) computes faster than DI2BBDF(3) and DI2BBDF(4). In fact, the number of back values involved in every formula will affect the computational time. Since DI2BBDF(4) has more back values than DI2BBDF(3) and DI2BBDF(2), it requires extra time to compute the approximated solutions.

#### 8. Conclusion

Research conducted in this paper shows the capability of constructing the diagonally implicit BBDF method with A-stable properties. From the results obtained via the numerical experiment, we can conclude that the proposed method serves the purpose of significant alternative numerical method for solving linear and nonlinear stiff IVPs occurring in the fields of engineering and physical sciences.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgment

This research was supported by the Institute for Mathematical Research (INSPEM), Department of Mathematics, Universiti Putra Malaysia (UPM), under Fundamental Grant Scheme (Project code 02-02-13-1380FR).