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

Authors: Iskandar Shah Mohd Zawawi; Zarina Bibi Ibrahim; Khairil Iskandar Othman
Journal: Mathematical Problems in Engineering
 (2015)

---

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


---

## Body

## 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:(1)
y
′
=
f
x
,
y
,
y
a
=
y
0
,
a
≤
x
≤
b
.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
λ
t
<
0, t
=
1,2
,
…
,
m;
(2)
max
⁡
t
Re
λ
t
≫
min
⁡
t
Re
λ
t, where λ
t are the eigenvalues of the Jacobian matrix, J
=
∂
f
/
∂
y.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 k-blocks and calculate the current block where each block contains r-point. The k-block r-point method is given by a matrix finite difference equation of the form:(2)
∑
j
=
0
k
A
j
y
n
+
j
=
h
∑
j
=
0
k
B
j
f
n
+
j
,where A
j and B
j are properly chosen r
×
r 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 r-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 r-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(3)
y
n
+
1
=
-
1
3
y
n
-
1
+
2
y
n
-
2
3
y
n
+
2
+
2
h
f
n
+
1
,
y
n
+
2
=
2
11
y
n
-
1
-
9
11
y
n
+
18
11
y
n
+
1
+
6
11
h
f
n
+
2
.
FI3BBDF(3). One has(4)
y
n
+
1
=
1
10
y
n
-
2
-
3
4
y
n
-
1
+
3
y
n
-
3
2
y
n
+
2
+
3
20
y
n
+
3
+
3
h
f
n
+
1
,
y
n
+
2
=
-
3
65
y
n
-
2
+
4
13
y
n
-
1
-
12
13
y
n
+
24
13
y
n
+
1
-
12
65
y
n
+
3
+
12
13
h
f
n
+
2
,
y
n
+
3
=
12
137
y
n
-
2
-
75
137
y
n
-
1
+
200
137
y
n
-
300
137
y
n
+
1
+
300
65
y
n
+
2
+
60
137
h
f
n
+
3
.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(5)
y
n
+
1
=
1
9
y
n
-
1
-
y
n
+
17
9
y
n
+
2
-
2
h
f
n
+
1
-
2
3
h
f
n
+
2
,
y
n
+
2
=
17
197
y
n
-
1
-
99
197
y
n
+
279
197
y
n
+
1
+
150
197
h
f
n
+
2
-
18
197
h
f
n
+
3
.
FI3BEBDF(3). One has(6)
y
n
+
1
=
-
1
80
y
n
-
2
+
1
80
y
n
-
1
-
3
4
y
n
+
25
16
y
n
+
2
+
3
40
y
n
+
3
-
3
2
h
f
n
+
1
-
3
4
h
f
n
+
2
,
y
n
+
2
=
-
3
25
y
n
-
2
+
y
n
-
1
-
4
y
n
+
12
y
n
+
2
-
197
25
y
n
+
3
+
12
h
f
n
+
2
+
12
5
h
f
n
+
3
,
y
n
+
3
=
394
14919
y
n
-
2
-
2925
14919
y
n
-
1
+
9600
14919
y
n
-
18700
14919
y
n
+
1
+
26550
14919
y
n
+
2
+
8820
14919
h
f
n
+
3
-
600
14919
h
f
n
+
4
.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 Definition2 by introducing the definition of diagonally implicit 2-point BBDF method as follows.Definition 3.
We consider thata
11
,
a
12
,
a
21, and a
22 are coefficients of y
n
+
1 and y
n
+
2 in the matrix form(7)
a
11
a
12
a
21
a
22
y
n
+
1
y
n
+
2
.Method (2) is defined to be diagonally implicit if a
12 is zero, whereas a
11 and a
22are 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 Section2, 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 aty
n
+
1 and y
n
+
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 P
k
(
x
) of degree k in terms of Lagrange polynomial is defined as follows:(8)
P
k
x
=
∑
j
=
0
k
L
k
,
j
x
f
x
n
+
1
-
j
,where(9)
L
k
,
j
x
=
∏
i
=
0
i
≠
j
k
(
x
-
x
n
+
1
-
i
)
(
x
n
+
1
-
j
-
x
n
+
1
-
i
)for each j
=
0,1
,
…
,
k.


### 2.1. Derivation of DI2BBDF(2)

DI2BBDF(2) will compute two approximated solutionsy
n
+
1 and y
n
+
2 simultaneously in each block using two back values. This formula is derived using interpolating points x
n
-
1
,
…
,
x
n
+
1 to obtain the first formula y
n
+
1 of DI2BBDF(2): (10)
P
x
=
x
-
x
n
x
-
x
n
+
1
x
n
-
1
-
x
n
x
n
-
1
-
x
n
+
1
y
n
-
1
+
x
-
x
n
-
1
x
-
x
n
+
1
x
n
-
x
n
-
1
x
n
-
x
n
+
1
y
n
+
x
-
x
n
-
1
x
-
x
n
x
n
+
1
-
x
n
-
1
x
n
+
1
-
x
n
y
n
+
1
.Replacing x
=
s
h
+
x
n
+
1 into (10) yields(11)
P
x
n
+
1
+
s
h
=
s
h
+
h
s
h
-
h
-
2
h
y
n
-
1
+
2
h
+
s
h
s
h
h
-
h
y
n
+
2
h
+
s
h
h
+
s
h
2
h
h
y
n
+
1
.Equation (11) is differentiated once with respect to s at the point x
=
x
n
+
1. Evaluating s
=
0 gives(12)
P
′
x
n
+
1
=
1
2
y
n
-
1
-
2
y
n
+
3
2
y
n
+
1
.The same technique is applied for the second point y
n
+
2 of DI2BBDF(2). This formula is derived using x
n
-
1
,
…
,
x
n
+
2 as the interpolating points and produces(13)
P
x
=
x
-
x
n
x
-
x
n
+
1
x
-
x
n
+
2
x
n
-
1
-
x
n
x
n
-
1
-
x
n
+
1
x
n
-
1
-
x
n
+
2
y
n
-
1
+
x
-
x
n
-
1
x
-
x
n
+
1
x
-
x
n
+
2
x
n
-
x
n
-
1
x
n
-
x
n
+
1
x
n
-
x
n
+
2
y
n
+
x
-
x
n
-
1
x
-
x
n
x
-
x
n
+
2
x
n
+
1
-
x
n
-
1
x
n
+
1
-
x
n
x
n
+
1
-
x
n
+
2
y
n
+
1
+
x
-
x
n
-
1
x
-
x
n
x
-
x
n
+
1
x
n
+
2
-
x
n
-
1
x
n
+
2
-
x
n
x
n
+
2
-
x
n
+
1
y
n
+
2
.Substituting x
=
s
h
+
x
n
+
2 into (13) yields(14)
P
x
n
+
2
+
s
h
=
s
h
+
2
h
s
h
+
h
s
h
-
h
-
2
h
-
3
h
y
n
-
1
+
s
h
+
3
h
s
h
+
h
s
h
h
-
h
-
2
h
y
n
+
s
h
+
3
h
s
h
+
2
h
s
h
2
h
h
-
h
y
n
+
1
+
(
s
h
+
3
h
)
(
s
h
+
2
h
)
s
h
+
h
(
3
h
)
(
2
h
)
h
y
n
+
2
.Differentiating (14) with respect to s at the point x
=
x
n
+
2 gives(15)
P
′
x
n
+
2
=
-
1
3
y
n
-
1
+
3
2
y
n
-
3
y
n
+
1
+
11
6
y
n
+
2
.Considering h
f
n
+
1
,
n
+
2
=
P
′
x
n
+
1
,
n
+
2, the corrector formula of DI2BBDF(2) is given as follows:(16)
y
n
+
1
=
-
1
3
y
n
-
1
+
4
3
y
n
+
2
3
h
f
n
+
1
,
y
n
+
2
=
2
11
y
n
-
1
-
9
11
y
n
+
18
11
y
n
+
1
+
6
11
h
f
n
+
2
.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 pointy
n
+
1 of DI2BBDF(3) is derived using interpolating points x
n
-
2
,
…
,
x
n
+
1, and we have(17)
P
x
=
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
1
x
n
-
2
-
x
n
+
1
x
n
-
2
-
x
n
x
n
-
2
-
x
n
-
1
y
n
-
2
+
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
2
x
n
-
1
-
x
n
+
1
x
n
-
1
-
x
n
x
n
-
1
-
x
n
-
2
y
n
-
1
+
x
-
x
n
+
1
x
-
x
n
-
1
x
-
x
n
-
2
x
n
-
x
n
+
1
x
n
-
x
n
-
1
x
n
-
x
n
-
2
y
n
+
x
-
x
n
x
-
x
n
-
1
x
-
x
n
-
2
x
n
+
1
-
x
n
x
n
+
1
-
x
n
-
1
x
n
+
1
-
x
n
-
2
y
n
+
1
.Substituting x
=
s
h
+
x
n
+
1 into (17) gives(18)
P
s
h
+
x
n
+
1
=
s
h
s
h
+
h
s
h
+
2
h
-
3
h
-
2
h
-
h
y
n
-
2
+
s
h
s
h
+
h
s
h
+
3
h
-
2
h
-
h
h
y
n
-
1
+
s
h
s
h
+
2
h
s
h
+
3
h
-
h
h
2
h
y
n
+
s
h
+
h
s
h
+
2
h
s
h
+
3
h
h
2
h
3
h
y
n
+
1
.Equation (18) is differentiated once with respect to s at the point x
=
x
n
+
1. Substituting s
=
0 will obtain(19)
P
′
x
n
+
1
=
11
6
y
n
+
1
-
1
3
y
n
-
2
+
3
2
y
n
-
1
-
3
y
n
.The derivation process continues for second point y
n
+
2 of DI2BBDF(3) using x
n
-
3
,
…
,
x
n
+
1 as the interpolating points. We have(20)
P
x
=
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
1
x
n
-
2
-
x
n
+
1
x
n
-
2
-
x
n
x
n
-
2
-
x
n
-
1
y
n
-
2
+
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
2
x
n
-
1
-
x
n
+
1
x
n
-
1
-
x
n
x
n
-
1
-
x
n
-
2
y
n
-
1
+
x
-
x
n
+
1
x
-
x
n
-
1
x
-
x
n
-
2
x
n
-
x
n
+
1
x
n
-
x
n
-
1
x
n
-
x
n
-
2
y
n
+
x
-
x
n
x
-
x
n
-
1
x
-
x
n
-
2
x
n
+
1
-
x
n
x
n
+
1
-
x
n
-
1
x
n
+
1
-
x
n
-
2
y
n
+
1
.Replacing x
=
s
h
+
x
n
+
2 into (20) produces(21)
P
s
h
+
x
n
+
2
=
s
h
s
h
+
h
s
h
+
2
h
s
h
+
3
h
-
4
h
-
3
h
-
2
h
-
h
y
n
-
2
+
s
h
s
h
+
h
s
h
+
3
h
s
h
+
4
h
-
3
h
-
2
h
-
h
h
y
n
-
1
+
s
h
s
h
+
h
s
h
+
2
h
s
h
+
4
h
-
2
h
-
h
h
2
h
y
n
+
s
h
s
h
+
2
h
s
h
+
3
h
s
h
+
4
h
-
h
h
2
h
3
h
y
n
+
1
+
s
h
+
h
h
s
h
+
2
h
s
h
+
3
h
s
h
+
4
h
h
2
h
3
h
4
h
y
n
+
2
.The resulting polynomial above is differentiated once with respect to s at the point x
=
x
n
+
2. Substituting s
=
0 will give(22)
P
′
x
n
+
2
=
25
12
y
n
+
2
+
1
4
y
n
-
2
-
4
3
y
n
-
1
+
3
y
n
-
4
y
n
+
1
.Considering h
f
n
+
1
,
n
+
2
=
h
P
′
x
n
+
1
,
n
+
2, the corrector formula of DI2BBDF(3) is given by(23)
y
n
+
1
=
2
11
y
n
-
2
-
9
11
y
n
-
1
+
18
11
y
n
+
6
11
h
f
n
+
1
,
y
n
+
2
=
-
3
25
y
n
-
2
+
16
25
y
n
-
1
-
36
25
y
n
+
48
25
y
n
+
1
+
12
25
h
f
n
+
2
.


### 2.3. Derivation of DI2BBDF(4)

The interpolating pointsx
n
-
3
,
…
,
x
n
+
1 are used to obtain the first point y
n
+
1 of DI2BBDF(4):(24)
P
x
=
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
1
x
-
x
n
-
2
x
n
-
3
-
x
n
+
1
x
n
-
3
-
x
n
x
n
-
3
-
x
n
-
1
x
n
-
3
-
x
n
-
2
y
n
-
3
+
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
1
x
-
x
n
-
3
x
n
-
2
-
x
n
+
1
x
n
-
2
-
x
n
x
n
-
2
-
x
n
-
1
x
n
-
2
-
x
n
-
3
y
n
-
2
+
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
2
x
-
x
n
-
3
x
n
-
1
-
x
n
+
1
x
n
-
1
-
x
n
x
n
-
1
-
x
n
-
2
x
n
-
1
-
x
n
-
3
y
n
-
1
+
x
-
x
n
+
1
x
-
x
n
-
1
x
-
x
n
-
2
x
-
x
n
-
3
x
n
-
x
n
+
1
x
n
-
x
n
-
1
x
n
-
x
n
-
2
x
n
-
x
n
-
3
y
n
+
x
-
x
n
x
-
x
n
-
1
x
-
x
n
-
2
x
-
x
n
-
3
x
n
+
1
-
x
n
x
n
+
1
-
x
n
-
1
x
n
+
1
-
x
n
-
2
x
n
+
1
-
x
n
-
3
y
n
+
1
.We define x
=
s
h
+
x
n
+
1 and produce(25)
P
s
h
+
x
n
+
1
=
s
h
s
h
+
h
s
h
+
2
h
s
h
+
3
h
-
4
h
-
3
h
-
2
h
-
h
y
n
-
3
+
s
h
s
h
+
h
s
h
+
2
h
s
h
+
4
h
-
3
h
-
2
h
-
h
h
y
n
-
2
+
s
h
s
h
+
h
s
h
+
3
h
s
h
+
3
h
-
2
h
-
h
h
2
h
y
n
-
1
+
s
h
s
h
+
2
h
s
h
+
3
h
s
h
+
3
h
-
h
h
2
h
3
h
y
n
+
s
h
+
h
s
h
+
2
h
s
h
+
3
h
s
h
+
3
h
h
2
h
3
h
4
h
y
n
+
1
.Differentiating (25) once with respect to s at the point x
=
x
n
+
1 and evaluating s
=
0 will produce(26)
P
′
x
n
+
1
=
25
12
y
n
+
1
+
1
4
y
n
-
3
-
4
3
y
n
-
2
+
3
y
n
-
1
-
4
y
n
.The derivation continues for the second point y
n
+
2 of formula by using x
n
-
3
,
…
,
x
n
+
2 as the interpolating points:(27)
P
x
=
x
-
x
n
+
2
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
1
x
-
x
n
-
2
x
n
-
3
-
x
n
+
2
x
n
-
3
-
x
n
+
1
x
n
-
3
-
x
n
x
n
-
3
-
x
n
-
1
x
n
-
3
-
x
n
-
2
y
n
-
3
+
x
-
x
n
+
2
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
1
x
-
x
n
-
3
x
n
-
2
-
x
n
+
2
x
n
-
2
-
x
n
+
1
x
n
-
2
-
x
n
x
n
-
2
-
x
n
-
1
x
n
-
2
-
x
n
-
3
y
n
-
2
+
x
-
x
n
+
2
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
2
x
-
x
n
-
3
x
n
-
1
-
x
n
+
2
x
n
-
1
-
x
n
+
1
x
n
-
1
-
x
n
x
n
-
1
-
x
n
-
2
x
n
-
1
-
x
n
-
3
y
n
-
1
+
x
-
x
n
+
2
x
-
x
n
+
1
x
-
x
n
-
1
x
-
x
n
-
2
x
-
x
n
-
3
x
n
-
x
n
+
2
x
n
-
x
n
+
1
x
n
-
x
n
-
1
x
n
-
x
n
-
2
x
n
-
x
n
-
3
y
n
+
x
-
x
n
+
2
x
-
x
n
x
-
x
n
-
1
x
-
x
n
-
2
x
-
x
n
-
3
x
n
+
1
-
x
n
+
2
x
n
+
1
-
x
n
x
n
+
1
-
x
n
-
1
x
n
+
1
-
x
n
-
2
x
n
+
1
-
x
n
-
3
y
n
+
1
+
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
1
x
-
x
n
-
2
x
-
x
n
-
3
x
n
+
2
-
x
n
+
1
x
n
+
2
-
x
n
x
n
+
2
-
x
n
-
1
x
n
+
2
-
x
n
-
2
x
n
+
2
-
x
n
-
3
y
n
+
2
.Substituting x
=
h
+
x
n
+
2 into (27) will produce(28)
P
s
h
+
x
n
+
2
=
s
h
s
h
+
h
s
h
+
2
h
s
h
+
3
h
s
h
+
4
h
-
5
h
-
4
h
-
3
h
-
2
h
-
h
y
n
-
3
+
s
h
s
h
+
h
s
h
+
2
h
s
h
+
3
h
s
h
+
5
h
-
4
h
-
3
h
-
2
h
-
h
h
y
n
-
2
+
s
h
s
h
+
h
s
h
+
2
h
s
h
+
4
h
s
h
+
5
h
-
3
h
-
2
h
-
h
h
2
h
y
n
-
1
+
s
h
s
h
+
h
s
h
+
3
h
s
h
+
4
h
s
h
+
5
h
-
2
h
-
h
h
2
h
3
h
y
n
+
s
h
s
h
+
2
h
s
h
+
3
h
s
h
+
4
h
s
h
+
5
h
-
h
h
2
h
3
h
4
h
y
n
+
1
+
s
h
+
h
s
h
+
2
h
s
h
+
3
h
s
h
+
4
h
s
h
+
5
h
h
2
h
3
h
4
h
5
h
y
n
+
2
.Equation (28) is differentiated once with respect to s followed by substituting s
=
0; we have(29)
P
′
x
n
+
2
=
137
60
y
n
+
2
-
1
5
y
n
-
3
+
5
4
y
n
-
2
-
10
3
y
n
-
1
+
5
y
n
-
5
y
n
+
1
.Therefore, the corrector formula of DI2BBDF(4) is obtained as follows:(30)
y
n
+
1
=
-
3
25
y
n
-
3
+
16
25
y
n
-
2
-
36
25
y
n
-
1
y
n
+
1
=
+
48
25
y
n
+
12
25
h
f
n
+
1
,
y
n
+
2
=
12
137
y
n
-
3
-
75
137
y
n
-
2
+
200
137
y
n
-
1
y
n
+
2
=
-
300
137
y
n
+
300
137
y
n
+
1
+
60
137
h
f
n
+
2
.


## 2.1. Derivation of DI2BBDF(2)

DI2BBDF(2) will compute two approximated solutionsy
n
+
1 and y
n
+
2 simultaneously in each block using two back values. This formula is derived using interpolating points x
n
-
1
,
…
,
x
n
+
1 to obtain the first formula y
n
+
1 of DI2BBDF(2): (10)
P
x
=
x
-
x
n
x
-
x
n
+
1
x
n
-
1
-
x
n
x
n
-
1
-
x
n
+
1
y
n
-
1
+
x
-
x
n
-
1
x
-
x
n
+
1
x
n
-
x
n
-
1
x
n
-
x
n
+
1
y
n
+
x
-
x
n
-
1
x
-
x
n
x
n
+
1
-
x
n
-
1
x
n
+
1
-
x
n
y
n
+
1
.Replacing x
=
s
h
+
x
n
+
1 into (10) yields(11)
P
x
n
+
1
+
s
h
=
s
h
+
h
s
h
-
h
-
2
h
y
n
-
1
+
2
h
+
s
h
s
h
h
-
h
y
n
+
2
h
+
s
h
h
+
s
h
2
h
h
y
n
+
1
.Equation (11) is differentiated once with respect to s at the point x
=
x
n
+
1. Evaluating s
=
0 gives(12)
P
′
x
n
+
1
=
1
2
y
n
-
1
-
2
y
n
+
3
2
y
n
+
1
.The same technique is applied for the second point y
n
+
2 of DI2BBDF(2). This formula is derived using x
n
-
1
,
…
,
x
n
+
2 as the interpolating points and produces(13)
P
x
=
x
-
x
n
x
-
x
n
+
1
x
-
x
n
+
2
x
n
-
1
-
x
n
x
n
-
1
-
x
n
+
1
x
n
-
1
-
x
n
+
2
y
n
-
1
+
x
-
x
n
-
1
x
-
x
n
+
1
x
-
x
n
+
2
x
n
-
x
n
-
1
x
n
-
x
n
+
1
x
n
-
x
n
+
2
y
n
+
x
-
x
n
-
1
x
-
x
n
x
-
x
n
+
2
x
n
+
1
-
x
n
-
1
x
n
+
1
-
x
n
x
n
+
1
-
x
n
+
2
y
n
+
1
+
x
-
x
n
-
1
x
-
x
n
x
-
x
n
+
1
x
n
+
2
-
x
n
-
1
x
n
+
2
-
x
n
x
n
+
2
-
x
n
+
1
y
n
+
2
.Substituting x
=
s
h
+
x
n
+
2 into (13) yields(14)
P
x
n
+
2
+
s
h
=
s
h
+
2
h
s
h
+
h
s
h
-
h
-
2
h
-
3
h
y
n
-
1
+
s
h
+
3
h
s
h
+
h
s
h
h
-
h
-
2
h
y
n
+
s
h
+
3
h
s
h
+
2
h
s
h
2
h
h
-
h
y
n
+
1
+
(
s
h
+
3
h
)
(
s
h
+
2
h
)
s
h
+
h
(
3
h
)
(
2
h
)
h
y
n
+
2
.Differentiating (14) with respect to s at the point x
=
x
n
+
2 gives(15)
P
′
x
n
+
2
=
-
1
3
y
n
-
1
+
3
2
y
n
-
3
y
n
+
1
+
11
6
y
n
+
2
.Considering h
f
n
+
1
,
n
+
2
=
P
′
x
n
+
1
,
n
+
2, the corrector formula of DI2BBDF(2) is given as follows:(16)
y
n
+
1
=
-
1
3
y
n
-
1
+
4
3
y
n
+
2
3
h
f
n
+
1
,
y
n
+
2
=
2
11
y
n
-
1
-
9
11
y
n
+
18
11
y
n
+
1
+
6
11
h
f
n
+
2
.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 pointy
n
+
1 of DI2BBDF(3) is derived using interpolating points x
n
-
2
,
…
,
x
n
+
1, and we have(17)
P
x
=
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
1
x
n
-
2
-
x
n
+
1
x
n
-
2
-
x
n
x
n
-
2
-
x
n
-
1
y
n
-
2
+
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
2
x
n
-
1
-
x
n
+
1
x
n
-
1
-
x
n
x
n
-
1
-
x
n
-
2
y
n
-
1
+
x
-
x
n
+
1
x
-
x
n
-
1
x
-
x
n
-
2
x
n
-
x
n
+
1
x
n
-
x
n
-
1
x
n
-
x
n
-
2
y
n
+
x
-
x
n
x
-
x
n
-
1
x
-
x
n
-
2
x
n
+
1
-
x
n
x
n
+
1
-
x
n
-
1
x
n
+
1
-
x
n
-
2
y
n
+
1
.Substituting x
=
s
h
+
x
n
+
1 into (17) gives(18)
P
s
h
+
x
n
+
1
=
s
h
s
h
+
h
s
h
+
2
h
-
3
h
-
2
h
-
h
y
n
-
2
+
s
h
s
h
+
h
s
h
+
3
h
-
2
h
-
h
h
y
n
-
1
+
s
h
s
h
+
2
h
s
h
+
3
h
-
h
h
2
h
y
n
+
s
h
+
h
s
h
+
2
h
s
h
+
3
h
h
2
h
3
h
y
n
+
1
.Equation (18) is differentiated once with respect to s at the point x
=
x
n
+
1. Substituting s
=
0 will obtain(19)
P
′
x
n
+
1
=
11
6
y
n
+
1
-
1
3
y
n
-
2
+
3
2
y
n
-
1
-
3
y
n
.The derivation process continues for second point y
n
+
2 of DI2BBDF(3) using x
n
-
3
,
…
,
x
n
+
1 as the interpolating points. We have(20)
P
x
=
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
1
x
n
-
2
-
x
n
+
1
x
n
-
2
-
x
n
x
n
-
2
-
x
n
-
1
y
n
-
2
+
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
2
x
n
-
1
-
x
n
+
1
x
n
-
1
-
x
n
x
n
-
1
-
x
n
-
2
y
n
-
1
+
x
-
x
n
+
1
x
-
x
n
-
1
x
-
x
n
-
2
x
n
-
x
n
+
1
x
n
-
x
n
-
1
x
n
-
x
n
-
2
y
n
+
x
-
x
n
x
-
x
n
-
1
x
-
x
n
-
2
x
n
+
1
-
x
n
x
n
+
1
-
x
n
-
1
x
n
+
1
-
x
n
-
2
y
n
+
1
.Replacing x
=
s
h
+
x
n
+
2 into (20) produces(21)
P
s
h
+
x
n
+
2
=
s
h
s
h
+
h
s
h
+
2
h
s
h
+
3
h
-
4
h
-
3
h
-
2
h
-
h
y
n
-
2
+
s
h
s
h
+
h
s
h
+
3
h
s
h
+
4
h
-
3
h
-
2
h
-
h
h
y
n
-
1
+
s
h
s
h
+
h
s
h
+
2
h
s
h
+
4
h
-
2
h
-
h
h
2
h
y
n
+
s
h
s
h
+
2
h
s
h
+
3
h
s
h
+
4
h
-
h
h
2
h
3
h
y
n
+
1
+
s
h
+
h
h
s
h
+
2
h
s
h
+
3
h
s
h
+
4
h
h
2
h
3
h
4
h
y
n
+
2
.The resulting polynomial above is differentiated once with respect to s at the point x
=
x
n
+
2. Substituting s
=
0 will give(22)
P
′
x
n
+
2
=
25
12
y
n
+
2
+
1
4
y
n
-
2
-
4
3
y
n
-
1
+
3
y
n
-
4
y
n
+
1
.Considering h
f
n
+
1
,
n
+
2
=
h
P
′
x
n
+
1
,
n
+
2, the corrector formula of DI2BBDF(3) is given by(23)
y
n
+
1
=
2
11
y
n
-
2
-
9
11
y
n
-
1
+
18
11
y
n
+
6
11
h
f
n
+
1
,
y
n
+
2
=
-
3
25
y
n
-
2
+
16
25
y
n
-
1
-
36
25
y
n
+
48
25
y
n
+
1
+
12
25
h
f
n
+
2
.


## 2.3. Derivation of DI2BBDF(4)

The interpolating pointsx
n
-
3
,
…
,
x
n
+
1 are used to obtain the first point y
n
+
1 of DI2BBDF(4):(24)
P
x
=
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
1
x
-
x
n
-
2
x
n
-
3
-
x
n
+
1
x
n
-
3
-
x
n
x
n
-
3
-
x
n
-
1
x
n
-
3
-
x
n
-
2
y
n
-
3
+
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
1
x
-
x
n
-
3
x
n
-
2
-
x
n
+
1
x
n
-
2
-
x
n
x
n
-
2
-
x
n
-
1
x
n
-
2
-
x
n
-
3
y
n
-
2
+
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
2
x
-
x
n
-
3
x
n
-
1
-
x
n
+
1
x
n
-
1
-
x
n
x
n
-
1
-
x
n
-
2
x
n
-
1
-
x
n
-
3
y
n
-
1
+
x
-
x
n
+
1
x
-
x
n
-
1
x
-
x
n
-
2
x
-
x
n
-
3
x
n
-
x
n
+
1
x
n
-
x
n
-
1
x
n
-
x
n
-
2
x
n
-
x
n
-
3
y
n
+
x
-
x
n
x
-
x
n
-
1
x
-
x
n
-
2
x
-
x
n
-
3
x
n
+
1
-
x
n
x
n
+
1
-
x
n
-
1
x
n
+
1
-
x
n
-
2
x
n
+
1
-
x
n
-
3
y
n
+
1
.We define x
=
s
h
+
x
n
+
1 and produce(25)
P
s
h
+
x
n
+
1
=
s
h
s
h
+
h
s
h
+
2
h
s
h
+
3
h
-
4
h
-
3
h
-
2
h
-
h
y
n
-
3
+
s
h
s
h
+
h
s
h
+
2
h
s
h
+
4
h
-
3
h
-
2
h
-
h
h
y
n
-
2
+
s
h
s
h
+
h
s
h
+
3
h
s
h
+
3
h
-
2
h
-
h
h
2
h
y
n
-
1
+
s
h
s
h
+
2
h
s
h
+
3
h
s
h
+
3
h
-
h
h
2
h
3
h
y
n
+
s
h
+
h
s
h
+
2
h
s
h
+
3
h
s
h
+
3
h
h
2
h
3
h
4
h
y
n
+
1
.Differentiating (25) once with respect to s at the point x
=
x
n
+
1 and evaluating s
=
0 will produce(26)
P
′
x
n
+
1
=
25
12
y
n
+
1
+
1
4
y
n
-
3
-
4
3
y
n
-
2
+
3
y
n
-
1
-
4
y
n
.The derivation continues for the second point y
n
+
2 of formula by using x
n
-
3
,
…
,
x
n
+
2 as the interpolating points:(27)
P
x
=
x
-
x
n
+
2
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
1
x
-
x
n
-
2
x
n
-
3
-
x
n
+
2
x
n
-
3
-
x
n
+
1
x
n
-
3
-
x
n
x
n
-
3
-
x
n
-
1
x
n
-
3
-
x
n
-
2
y
n
-
3
+
x
-
x
n
+
2
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
1
x
-
x
n
-
3
x
n
-
2
-
x
n
+
2
x
n
-
2
-
x
n
+
1
x
n
-
2
-
x
n
x
n
-
2
-
x
n
-
1
x
n
-
2
-
x
n
-
3
y
n
-
2
+
x
-
x
n
+
2
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
2
x
-
x
n
-
3
x
n
-
1
-
x
n
+
2
x
n
-
1
-
x
n
+
1
x
n
-
1
-
x
n
x
n
-
1
-
x
n
-
2
x
n
-
1
-
x
n
-
3
y
n
-
1
+
x
-
x
n
+
2
x
-
x
n
+
1
x
-
x
n
-
1
x
-
x
n
-
2
x
-
x
n
-
3
x
n
-
x
n
+
2
x
n
-
x
n
+
1
x
n
-
x
n
-
1
x
n
-
x
n
-
2
x
n
-
x
n
-
3
y
n
+
x
-
x
n
+
2
x
-
x
n
x
-
x
n
-
1
x
-
x
n
-
2
x
-
x
n
-
3
x
n
+
1
-
x
n
+
2
x
n
+
1
-
x
n
x
n
+
1
-
x
n
-
1
x
n
+
1
-
x
n
-
2
x
n
+
1
-
x
n
-
3
y
n
+
1
+
x
-
x
n
+
1
x
-
x
n
x
-
x
n
-
1
x
-
x
n
-
2
x
-
x
n
-
3
x
n
+
2
-
x
n
+
1
x
n
+
2
-
x
n
x
n
+
2
-
x
n
-
1
x
n
+
2
-
x
n
-
2
x
n
+
2
-
x
n
-
3
y
n
+
2
.Substituting x
=
h
+
x
n
+
2 into (27) will produce(28)
P
s
h
+
x
n
+
2
=
s
h
s
h
+
h
s
h
+
2
h
s
h
+
3
h
s
h
+
4
h
-
5
h
-
4
h
-
3
h
-
2
h
-
h
y
n
-
3
+
s
h
s
h
+
h
s
h
+
2
h
s
h
+
3
h
s
h
+
5
h
-
4
h
-
3
h
-
2
h
-
h
h
y
n
-
2
+
s
h
s
h
+
h
s
h
+
2
h
s
h
+
4
h
s
h
+
5
h
-
3
h
-
2
h
-
h
h
2
h
y
n
-
1
+
s
h
s
h
+
h
s
h
+
3
h
s
h
+
4
h
s
h
+
5
h
-
2
h
-
h
h
2
h
3
h
y
n
+
s
h
s
h
+
2
h
s
h
+
3
h
s
h
+
4
h
s
h
+
5
h
-
h
h
2
h
3
h
4
h
y
n
+
1
+
s
h
+
h
s
h
+
2
h
s
h
+
3
h
s
h
+
4
h
s
h
+
5
h
h
2
h
3
h
4
h
5
h
y
n
+
2
.Equation (28) is differentiated once with respect to s followed by substituting s
=
0; we have(29)
P
′
x
n
+
2
=
137
60
y
n
+
2
-
1
5
y
n
-
3
+
5
4
y
n
-
2
-
10
3
y
n
-
1
+
5
y
n
-
5
y
n
+
1
.Therefore, the corrector formula of DI2BBDF(4) is obtained as follows:(30)
y
n
+
1
=
-
3
25
y
n
-
3
+
16
25
y
n
-
2
-
36
25
y
n
-
1
y
n
+
1
=
+
48
25
y
n
+
12
25
h
f
n
+
1
,
y
n
+
2
=
12
137
y
n
-
3
-
75
137
y
n
-
2
+
200
137
y
n
-
1
y
n
+
2
=
-
300
137
y
n
+
300
137
y
n
+
1
+
60
137
h
f
n
+
2
.


## 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-planeRe
h
λ
<
0.We consider the simplest test equation(31)
f
=
y
′
=
λ
y
,where the eigenvalues λ
i, i
=
1
,
2
,
…
,
s, satisfy λ
i
<
0 to analyze the stiff stability of the method. Substituting (31) into (16) will produce the following matrices form:(32)
1
0
-
18
11
1
y
n
+
1
y
n
+
2
=
-
1
3
4
3
2
11
-
9
11
y
n
-
1
y
n
+
h
2
3
0
0
6
11
f
n
+
1
f
n
+
2
,where(33)
A
=
1
-
2
3
λ
h
0
-
18
11
1
-
6
11
λ
h
,
B
=
-
1
3
4
3
2
11
-
9
11
,
C
=
0
0
0
0
.Let h
^
=
h
λ and evaluating the determinant of |
A
t
-
B
| from (32), the stability polynomial ρ
t
,
h
^ of DI2BBDF(2) is obtained as follows:(34)
ρ
t
,
h
^
=
1
33
-
34
t
33
-
8
h
^
t
11
+
t
2
-
40
h
^
t
2
33
+
4
h
^
2
t
2
11
.The similar approach is applied to obtain the stability polynomial of DI2BBDF(3) and DI2BBDF(4). Substituting (31) into (23), we have(35)
1
0
-
48
25
1
y
n
+
1
y
n
+
2
=
0
2
11
0
-
3
25
y
n
-
3
y
n
-
2
+
-
9
11
18
11
16
25
-
36
25
y
n
-
1
y
n
+
h
6
11
0
0
12
25
f
n
+
1
f
n
+
2
,where(36)
A
=
1
-
6
11
λ
h
0
-
48
25
1
-
12
25
λ
h
,
B
=
-
9
11
18
11
16
25
-
36
25
,
C
=
0
2
11
0
3
25
.Substituting (31) into (30) will give(37)
1
0
-
300
137
1
y
n
+
1
y
n
+
2
=
-
3
25
16
25
12
137
-
75
137
y
n
-
3
y
n
-
2
+
-
36
25
48
25
200
137
-
300
137
y
n
-
1
y
n
+
h
12
25
0
0
60
137
f
n
+
1
f
n
+
2
,where(38)
A
=
1
-
12
25
λ
h
0
-
300
137
1
-
60
137
λ
h
,
B
=
-
36
25
48
25
200
137
-
300
137
,
C
=
-
3
25
16
25
12
137
-
75
137
.We compute the determinant ofA
t
2
-
B
t
-
C from (35) and (37) to obtain the stability polynomials ρ
t
,
h
^ of DI2BBDF(3) and DI2BBDF(4), respectively. Consider(39)
ρ
t
,
h
^
=
36
275
-
86
t
275
-
243
t
2
275
-
324
h
^
t
2
275
-
162
t
3
275
+
36
h
^
t
3
275
+
t
4
-
282
h
^
t
4
275
+
72
h
^
2
t
4
275
,
(40)
ρ
t
,
h
^
=
48
137
-
176
t
3425
-
387
t
2
685
-
1152
h
^
t
2
685
-
2514
t
3
3425
-
216
h
^
t
3
685
+
t
4
-
3144
h
^
t
4
3425
+
144
h
^
2
t
4
685
.Next, the boundary of the stability region will be determined by substitutingt
=
e
i
θ into (34), (39), and (40). The graphs of stability region for all formulas are given in Figures 1, 2, and 3.Figure 1
Graph of stability region for DI2BBDF(2).Figure 2
Graph of stability region for DI2BBDF(3).Figure 3
Graph of stability region for DI2BBDF(4).In Figures1, 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).Figure 4
Graph of stability region for DI2BBDF(2), DI2BBDF(3), and DI2BBDF(4).Clearly, the unstable region becomes larger when the order of the method increases. From  Definition4, 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:(41)
y
n
+
1
=
α
1
y
n
+
2
+
β
1
h
f
n
+
1
+
μ
1
,
y
n
+
2
=
α
2
y
n
+
1
+
β
2
h
f
n
+
2
+
μ
2
,where μ
1 and μ
2 are the back values. Equation (41) is transformed into matrix form as follows:(42)
1
0
0
1
y
n
+
1
y
n
+
2
=
0
α
1
α
2
0
y
n
+
1
y
n
+
2
+
h
β
1
0
0
β
2
f
n
+
1
f
n
+
2
+
μ
1
μ
2
,where (43)
I
=
1
0
0
1
,
Y
=
y
n
+
1
y
n
+
2
,
B
=
β
1
0
0
β
2
,
F
=
f
n
+
1
f
n
+
2
,
μ
=
μ
1
μ
2
.Implementing Newton’s method to (42) produces(44)
F
^
=
I
-
A
Y
-
h
B
F
+
μ
.Iteration for (44) is given by(45)
y
n
+
1
,
n
+
2
i
+
1
-
y
n
+
1
,
n
+
2
i
=
F
y
n
+
1
,
n
+
2
i
·
F
′
y
n
+
1
,
n
+
2
i
-
1
,where y
n
+
1
(
i
+
1
) denotes the (
i
+
1
)th iteration. Equation (45) can be rearranged into the following equation:(46)
F
′
y
n
+
1
,
n
+
2
i
y
n
+
1
,
n
+
2
i
+
1
-
y
n
+
1
,
n
+
2
i
=
F
y
n
+
1
,
n
+
2
i
.Replacing (44) into (46) yields(47)
I
-
A
Y
-
h
B
∂
f
∂
y
y
n
+
1
,
n
+
2
i
y
n
+
1
,
n
+
2
i
+
1
-
y
n
+
1
,
n
+
2
i
=
I
-
A
y
n
+
1
,
n
+
2
i
-
h
B
y
n
+
1
,
n
+
2
i
+
μ
,where ∂
f
/
∂
y
y
n
+
1
,
n
+
2
i is the Jacobian matrix of f with respect to y.


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

Formula (16) is written in the form of (41):(48)
F
1
=
y
n
+
1
-
2
3
h
f
n
+
1
+
μ
1
,
F
2
=
y
n
+
2
-
18
11
y
n
+
1
-
6
11
h
f
n
+
2
+
μ
2
.Applying (48) into (47) will yield matrix form as follows:(49)
1
-
2
3
h
∂
f
n
+
1
∂
y
n
+
1
0
0
1
-
6
11
h
∂
f
n
+
2
∂
y
n
+
2
e
n
+
1
(
i
+
1
)
e
n
+
2
(
i
+
1
)
=
-
1
0
18
11
-
1
y
n
+
1
i
y
n
+
2
i
+
h
2
3
0
0
6
11
f
n
+
1
i
f
n
+
2
i
+
μ
1
μ
2
,with e
n
+
1
,
n
+
2
i
+
1
=
y
n
+
1
,
n
+
2
i
+
1
-
y
n
+
1
,
n
+
2
i being the increment.


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

Rearranging formula (23) in the form of (41), we have(50)
F
1
=
y
n
+
1
-
6
11
h
f
n
+
1
+
μ
1
,
F
2
=
y
n
+
2
-
48
25
y
n
+
1
-
12
25
h
f
n
+
2
+
μ
2
.Replacing (50) into (47), this will yield matrix form as follows:(51)
1
-
6
11
h
∂
f
n
+
1
∂
y
n
+
1
0
0
1
-
12
25
h
∂
f
n
+
2
∂
y
n
+
2
e
n
+
1
(
i
+
1
)
e
n
+
2
(
i
+
1
)
=
-
1
0
48
25
-
1
y
n
+
1
i
y
n
+
2
i
+
h
6
11
0
0
12
25
f
n
+
1
i
f
n
+
2
i
+
μ
1
μ
2
,with e
n
+
1
,
n
+
2
i
+
1
=
y
n
+
1
,
n
+
2
i
+
1
-
y
n
+
1
,
n
+
2
i being the increment.


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

We rewrite formula (30) in the form of (41) and obtain(52)
F
1
=
y
n
+
1
-
12
25
h
f
n
+
1
+
μ
1
,
F
2
=
y
n
+
2
-
300
137
y
n
+
1
-
60
137
h
f
n
+
2
+
μ
2
.Substituting (52) into (47) will produce(53)
1
-
12
25
h
∂
f
n
+
1
∂
y
n
+
1
0
0
1
-
60
137
h
∂
f
n
+
2
∂
y
n
+
2
e
n
+
1
(
i
+
1
)
e
n
+
2
(
i
+
1
)
=
-
1
0
300
137
-
1
y
n
+
1
i
y
n
+
2
i
+
h
12
25
0
0
60
137
f
n
+
1
i
f
n
+
2
i
+
μ
1
μ
2
,with e
n
+
1
,
n
+
2
i
+
1
=
y
n
+
1
,
n
+
2
i
+
1
-
y
n
+
1
,
n
+
2
i 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 functionf, given x and y. The approximation calculations for y
n
+
1 and y
n
+
2 in PECE are as follows.(1)
P (Predict):y
n
+
1
,
n
+
2
p.
(2)
E (Evaluate):y
n
+
1
,
n
+
2
′
=
f
x
n
+
1
,
n
+
2
,
y
n
+
1
,
n
+
2
p.
(3)
C (Correct):y
n
+
1
,
n
+
2
c.
(4)
E (Evaluate):y
n
+
1
,
n
+
2
′
=
f
x
n
+
1
,
n
+
2
,
y
n
+
1
,
n
+
2
c.To approximate the solution to y
n
+
1
,
n
+
2
c, we apply the two-stage Newton type iteration. The iteration process for DI2BBDF(
2
) is done as follows.(1)
Compute the values fore
n
+
1
,
n
+
2
i
+
1
=
A
-
1
B, where(54)
A
=
1
-
2
3
h
∂
f
n
+
1
∂
y
n
+
1
0
0
1
-
6
11
h
∂
f
n
+
2
∂
y
n
+
2
,
B
=
-
y
n
+
1
i
+
2
3
h
f
n
+
1
i
+
μ
1
0
18
11
y
n
+
1
i
-
y
n
+
2
i
+
6
11
h
f
n
+
2
i
+
μ
2
.
(2)
Calculate the corrected value fory
n
+
1
,
n
+
2
i
+
1 with the value e
n
+
1
,
n
+
2
i
+
1 from Step (
1
).
(3)
Solvee
n
+
1
,
n
+
2
i
+
1
=
A
-
1
B for second stage iteration.
(4)
The final values fory
n
+
1
,
n
+
2
i
+
1 are obtained from the second stage iteration of e
n
+
1
,
n
+
2
i
+
1.The similar iteration process is applied for DI2BBDF(3) and DI2BBDF(4).


## 4.1. Newton Iteration of DI2BBDF(2)

Formula (16) is written in the form of (41):(48)
F
1
=
y
n
+
1
-
2
3
h
f
n
+
1
+
μ
1
,
F
2
=
y
n
+
2
-
18
11
y
n
+
1
-
6
11
h
f
n
+
2
+
μ
2
.Applying (48) into (47) will yield matrix form as follows:(49)
1
-
2
3
h
∂
f
n
+
1
∂
y
n
+
1
0
0
1
-
6
11
h
∂
f
n
+
2
∂
y
n
+
2
e
n
+
1
(
i
+
1
)
e
n
+
2
(
i
+
1
)
=
-
1
0
18
11
-
1
y
n
+
1
i
y
n
+
2
i
+
h
2
3
0
0
6
11
f
n
+
1
i
f
n
+
2
i
+
μ
1
μ
2
,with e
n
+
1
,
n
+
2
i
+
1
=
y
n
+
1
,
n
+
2
i
+
1
-
y
n
+
1
,
n
+
2
i being the increment.


## 4.2. Newton Iteration of DI2BBDF(3)

Rearranging formula (23) in the form of (41), we have(50)
F
1
=
y
n
+
1
-
6
11
h
f
n
+
1
+
μ
1
,
F
2
=
y
n
+
2
-
48
25
y
n
+
1
-
12
25
h
f
n
+
2
+
μ
2
.Replacing (50) into (47), this will yield matrix form as follows:(51)
1
-
6
11
h
∂
f
n
+
1
∂
y
n
+
1
0
0
1
-
12
25
h
∂
f
n
+
2
∂
y
n
+
2
e
n
+
1
(
i
+
1
)
e
n
+
2
(
i
+
1
)
=
-
1
0
48
25
-
1
y
n
+
1
i
y
n
+
2
i
+
h
6
11
0
0
12
25
f
n
+
1
i
f
n
+
2
i
+
μ
1
μ
2
,with e
n
+
1
,
n
+
2
i
+
1
=
y
n
+
1
,
n
+
2
i
+
1
-
y
n
+
1
,
n
+
2
i being the increment.


## 4.3. Newton Iteration of DI2BBDF(4)

We rewrite formula (30) in the form of (41) and obtain(52)
F
1
=
y
n
+
1
-
12
25
h
f
n
+
1
+
μ
1
,
F
2
=
y
n
+
2
-
300
137
y
n
+
1
-
60
137
h
f
n
+
2
+
μ
2
.Substituting (52) into (47) will produce(53)
1
-
12
25
h
∂
f
n
+
1
∂
y
n
+
1
0
0
1
-
60
137
h
∂
f
n
+
2
∂
y
n
+
2
e
n
+
1
(
i
+
1
)
e
n
+
2
(
i
+
1
)
=
-
1
0
300
137
-
1
y
n
+
1
i
y
n
+
2
i
+
h
12
25
0
0
60
137
f
n
+
1
i
f
n
+
2
i
+
μ
1
μ
2
,with e
n
+
1
,
n
+
2
i
+
1
=
y
n
+
1
,
n
+
2
i
+
1
-
y
n
+
1
,
n
+
2
i 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 functionf, given x and y. The approximation calculations for y
n
+
1 and y
n
+
2 in PECE are as follows.(1)
P (Predict):y
n
+
1
,
n
+
2
p.
(2)
E (Evaluate):y
n
+
1
,
n
+
2
′
=
f
x
n
+
1
,
n
+
2
,
y
n
+
1
,
n
+
2
p.
(3)
C (Correct):y
n
+
1
,
n
+
2
c.
(4)
E (Evaluate):y
n
+
1
,
n
+
2
′
=
f
x
n
+
1
,
n
+
2
,
y
n
+
1
,
n
+
2
c.To approximate the solution to y
n
+
1
,
n
+
2
c, we apply the two-stage Newton type iteration. The iteration process for DI2BBDF(
2
) is done as follows.(1)
Compute the values fore
n
+
1
,
n
+
2
i
+
1
=
A
-
1
B, where(54)
A
=
1
-
2
3
h
∂
f
n
+
1
∂
y
n
+
1
0
0
1
-
6
11
h
∂
f
n
+
2
∂
y
n
+
2
,
B
=
-
y
n
+
1
i
+
2
3
h
f
n
+
1
i
+
μ
1
0
18
11
y
n
+
1
i
-
y
n
+
2
i
+
6
11
h
f
n
+
2
i
+
μ
2
.
(2)
Calculate the corrected value fory
n
+
1
,
n
+
2
i
+
1 with the value e
n
+
1
,
n
+
2
i
+
1 from Step (
1
).
(3)
Solvee
n
+
1
,
n
+
2
i
+
1
=
A
-
1
B for second stage iteration.
(4)
The final values fory
n
+
1
,
n
+
2
i
+
1 are obtained from the second stage iteration of e
n
+
1
,
n
+
2
i
+
1.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 has(55)
y
′
=
-
10
y
+
10
.Exact solution is(56)
y
x
=
1
+
e
-
10
x
.Initial condition is(57)
y
0
=
2
,
0
≤
x
≤
10(source: Ibrahim [8]).Problem 2 (linear). One has(58)
y
1
′
=
-
5
y
1
-
4
y
2
,
y
2
′
=
-
4
y
1
-
5
y
2
.Exact solutions are(59)
y
1
x
=
1
2
e
-
9
x
+
1
2
e
-
x
,
y
2
x
=
1
2
e
-
9
x
-
1
2
e
-
x
.Initial conditions are(60)
y
1
0
=
1
,
y
2
0
=
0
,
0
≤
x
≤
20(source: Musa et al. [13]).Problem 3 (linear). One has(61)
y
′
=
-
10
x
y
.Exact solution is(62)
y
x
=
e
-
5
x
2
.Initial condition is(63)
y
0
=
1
,
0
≤
x
≤
10(source: Ibrahim et al. [10]).Problem 4 (linear). One has(64)
y
1
′
=
9
y
1
+
24
y
2
+
5
cos
⁡
⁡
x
-
1
3
sin
⁡
x
,
y
2
′
=
-
24
y
1
-
51
y
2
-
9
cos
⁡
⁡
x
+
1
3
sin
⁡
x
.Exact solutions are(65)
y
1
x
=
2
e
-
3
x
-
e
-
39
x
+
1
3
cos
⁡
⁡
x
,
y
2
x
=
-
e
-
3
x
+
2
e
-
39
x
-
1
3
cos
⁡
⁡
x
.Initial conditions are(66)
y
1
0
=
4
3
,
y
2
0
=
2
3
,
0
≤
x
≤
10(source: Musa et al. [14]).Problem 5 (nonlinear). One has(67)
y
′
=
-
y
3
2
.Exact solution is(68)
y
x
=
1
1
+
x
.Initial condition is(69)
y
0
=
1
,
0
≤
x
≤
4(source: Musa et al. [14]).Problem 6 (nonlinear). One has(70)
y
′
=
y
(
1
-
y
)
2
y
-
1
.Exact solution is(71)
y
x
=
1
2
+
1
4
-
5
36
e
-
x
.Initial condition is(72)
y
0
=
5
6
,
0
≤
x
≤
1(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 consider10
-
2, 10
-
4, and 10
-
6 as the step size, h. 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(h) are illustrated as shown in Figures 5–10. The following notations are used in Tables 1–6:h: 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.Table 1
The accuracy for problem 1.
h
Methods
MAXE
TIME
10
-
2
FI2BBDF(3)
5.97499
E
-
2
1.36773
E
-
3
FI2BEBDF(3)
5.67155
E
-
2
1.45525
E
-
3
DI2BBDF(2)
1.10568
E
-
2
4.72498
E
-
3
DI2BBDF(3)
1.66455
E
-
2
4.77409
E
-
3
DI2BBDF(4)
2.16342
E
-
2
4.83298
E
-
3
10
-
4
FI2BBDF(3)
7.34012
E
-
4
1.25797
E
-
1
FI2BEBDF(3)
3.23640
E
-
5
1.36259
E
-
1
DI2BBDF(2)
1.18355
E
-
6
1.67260
E
-
1
DI2BBDF(3)
1.73430
E
-
6
1.69325
E
-
1
DI2BBDF(4)
2.25966
E
-
6
1.72995
E
-
1
10
-
6
FI2BBDF(3)
7.35741
E
-
6
1.25700
E
1
FI2BEBDF(3)
3.50090
E
-
7
1.36230
E
1
DI2BBDF(2)
1.18742
E
-
10
1.49159
E
1
DI2BBDF(3)
1.73634
E
-
10
1.54493
E
1
DI2BBDF(4)
2.26000
E
-
10
1.54620
E
1Table 2
The accuracy for problem 2.
h
Methods
MAXE
TIME
10
-
2
FI2BBDF(3)
2.71597
E
-
2
5.96467
E
-
3
FI2BEBDF(3)
1.49811
E
-
2
6.61025
E
-
3
DI2BBDF(2)
4.56073
E
-
3
1.04480
E
-
2
DI2BBDF(3)
6.86126
E
-
3
1.08421
E
-
2
DI2BBDF(4)
8.93279
E
-
3
1.11089
E
-
2
10
-
4
FI2BBDF(3)
3.40437
E
-
4
5.90430
E
-
1
FI2BEBDF(3)
2.62649
E
-
5
6.54305
E
-
1
DI2BBDF(2)
4.85404
E
-
7
3.73210
E
-
1
DI2BBDF(3)
7.11155
E
-
7
3.86556
E
-
1
DI2BBDF(4)
9.26476
E
-
7
3.87070
E
-
1
10
-
6
FI2BBDF(3)
3.41146
E
-
6
5.90198
E
1
FI2BEBDF(3)
2.44097
E
-
7
6.54106
E
1
DI2BBDF(2)
4.86843
E
-
11
3.48451
E
1
DI2BBDF(3)
7.11903
E
-
11
3.60490
E
1
DI2BBDF(4)
9.26598
E
-
11
3.60649
E
1Table 3
The accuracy for problem 3.
h
Methods
MAXE
TIME
10
-
2
FI3BBDF(3)
3.56692
E
-
2
1.41642
E
-
3
FI3BEBDF(3)
1.24084
E
-
2
3.02167
E
-
3
DI2BBDF(2)
1.17256
E
-
3
4.48489
E
-
3
DI2BBDF(3)
1.72278
E
-
3
4.81200
E
-
3
DI2BBDF(4)
2.25003
E
-
3
4.85301
E
-
3
10
-
4
FI3BBDF(3)
4.35640
E
-
4
1.30103
E
-
1
FI3BEBDF(3)
7.06516
E
-
5
2.08115
E
-
1
DI2BBDF(2)
1.18749
E
-
7
1.60304
E
-
1
DI2BBDF(3)
1.73636
E
-
7
1.77445
E
-
1
DI2BBDF(4)
2.26000
E
-
7
1.78064
E
-
1
10
-
6
FI3BBDF(3)
4.34625
E
-
6
1.29775
E
1
FI3BEBDF(3)
7.03257
E
-
7
2.07609
E
0
DI2BBDF(2)
1.24989
E
-
11
1.50739
E
1
DI2BBDF(3)
1.73636
E
-
11
1.55073
E
1
DI2BBDF(4)
2.25999
E
-
11
1.60216
E
1Table 4
The accuracy for problem 4.
h
Methods
MAXE
TIME
10
-
2
FI3BBDF(3)
6.62694
E
99
7.01775
E
-
3
FI3BEBDF(3)
1.68449
E
-
1
1.37847
E
-
2
DI2BBDF(2)
2.80830
E
-
1
1.06871
E
-
2
DI2BBDF(3)
4.00389
E
-
1
1.24969
E
-
2
DI2BBDF(4)
4.77342
E
-
1
1.26502
E
-
2
10
-
4
FI3BBDF(3)
8.45376
E
-
3
6.76737
E
-
1
FI3BEBDF(3)
6.95725
E
-
3
1.35671
E
+
0
DI2BBDF(2)
3.56255
E
-
5
4.32359
E
-
1
DI2BBDF(3)
5.24241
E
-
5
4.42599
E
-
1
DI2BBDF(4)
6.85105
E
-
5
4.56207
E
-
1
10
-
6
FI3BBDF(3)
8.54545
E
-
5
6.76819
E
1
FI3BEBDF(3)
7.18859
E
-
5
1.35862
E
2
DI2BBDF(2)
3.60130
E
-
9
4.00010
E
1
DI2BBDF(3)
5.26673
E
-
9
4.00192
E
1
DI2BBDF(4)
6.85529
E
-
9
4.10322
E
1Table 5
The accuracy for problem 5.
h
Methods
MAXE
TIME
10
-
2
FI3BBDF(3)
5.25481
E
-
3
5.78083
E
-
4
FI3BEBDF(3)
3.88475
E
-
3
9.16583
E
-
4
DI2BBDF(2)
8.56253
E
-
5
2.29836
E
-
4
DI2BBDF(3)
1.21629
E
-
4
1.89066
E
-
4
DI2BBDF(4)
1.65094
E
-
4
2.30074
E
-
4
10
-
4
FI3BBDF(3)
5.51558
E
-
5
4.83860
E
-
2
FI3BEBDF(3)
4.27506
E
-
5
7.74877
E
-
2
DI2BBDF(2)
8.89803
E
-
9
1.17438
E
-
2
DI2BBDF(3)
1.30124
E
-
8
1.31578
E
-
2
DI2BBDF(4)
1.69456
E
-
8
1.26071
E
-
2
10
-
6
FI3BBDF(3)
5.51803
E
-
7
5.47768
E
+
0
FI3BEBDF(3)
4.27924
E
-
7
7.77643
E
+
0
DI2BBDF(2)
1.26115
E
-
11
1.19318
E
+
0
DI2BBDF(3)
6.27106
E
-
11
1.33892
E
+
0
DI2BBDF(4)
1.04222
E
-
11
1.26146
E
+
0Table 6
The accuracy for problem 6.
h
Methods
MAXE
TIME
10
-
2
FI3BBDF(3)
2.18555
E
-
3
2.83084
E
-
4
FI3BEBDF(3)
1.52719
E
-
3
4.28500
E
-
4
DI2BBDF(2)
3.87041
E
-
5
5.38826
E
-
5
DI2BBDF(3)
5.47958
E
-
5
5.31673
E
-
5
DI2BBDF(4)
7.47696
E
-
5
6.18877
E
-
5
10
-
4
FI3BBDF(3)
2.29812
E
-
5
1.48627
E
-
2
FI3BEBDF(3)
1.78708
E
-
5
2.32184
E
-
2
DI2BBDF(2)
4.01638
E
-
9
3.49712
E
-
3
DI2BBDF(3)
5.87347
E
-
9
3.75295
E
-
3
DI2BBDF(4)
7.64927
E
-
9
3.64804
E
-
3
10
-
6
FI3BBDF(3)
2.29936
E
-
7
1.46660
E
+
0
FI3BEBDF(3)
1.78973
E
-
7
2.28775
E
+
0
DI2BBDF(2)
9.90807
E
-
12
3.38157
E
-
1
DI2BBDF(3)
3.77716
E
-
11
3.69837
E
-
1
DI2BBDF(4)
4.62796
E
-
12
3.56802
E
-
1Figure 5
Graph of log(MAXE) versus log(h) for problem 1.Figure 6
Graph of log(MAXE) versus log(h) for problem 2.Figure 7
Graph of log(MAXE) versus log(h) for problem 3.Figure 8
Graph of log(MAXE) versus log(h) for problem 4.Figure 9
Graph of log(MAXE) versus log(h) for problem 5.Figure 10
Graph of log(MAXE) versus log(h) for problem 6.


## 7. Discussion

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


### 7.1. Maximum Error

From the numerical results in Tables1–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, h, 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 Tables1-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.


## 7.1. Maximum Error

From the numerical results in Tables1–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, h, 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 Tables1-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.



---

*Source: 179231-2015-04-19.xml*