Research Article  Open Access
Derivation of Diagonally Implicit Block Backward Differentiation Formulas for Solving Stiff Initial Value Problems
Abstract
The diagonally implicit 2point 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 reallife 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 2point 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 2point 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 hasFI3BBDF(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 Astable and suitable for solving stiff problems. Nasir et al. [12] extended the order of formula (3) which is called fifth order 2point 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 2point 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 hasFI3BEBDF(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 semiimplicit 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 upperdiagonal entries are zero.
From this motivation, we established Definition 2 by introducing the definition of diagonally implicit 2point 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 twopoint 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 Astable if its region of absolute stability contains the whole of the lefthand halfplane .
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 Astable, while DI2BBDF(3) and DI2BBDF(4) are almost Astable 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 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 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 twostage 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 2point block backward differentiation formulas of order 2;â€‰DI2BBDF(3): diagonally implicit 2point block backward differentiation formulas of order 3;â€‰DI2BBDF(4): diagonally implicit 2point block backward differentiation formulas of order 4;â€‰FI2BBDF(3): fully implicit 2point block backward differentiation formulas of order 3;â€‰FI2BEBDF(3): fully implicit 2point block extended backward differentiation formulas of order 3;â€‰FI3BBDF(3): fully implicit 3point block backward differentiation formulas of order 3;â€‰FI3BEBDF(3): fully implicit 3point block extended backward differentiation formulas of order 3.
