# Diagonally Implicit Block Backward Differentiation Formulas for Solving Ordinary Differential Equations

Authors: I. S. M. Zawawi; Z. B. Ibrahim; F. Ismail; Z. A. Majid
Journal: International Journal of Mathematics and Mathematical Sciences
 (2012)

---

## Abstract

This paper focuses on the derivation of diagonally implicit two-point block backward differentiation formulas (DI2BBDF) for solving first-order initial value problem (IVP) with two fixed points. The method approximates the solution at two points simultaneously. The implementation and the stability of the proposed method are also discussed. A performance of the DI2BBDF is compared with the existing methods.


---

## Body

## 1. Introduction

Most real life problems that arise in various fields of study such as engineering or science are constructed as mathematical models before they are solved. These models often lead to differential equations. A differential equation can simply be defined as an equation that contains a derivative. In this paper, we are concerned with the numerical solution of IVP with two fixed points of the general form(1.1)y′=(y-v1)(y-v2)g(y),
given initial values y(xn)=yn, where v1<v2∈ℝ and g(y)≠0 is a bounded real-valued function with continuous derivatives. See; [1] how the problem may be simplified and we can assume that the fixed points are y1(x)=0 and y2(x)=1. The aim of this paper is to introduce a new block method which is called diagonally implicit two point block backward differentiation formulas. Basically, block method is a method which is used to compute previous  k  blocks and to calculate the current block where each block contain r points. Ibrahim et al. [2] gave the general form for r point k block as
(1.2)∑j=0kAjyn+j=h∑j=0kBjfn+j,
where  Aj  and Bj are r  by  r  matrices.The rapid growth of the studies on the block methods for solving ODEs contribute to the competition in developing and deriving an accurate method for solving many types of ODEs. Among the earliest research on block methods was proposed by Shampine and Watts [3] with block implicit one step method, Cash [4] with modified extended backward differential formulas, Chu and Hamilton [5] with multi-block methods, Majid and Suleiman [6] with fully implicit block method for solving ODEs. The most recent study of block method is on computing block approximation yn+1,yn+2…yn+k  known as block backward differentiation formulas (BBDF) presented by Ibrahim et al. [7]. Diagonally implicit for multistep method are discussed by a few researchers such as Majid and Suleiman [8] with 4-point diagonally implicit block method and Alexander [9] with diagonally implicit Runge-Kutta method for solving stiff ODEs. In this paper, we interested to compare the accuracy of the proposed method with the one-step predictor-corrector method (1SPCM) and P-C ADAMS method studied by Aguiar and Ramos [1].


## 2. Formulation of Diagonally Implicit BBDF

In Figure1, diagonally implicit 2-point BBDF (DI2BBDF) will create two new equally spaced solution values simultaneously. Majid and Suleiman [8] stated that the method is called diagonally implicit because the coefficients of the upper triangular matrix entries are zero.Figure 1
Diagonally Implicit 2-point block method of constant step size.We will derive the formula to compute the approximate values ofyn+1 and yn+2 simultaneously with two previous back values xn-1 and xn. The formula was derived by using Lagrange interpolation polynomial which compute the solutions at xn+1 and xn+2 separately. Consider the polynomial Pk(x) of degree k which interpolates the valuesyn, yn-1,⋯, yn-k+1 of a function f at interpolating points xn, xn-1,⋯, xn-k+1  in terms of Lagrange polynomial defined as follows:
(2.1)Pk(x)=∑j=0kLk,j(x)f(xn+1-j),
where
(2.2)Lk,j(x)=∏i=0i≠jk(x-xn+1-i)(xn+1-j-xn+1-i)
for each j=0,1,…,k.Foryn+1, we defined s=(x-xn+1)/h and replace f(x,y) by polynomial (2.1) which interpolates only the values yn-1, yn, and yn+1 at the interpolating points xn-1,xn,…,xn+1.
(2.3)P(xn+1+sh)=(xn+1+sh-xn)(xn+1+sh-xn+1)(-h)(-2h)yn-1+(xn+1+sh-xn-1)(xn+1+sh-xn+1)(h)(-h)yn+(xn+1+sh-xn-1)(xn+1+sh-xn)(2h)(h)yn+1.
Thus, differentiating the resulting polynomial once with respect to s  at the point x=xn+1 and evaluating at s=0  gives the following:
(2.4)P′(x)=P′(xn+1)=hfn+1=12yn-1-2yn+32yn+1.Then, we interpolate the valuesyn-1, yn, yn+1, and yn+2 at the interpolating points xn-1, xn, xn+1,…,xn+2 and differentiating the resulting polynomial once with respect to s at the point x=xn+2  and substituting s=0 yield
(2.5)P'(x)=P'(xn+2)=hfn+2=-13yn-1+32yn-3yn+1+116yn+2.Therefore, the corrector formulae atxn+1   and  xn+2  are given by
(2.6)yn+1=-13yn-1+43yn+23hfn+1,yn+2=211yn-1-911yn+1811yn+1+611hfn+2.


## 3. Stability of Diagonally Implicit BBDF

Many numerical methods depend on speed of convergence, computational expense and accuracy for the solution of initial value problems (IVPs). In this section, our aim is to identify the properties of linear stability of DI2BBDF. According to Shampine and Watts [3], the most critical limitation of block method is due to the existence of stability problems. The basic definition of zero stable for block method described by Ibrahim et al. [10]. The block method is said to be zero-stable if the roots Rj,j=1(1)k, of the first characteristic polynomial ρ(R)=det[∑i=0kAiRk-i]=0,A0=-I, satisfy |Rj|≤1.  If one of the roots is +1, we call this root the principal root of ρ(R).The linear stability properties are determined through application of the standard linear test problem(3.1)y′=λy,λ<0,λcomplex.
Hence, we will construct the formulas that have been derived in the previous section, which is (2.6). Application of the standard linear test problem (3.1) to (2.6) gives the following
(3.2)yn+1-23λhyn+1=-13yn-1+43yn.yn+2-611λhyn+2=211yn-1-911yn+1811yn+1.
Setting h^=λh, we write (3.2)  in the matrix coefficient specified as
(3.3)A=[(1-23h^)00(1-611h^)],Ym=[yn+1yn+2],B=[-1343211-911],Ym-1=[yn-1yn].
The first characteristic polynomial of (3.2) is given by
(3.4)ρ(t)=det[tA-B]=t2-4033t2h^-3433t+411t2h^2-811th^+13.By solvingρ(t)=0  and h^=0, we can determine the zero stable and yield t=1 and t=1/33. Thus the diagonally implicit BBDF is zero stable. Since one of the roots is +1, we call this root as principal root.


## 4. Implementation of the Method

In this section, we were obtained the calculation ofyn+1 and yn+2 using the Newton iteration. The DI2BBDF can be written in general form as
(4.1)F1=yn+1-23hfn+1-μ1F2=yn+2-1811yn+1-611hfn+2-μ2,
where μ1 and μ2 are the  backvalues.To specify the iteration, the following notation is introduced.yn+1(i+1) will denote the (i+1)th iterative value of yn+1, and let  en+1(i+1)=yn+1(i+1)-yn+1(i) and en+2(i+1)=yn+2(i+1)-yn+2(i). Newton iteration takes the form as below:
(4.2)yn+1(i+1)=yn+1(i)-(F1[yn+1(i)])(F′1[yn+1(i)])-1,yn+2(i+1)=yn+2(i)-(F2[yn+2(i)])(F2′[yn+2(i)])-1.SubstitutingF1F1′, F2 and F2′ take the form
(4.3)yn+1(i+1)-yn+1(i)=-(yn+1(i)-23hfn+1(i)-μ1)(1-23h∂fn+1∂yn+1)-1,(1-23h∂fn+1∂yn+1)(yn+1(i+1)-yn+1(i))=-y(i)n+1+23hfn+1(i)+μ1.Similarly for the second point,yn+2 gives
(4.4)yn+2(i+1)-yn+2(i)=-(yn+2(i)-1811yn+1(i)-611hfn+2(i)-μ2)(1-611h∂fn+2∂yn+2)-1,(yn+2(i+1)-yn+2(i))(1-611h∂fn+2∂yn+2)=-yn+2(i)+1811y(i)n+1+611hfn+2(i)+μ2.Hence, equation (4.3) and (4.4) in matrix form
(4.5)[(1-23h∂fn+1∂yn+1)01811(1-611h∂fn+2∂yn+2)][en+1(i+1)en+2(i+1)]=[-101811-1][yn+1(i)yn+2(i)]+h[2300611][fn+1(i)fn+2(i)]+[μ1μ2].Therefore, the approximation values ofyn+1and yn+2 are
(4.6)yn+1(i+1)=yn+1(i)+en+1,(i+1)yn+2(i+1)=yn+2(i)+en+2(i+1).


## 5. Numerical Example

In this section, an autonomous problem with two fixed points considered for the purpose of validating the numerical results:(5.1)y′(x)=y(x)(y(x)-1)y(x)-2.
The solution is  y(0)=y0.Exact solution for0<y0<1 is
(5.2)y(x)=ex/2(ex/2y02-y0exy02-4y0+4)2(y0-1).
See  Aguiar and Ramos [1].We have tested the problem fory0=0.1 and consider the step size, h=2-j where j=2,3,4,5,6,7,8 with interval [0,20]. The accuracy for all methods in terms of maximum error is shown in Table 1.Table 1
Numerical results for problem (5.1) for different step sizes, h=2-j.
j
Method
MAXE
1SPCM
4.3821
E
-
5
2
P-C ADAMS
1.2533
E
-
3
DI2BBDF
7.4651
E
-
3
1SPCM
1.1130
E
-
5
3
P-C ADAMS
3.2000
E
-
4
DI2BBDF
1.9472
E
-
3
1SPCM
2.8070
E
-
6
4
P-C ADAMS
8.0949
E
-
5
DI2BBDF
4.9778
E
-
4
1SPCM
7.0478
E
-
7
5
P-C ADAMS
2.0362
E
-
5
DI2BBDF
1.2576
E
-
4
1SPCM
1.7656
E
-
7
6
P-C ADAMS
5.1065
E
-
6
DI2BBDF
3.1612
E
-
5
1SPCM
4.4189
E
-
8
7
P-C ADAMS
1.2786
E
-
6
DI2BBDF
7.9239
E
-
6
1SPCM
1.1053
E
-
8
8
P-C ADAMS
3.1991
E
-
7
DI2BBDF
1.9836
E
-
6The notation used in the table takes the following meaning:h: step size,MAXE: maximum error,DI2BBDF: diagonally implicit 2-point block backward differentiation formulas method,1SPCM: one-step predictor-corrector method,P-C ADAMS: explicit Euler as a predictor and trapezoidal rule as a corrector.The maximum error is defined as follows: Maximum error =max1≤i≤NS(|yi-y(x)i|)  where NS is the number of steps.In Figure2, the numerical results obtained with the proposed method and existing methods are plotted. The exact and approximate solutions using h=0.4 are also plotted in Figure 3.Figure 2
The graph oflog10 (MAXE) against  j.Figure 3
The graph of exact and approximate solutions againstx∈[0,20]  using h=0.4.


## 6. Conclusion

In this work, we have applied DI2BBDF as a numerical method for solving ordinary initial value problem with two fixed points. By reducing the step sizeh, we observe that the value of maximum errors for DI2BBDF gets smaller. Furthermore, by comparing our results with the solutions obtained by the existing methods, we can conclude that DI2BBDF yields less accurate approximations. However, for the case of high dimension problem, DI2BBDF still performs well although as known that block backward differentiation formula is more effective for solving stiff initial value problems. In Figure 3, we can observe that both exact and approximate solutions of DI2BBDF are still increase towards the stable fixed point y=1. For future research, we map out the strategy to improve the performance of this method by extending the order of the formula and modify the formula in a variable-step size for solving differential initial value problems with two fixed points.



---

*Source: 767328-2012-11-13.xml*