# Block Backward Differentiation Formulas for Fractional Differential Equations

Authors: T. A. Biala; S. N. Jator
Journal: International Journal of Engineering Mathematics
 (2015)

---

## Abstract

This paper concerns the numerical approximation of Fractional Initial Value Problems (FIVPs). This is achieved by constructingk-step continuous BDFs. These continuous schemes are developed via the interpolation and collocation approach and are used to obtain the discrete k-step BDF and (k
-
1) additional methods which are applied as numerical integrators in a block-by-block mode for the integration of FIVP. The properties of the methods are established and regions of absolute stability of the methods are plotted in the complex plane. Numerical tests including large systems arising form the semidiscretization of one-dimensional fractional Burger’s equation show that the methods are highly accurate and efficient.


---

## Body

## 1. Introduction

In what follows, we consider the FIVP of the following form:(1)
D
x
0
α
c
y
x
=
f
x
,
y
x
,
y
x
0
=
y
0
,where 0
<
α
<
1 is the fractional order and D
c
x
0
α (in the sequel we will simply use D
α) denotes the Caputo α derivative operator which is defined as(2)
D
α
y
x
=
1
Γ
1
-
α
∫
x
0
x
x
-
s
-
α
y
′
s
d
s
.The mathematical modeling of several physical phenomena results in fractional differential equations of form (1) and it plays an important role in various branches of science and engineering. Applications of FDEs are found in chemistry, electronics, circuit theory, seismology, signal processing, control theory, and so on. Also, these FDEs serve as a generalization of their corresponding ordinary differential equations (ODEs). For a brief history and introduction to fractional calculus, we refer the reader to [1–3].We have adopted the Caputo’s definition of derivatives of noninteger order (which is a modification of the Riemann-Liouville definition) since it can be coupled with initial conditions having a clear physical meaning. The existence and the uniqueness of the solution of (1) have been given in Diethelm and Ford [4].The development of well suited methods for numerically approximating FDEs has received great attention over the past few decades. This is due to the occurence of FDEs in several models. Several methods have been proposed and analyzed for the numerical approximation of this important class of problems (see Lubich [5–7], Garrappa [8], Galeone and Garrappa [9–11], and the references therein). These authors have independently developed Fractional Linear Multistep Methods (FLMMs) using convolution quadratures. Lubich [6] proposed formulas of the following form:(3)
y
n
=
f
t
n
+
h
α
∑
j
=
0
n
ω
n
-
j
α
g
t
j
,
y
j
+
h
α
∑
j
=
0
m
ω
n
j
g
t
j
,
y
j
,
n
h
∈
I
,where ω
n
α and ω
n
j are the convolution and starting quadrature weights, respectively, and are independent of the stepsize h.One major difficulty in the FLMMs (3) is in evaluating the convolution weights ω
n
α. Most of the methods rely on the J.C.P Miller formula for the computation of these weights. In order to avoid this major drawback, we give a different approach in the construction of the FLMMs. This approach is based on interpolation and collocation as was discussed by Onumanyi et al. [12].The main aim of this paper is to present and investigate a class of fractional BDF methods which generalize the BDF methods for ODEs. The fractional BDF methods are developed using the interpolation and collocation approach and the regions of absolute stability of the methods are plotted via the boundary locus method.The paper is organized as follows: in Section2, we discuss the development of the fractional BDF methods. Section 3 details the convergence of the methods while, in Section 4, we give the stability properties and implementation of the methods. In Section 5, we give five numerical examples to elucidate our theoretical results. Finally, we give some concluding remarks in Section 6.


## 2. Fractional BDFs

We will construct ak-step Continuous Fractional BDF (CFBDF) which will be used to obtain the discrete fractional BDF (FBDF). The CFBDF has the following general form:(4)
U
x
=
∑
j
=
0
k
-
1
γ
j
x
y
n
+
j
+
h
α
β
k
x
f
n
+
k
,where γ
j
(
x
) and β
k
(
x
) are continuous coefficients. We assume that y
n
+
j
=
U
(
x
n
+
j
h
) is the numerical approximation to the analytical solution y
(
x
n
+
j
) and f
n
+
j
=
D
α
U
(
x
n
+
j
h
) is an approximation to D
α
y
(
x
n
+
j
). The CFBDF is constructed from its equivalent form by requiring that the exact solution y
(
x
) is locally approximated by function (4) on the interval [
x
n
,
x
n
+
k
]. Next, we discuss the construction of the CFBDF in the following theorem.Theorem 1.
Let (4) satisfy the following conditions:(5)
U
x
n
+
j
=
y
n
+
j
,
j
=
0
1
k
-
1
,
D
α
U
x
n
+
k
=
f
n
+
k
. Then continuous representation (4) is equivalent to(6)
U
x
=
∑
j
=
0
k
det
⁡
W
j
det
⁡
W
P
j
x
, where we define matrix W as(7)
W
=
P
0
x
n
⋯
P
k
x
n
P
0
x
n
+
1
⋯
P
k
x
n
+
1
⋮
⋮
⋮
P
0
x
n
+
k
-
1
⋯
P
k
x
n
+
k
-
1
D
α
P
0
x
n
+
k
⋯
D
α
P
k
x
n
+
k
.
W
j is obtained by replacing the jth column of W by V where T denotes the transpose, P
j
(
x
)
=
x
j, j
=
0
(
1
)
k, are basis functions, and V is a vector given by (8)
V
=
y
n
,
y
n
+
1
,
y
n
+
2
,
…
,
y
n
+
k
-
1
,
f
n
+
k
T
.Proof.
We require that method (4) be defined by the assumed polynomial basis functions:(9)
γ
j
x
=
∑
i
=
0
k
γ
i
+
1
,
j
P
i
x
,
j
=
0
1
k
-
1
,
h
α
β
k
x
=
∑
i
=
0
k
h
α
β
i
+
1
,
k
P
i
x
, where γ
i
+
1
,
j and h
α
β
i
+
1
,
k are coefficients to be determined.
Substituting (9) into (4), we have (10)
U
x
=
∑
i
=
0
k
∑
j
=
0
k
-
1
γ
i
+
1
,
j
P
i
x
y
n
+
j
+
∑
i
=
0
k
h
α
β
i
+
1
,
k
P
i
x
f
n
+
k which may be written as (11)
U
x
=
∑
i
=
0
k
∑
j
=
0
k
-
1
γ
i
+
1
,
j
y
n
+
j
+
h
α
β
i
+
1
,
k
f
n
+
k
P
i
x and expressed as(12)
U
x
=
∑
i
=
0
k
l
i
P
i
x
, where (13)
l
i
=
∑
j
=
0
k
-
1
γ
i
+
1
,
j
y
n
+
j
+
h
α
β
i
+
1
,
k
f
n
+
k
. By imposing condition (5) on (12), we obtain a system of (
k
+
1
) equations, which can be expressed as W
=
L
V where L
=
l
0
,
l
1
,
…
,
l
k
T is a vector of (
k
+
1
) undetermined coefficients. Using Crammer’s rule, the elements of L can be obtained and are given by (14)
l
i
=
det
⁡
W
j
det
⁡
W
,
j
=
0
1
k
, where W
j is obtained by replacing the jth column of W by V. We rewrite (12) using the newly found elements of L as(15)
U
x
=
∑
j
=
0
k
det
⁡
W
j
det
⁡
W
P
j
x
.Remark 2.
Continuous scheme (4) which is equivalent to (6) is evaluated at x
n
+
k to obtain the k-step FBDF of the following form:(16)
y
n
+
k
-
∑
j
=
0
k
-
1
γ
k
j
y
n
+
j
=
h
α
β
k
k
f
n
+
k
.Also, we emphasize that continuous scheme (6) is used to obtain U
′
(
x
) and evaluated at x
n
+
i,  i
=
1
(
1
)
(
k
-
1
), to obtain(17)
h
α
f
n
+
1
-
∑
j
=
0
k
-
1
γ
1
j
y
n
+
j
=
h
α
β
1
k
f
n
+
k
h
α
f
n
+
2
-
∑
j
=
0
k
-
1
γ
2
j
y
n
+
j
=
h
α
β
2
k
f
n
+
k
⋮
h
α
f
n
+
k
-
1
-
∑
j
=
0
k
-
1
γ
k
-
1
j
y
n
+
j
=
h
α
β
k
-
1
k
f
n
+
k which, together with (16), forms the Block FBDF (BFBDF) which may be written in the following form:(18)
A
Y
¯
-
h
α
B
F
Y
¯
+
C
=
0
, where Y
¯
=
y
n
+
1
,
y
n
+
2
,
…
,
y
n
+
k
T, F
(
Y
¯
)
=
f
n
+
1
,
f
n
+
2
,
…
,
f
n
+
k
T, A and B are the coefficients of the formulas in (16) and (17), and C is a vector of initial conditions.Remark 3.
We note that it is possible to construct method (6) using other bases such exponential and trigonometric functions. However, (6) is constructed using polynomial basis functions, since methods produced via polynomial basis functions are easier to analyze. In fact, using other bases for the construction of (6) has the disadvantage of introducing additional parameters which makes the analysis of the methods produced more cumbersome. Moreover, the polynomial basis functions are appropriate for the construction of this class of methods since other bases can be written as polynomials in x via the Taylor series expansion.


## 3. Convergence of Methods

In this section, we will discuss the convergence of the methods in the following theorem.Theorem 4.
Letf
(
x
,
y
) be Lipschitz continuous with respect to y in a region D defined by a
≤
x
≤
b and -
∞
<
y
<
∞, where a and b are finite. Let (16) and (17) be constructed in such a way that β
i
k
=
O
(
n
α
-
1
); then (18) is said to be convergent if, for all initial value problems (1), we have that (19)
Y
¯
-
Y
≤
K
·
x
α
-
1
h
p
,
x
=
n
h
.for all x
∈
[
a
,
b
] and where constant K does not depend on h and p is the order of the method and h
>
0 is sufficiently small.Proof.
The exact form of the system formed by (16) and (17) is given by(20)
A
Y
-
h
α
B
F
Y
+
C
+
L
h
=
0
, where L
(
h
) is the truncation error vector of the formula in (18), Y
=
y
(
x
n
+
1
)
,
y
(
x
n
+
2
)
,
…
,
y
(
x
n
+
k
)
T, and F
(
Y
)
=
f
(
x
n
+
1
,
y
(
x
n
+
1
)
)
,
f
(
x
n
+
2
,
y
(
x
n
+
2
)
)
,
…
,
f
(
x
n
+
k
,
y
(
x
n
+
k
)
)
T.
The approximate form of the system is given by(21)
A
Y
¯
-
h
α
B
F
Y
¯
+
C
=
0
, where Y
¯ is the approximate solution of vector Y.
Subtracting (20) from (21), we obtain the following error system: (22)
A
E
=
h
α
B
F
Y
¯
-
F
Y
+
L
h
, where E
=
Y
¯
-
Y. Let K
1 be the Lipschitz constant of f; then (23)
A
E
≤
h
α
B
·
K
1
·
E
+
L
h
,
E
≤
A
-
h
α
B
·
K
1
-
1
L
h
. Let P
=
-
h
α
B
·
K
1, so that (24)
E
≤
A
+
P
-
1
L
h
=
O
n
h
α
-
1
·
K
h
p
.
(
(
A
+
P
)
-
1 exists since it is a monotone matrix (see [13])).Definition 5.
BFAMM (18) is said to be (a)
stable if and only if, for anyy
0
∈
R
m, there exists K
2
>
0 such that the solution of (1) satisfies y
(
x
n
)
≤
K
2 for n
≥
1,
(b)
asymptotically stable if and only if, for anyy
0
∈
R
m, the solution of (1) satisfies y
(
x
n
)
→
0 as n
→
∞.


## 4. Stability Properties of the Methods

To study the stability properties of BFBDF (18), we consider the following linear test problem:(25)
D
α
y
x
=
λ
y
x
,
λ
∈
C
,
0
<
α
<
1
,
y
x
0
=
y
0
, whose exact solution can be expressed in terms of the Mittag-Leffler function E
α
(
x
)
=
∑
k
=
0
∞
x
k
/
Γ
(
α
k
+
1
), as y
(
x
)
=
E
α
(
λ
(
x
-
x
0
)
α
)
y
0.We rewrite (18) in the following form:(26)
A
0
Y
γ
+
A
1
Y
γ
-
1
=
B
0
F
γ
, where Y
γ
=
y
n
+
1
,
y
n
+
2
,
…
,
y
n
+
k
T, Y
γ
-
1
=
y
n
-
k
-
1
,
…
,
y
n
-
1
,
y
n
T, and F
γ
=
f
n
+
1
,
f
n
+
2
,
…
,
f
n
+
k
T. Applying (26) to the linear test equation, we have the following linear recurrence relation:(27)
Y
γ
=
M
q
Y
γ
-
1
,
q
=
λ
h
α
, where(28)
M
q
=
A
0
-
q
B
0
-
1
A
1
, where M
(
q
) is the amplification matrix which determines the stability of the method.We are interested in investigating the values ofh
α
λ for which the numerical solution of (1) given by (18) asymptotically vanishes as the true solution. We give below the following definitions.Definition 6.
Stability domainΩ of BFBDF (18) is the set of all q
=
h
α
λ
∈
C such that linear recurrence (18) is asymptotically stable.Definition 7.
The stability domain of BFBDF (18) is given by the following set: (29)
Ω
=
q
:
ρ
q
<
1
,where ρ
(
q
) is the spectral radius of M
(
q
).Figures1, 2, 3, and 4 show the plot of the stability domain Ω of BFBDF (18) for some α for k
=
1 to 4
.Figure 1
The region of stability of the BFBDF fork
=
1 is to the left of the dividing line and is symmetric about the real axis.
(a)
α
=
0.25
(b)
α
=
0.5
(c)
α
=
0.75Figure 2
The region of stability of the BFBDF fork
=
2 is to the left of the dividing line and is symmetric about the real axis; the square and plus symbols represent the zeros and poles of the spectral radius of M
(
q
), respectively.
(a)
α
=
0.25
(b)
α
=
0.5
(c)
α
=
0.75Figure 3
The region of stability of the BFBDF fork
=
3 is to the left of the dividing line and is symmetric about the real axis; the square and plus symbols represent the zeros and poles of the spectral radius of M
(
q
), respectively.
(a)
α
=
0.25
(b)
α
=
0.5
(c)
α
=
0.75Figure 4
The region of stability of the BFBDF fork
=
4 is to the left of the dividing line and is symmetric about the real axis; the square and plus symbols represent the zeros and poles of spectral radius of M
(
q
), respectively.
(a)
α
=
0.25
(b)
α
=
0.5
(c)
α
=
0.75
(d)
α
=
0.80
(e)
α
=
0.85Remark 8.
The 1-step BFBDF isA-stable for all values of α
∈
(
0,1
) while the 2-step and 3-step BFBDF are A-stable for values of α
≥
0.75. Also, the 4-step BFBDF is A-stable for α
≥
0.80.


### 4.1. Implementation

Conventionally, FBDF (16) requires k initial conditions which are usually provided using a one-step method (like the Runge-Kutta method). However, in this paper, we construct additional methods from continuous scheme (15) which are implemented together with FBDF (16) to obtain approximate solutions to the exact solution of (1) without requiring starting values and predictors. For instance, if n
=
0 and k
=
2, then y
1
,
y
2
T are simultaneously obtained over the interval [
x
0
,
x
2
] as y
0 is known from the IVP. Similarly, if n
=
1 and k
=
3, then y
4
,
y
5
,
y
6
T are simultaneously obtained over the interval [
x
3
,
x
6
] and y
3 is known from the previous block and so on, until we reach the final subinterval [
x
N
-
k
,
x
N
].


## 4.1. Implementation

Conventionally, FBDF (16) requires k initial conditions which are usually provided using a one-step method (like the Runge-Kutta method). However, in this paper, we construct additional methods from continuous scheme (15) which are implemented together with FBDF (16) to obtain approximate solutions to the exact solution of (1) without requiring starting values and predictors. For instance, if n
=
0 and k
=
2, then y
1
,
y
2
T are simultaneously obtained over the interval [
x
0
,
x
2
] as y
0 is known from the IVP. Similarly, if n
=
1 and k
=
3, then y
4
,
y
5
,
y
6
T are simultaneously obtained over the interval [
x
3
,
x
6
] and y
3 is known from the previous block and so on, until we reach the final subinterval [
x
N
-
k
,
x
N
].


## 5. Numerical Examples

We validate our theoretical results from the previous sections by considering the following examples, which were solved using the BFBDF using a written code in Mathematica. The maximum errors are obtained for different step sizes in the interval of integration. We have solved two scalar examples, one system, one linear, and one nonlinear heat-type fractional differential equation.Example 1.
We consider the following problem given in [14]:(30)
D
α
y
x
=
-
y
+
x
2
-
x
+
2
x
2
-
α
Γ
3
-
α
+
x
1
-
α
Γ
2
-
α
,
0
<
α
<
1
,
0
≤
x
≤
1
,
y
0
=
0
,
Exact:
y
x
=
x
2
-
x
.Tables1, 2, and 3 show the numerical results for Example 1 using different values of α. It is evident from the table that the BFBDF performs favorably well with smaller number of steps.Table 1
Maximum errors using the BFBDFk
=
2 for Example 1.
N
α
=
0.25
α
=
0.50
α
=
0.75
10
76.661
e
-
16
8.882
e
-
16
3.1091
e
-
15
20
4.441
e
-
16
9.992
e
-
16
1.998
e
-
15
40
9.437
e
-
16
1.499
e
-
15
9.992
e
-
16
80
2.220
e
-
15
1.499
e
-
15
3.386
e
-
15Table 2
Maximum errors using the BFBDFk
=
3 for Example 1.
N
α
=
0.25
α
=
0.50
α
=
0.75
10
8.882
e
-
16
3.071
e
-
15
2.442
e
-
15
20
2.648
e
-
15
2.387
e
-
15
2.054
e
-
15
40
1.724
e
-
15
6.085
e
-
15
2.866
e
-
15
80
3.109
e
-
15
1.263
e
-
14
1.379
e
-
14Table 3
Maximum errors using the BFDBFk
=
4 for Example 1.
N
α
=
0.25
α
=
0.50
α
=
0.75
10
3.997
e
-
15
9.869
e
-
15
2.663
e
-
12
20
3.425
e
-
15
8.040
e
-
14
2.988
e
-
12
40
8.527
e
-
14
5.887
e
-
13
3.077
e
-
12
80
1.337
e
-
12
4.041
e
-
12
9.291
e
-
12Example 2.
We also consider the following FIVP with variable coefficients:(31)
D
α
y
x
-
x
y
=
Γ
2
α
+
1
Γ
α
+
1
x
α
-
x
-
x
2
α
+
1
,
0
<
α
<
1
,
0
≤
x
≤
1
,
y
0
=
1
,
Exact:
y
x
=
1
+
x
2
α
.This example was chosen to show the performance of the BFBDF on FIVP with variable coefficients. Tables4, 5, and 6 show the numerical results for Example 2 using different values of α. It is evident from the table that the BFBDF performs favorably well with smaller number of steps.Table 4
Maximum errors using the BFDBFk
=
2 for Example 2.
N
α
=
0.25
α
=
0.50
α
=
0.75
10
4.723
e
-
01
1.243
e
-
15
2.563
e
-
02
20
5.122
e
-
01
6.661
e
-
15
2.900
e
-
02
40
5.281
e
-
01
8.660
e
-
15
2.931
e
-
02
80
5.320
e
-
02
1.776
e
-
14
2.894
e
-
02Table 5
Maximum errors using the BFBDFk
=
3 for Example 2.
N
α
=
0.25
α
=
0.50
α
=
0.75
10
8.047
e
-
02
1.688
e
-
14
7.124
e
-
03
20
1.322
e
-
01
5.240
e
-
14
9.719
e
-
03
40
1.587
e
-
01
6.239
e
-
14
1.036
e
-
02
80
1.796
e
-
01
4.352
e
-
13
1.074
e
-
02Table 6
Maximum errors using the BFBDFk
=
4 for Example 2.
N
α
=
0.25
α
=
0.50
α
=
0.75
10
3.761
e
-
01
1.165
e
-
11
5.805
e
-
03
20
9.320
e
-
02
1.075
e
-
11
5.008
e
-
03
40
1.123
e
-
01
5.349
e
-
11
5.606
e
-
03
80
1.242
e
-
01
2.826
e
-
10
5.740
e
-
03Example 3.
We also consider the system of FIVP given in [15]:(32)
D
α
y
1
x
=
y
1
x
+
y
2
x
,
D
α
y
2
x
=
-
y
1
x
+
y
2
x
,
0
<
α
<
1
,
0
≤
x
≤
1
,
y
1
0
=
0
,
y
2
0
=
1
,
Exact:
y
1
x
=
e
x
sin
⁡
x
,
y
2
x
=
e
x
cos
⁡
x
for
α
=
1
.This example was chosen to show that the BFBDF performs well on a system as demonstrated forα
=
1. We note that for α
=
0.85 the solutions produced by the BFBDF are in agreement with the expected behavior of the exact solution. Details of the results are displayed in Table 7 and Figures 5 and 6.Table 7
Maximum errors using the BFBDFα
=
1 for Example 3.
N
k
=
2
k
=
3
k
=
4
5
6.647
e
-
02
6.033
e
-
03
3.991
e
-
03
10
5.276
e
-
03
7.498
e
-
04
1.382
e
-
04
20
9.839
-
04
9.430
e
-
05
2.814
e
-
06
40
2.073
e
-
04
1.174
e
-
05
1.727
e
-
07
80
4.717
e
-
05
1.446
e
-
06
6.489
e
-
06Figure 5
Graphical evidence for Example3, N
=
10.
(a)
Exact forα
=
1
(b)
BFBDF forα
=
1, k
=
2
(c)
BFBDF forα
=
1, k
=
3
(d)
BFBDF forα
=
1, k
=
4Figure 6
Graphical evidence for Example3, N
=
10.
(a)
BFBDF forα
=
0.85, k
=
2
(b)
BFBDF forα
=
0.85, k
=
3
(c)
BFBDF forα
=
0.85, k
=
4Example 4.
We also consider the following one-dimensional fractional heat-like problem given in [16]:(33)
∂
α
u
∂
t
α
=
x
2
2
∂
2
u
∂
x
2
,
0
<
α
<
1
,
0
≤
x
≤
1subject to the initial/boundary conditions u
(
x
,
0
)
=
x
2, u
(
0
,
t
)
=
0, u
(
1
,
t
)
=
e
t, and t
≥
0. The exact solution u
(
x
,
t
)
=
x
2
(
1
+
t
α
/
Γ
(
α
+
1
)
+
t
2
α
/
Γ
(
2
α
+
1
)
+
t
3
α
/
Γ
(
3
α
+
1
)
+
⋯
)
.In order to solve this PDE using the BFBDF, we carry out the semidiscretization of the spatial variablex using the second-order finite difference method to obtain the following first-order system in the second variable t:(34)
∂
α
u
m
∂
t
α
-
x
m
2
u
m
+
1
-
2
u
m
+
u
m
-
1
Δ
x
2
=
g
m
,
m
=
1
,
…
,
M
-
1
,
u
x
m
,
0
=
x
m
2
, where Δ
x
=
(
b
-
a
)
/
M, x
m
=
a
+
m
Δ
x, m
=
0,1
,
…
,
M, u
=
u
1
t
,
…
,
u
M
t
T, g
=
g
1
t
,
…
,
g
m
t
T, u
m
(
t
)
≈
u
(
x
m
,
t
), and g
m
(
t
)
≈
g
(
x
m
,
t
)
=
0, which can be written in the following form:(35)
u
α
=
f
t
,
u
,subject to the boundary conditions u
(
t
0
)
=
u
0, where f
(
t
,
u
)
=
A
u
+
g andA is M
×
M matrix arising from the semidiscretized system and g is a vector of constants.This example was chosen to show that the BFBDF performs well on a one-dimensional fractional heat-like problem with variable coefficients as demonstrated forα
=
1. We note that for α
=
0.75 the solutions produced by the BFBDF are in agreement with the expected behavior of the exact solution. Details of the results are displayed in Figures 7 and 8.Figure 7
Graphical evidence for Example4, N
=
10.
(a)
Error forα
=
1, k
=
2
(b)
Exact forα
=
1
(c)
BFBDF forα
=
1, k
=
2
(d)
Exact forα
=
0.75
(e)
BFBDF forα
=
0.75, k
=
2Figure 8
Graphical evidence for Example4, N
=
10.
(a)
Error forα
=
1, k
=
3
(b)
Exact forα
=
1
(c)
BFBDF forα
=
1, k
=
3
(d)
Exact forα
=
0.75
(e)
BFBDF forα
=
0.75
,
k
=
3Example 5.
We consider the following fractional Burger’s equation given in [17]:(36)
∂
α
u
∂
t
α
+
∂
u
2
∂
x
-
∂
2
u
∂
x
2
=
0
,
0
<
α
≤
1
,
0
<
x
<
1subject to the initial/boundary conditions u
(
x
,
0
)
=
x, u
(
0
,
t
)
=
0, u
(
1
,
t
)
=
1
/
1
+
t, and t
≥
0. The exact solution u
(
x
,
t
)
=
x
/
1
+
t,  α
=
1
.In order to solve this PDE using the BFBDF, we carry out the semidiscretization of the spatial variablex using the second-order finite difference method to obtain the following first-order system in the second variable t:(37)
∂
α
u
m
∂
t
α
-
x
m
2
u
m
+
1
-
2
u
m
+
u
m
-
1
Δ
x
2
=
g
m
,
m
=
1
,
…
,
M
-
1
,
u
x
m
,
0
=
x
m
2
, where Δ
x
=
(
b
-
a
)
/
M, x
m
=
a
+
m
Δ
x, m
=
0,1
,
…
,
M, u
=
[
u
1
(
t
)
,
…
,
u
M
(
t
)
]
T, g
=
[
g
1
(
t
)
,
…
,
g
m
(
t
)
]
T, u
m
(
t
)
≈
u
(
x
m
,
t
), and g
m
(
t
)
≈
g
(
x
m
,
t
)
=
0, which can be written in the following form:(38)
u
α
=
f
t
,
u
,subject to the boundary conditions u
(
t
0
)
=
u
0, where f
(
t
,
u
)
=
A
u
+
g and A is M
×
M matrix arising from the semidiscretized system and g is a vector of constants.This example was chosen to show that the BFBDF performs well on the Burger’s equation as demonstrated forα
=
1. We note that for α
=
0.75 the solutions produced by the BFBDF are in agreement with the expected behavior of the exact solution. Details of the results are displayed in Figures 9 and 10.Figure 9
Graphical evidence for Example5, N
=
10.
(a)
Error forα
=
1, k
=
2
(b)
Exact forα
=
1
(c)
BFBDF forα
=
1, k
=
2
(d)
BFBDF forα
=
0.75, k
=
1
(e)
BFBDF forα
=
0.75, k
=
2Figure 10
Graphical evidence for Example5, N
=
10.
(a)
Error forα
=
1, k
=
3
(b)
Exact forα
=
1
(c)
BFBDF forα
=
1, k
=
3
(d)
BFBDF forα
=
0.75, k
=
3
(e)
BFBDF forα
=
0.75
,
k
=
4


### 5.1. Block versus Predictor-Corrector Implementations

In order to demonstrate the superiority of the block implementation of the BFBDF over its predictor-corrector (PC) implementation, we have solved Example1 for k
=
1,2
,
3,4 and the results are displayed in Table 8. We constructed Fractional Explicit Adams Methods and used them as predictors for the FBDF. In Table 8, we observe that as the step length decreases the PC mode gives inaccurate approximations, while the BFBDF still retains its high accuracy. With these observations, we conclude that the BFBDF can be used as a general purpose method for the integration of fractional differential equations.Table 8
Comparison of exact errors for Example1 (α
=
0.75).
N
k
=
1
k
=
2
k
=
3
k
=
4
Block
PC
Block
PC
Block
PC
Block
PC
10
9.434
e
-
02
8.054
e
-
02
3.107
e
-
15
1.611
e
-
12
2.442
e
-
15
6.158
e
-
16
2.663
e
-
12
5.480
e
-
10
20
1.209
e
-
01
1.143
e
-
01
1.887
e
-
15
4.733
e
-
10
2.054
e
-
15
1.409
e
-
15
2.988
e
-
12
4.590
e
-
07
40
1.344
e
-
01
1.306
e
-
01
1.110
e
-
15
7.388
e
-
09
2.866
e
-
15
5.560
e
-
14
3.077
e
-
12
6.973
e
-
05
80
1.412
e
-
01
1.389
e
-
01
3.776
e
-
15
3.203
e
-
06
1.379
e
-
14
4.492
e
-
13
7.365
e
-
12
5.914
e
-
03


## 5.1. Block versus Predictor-Corrector Implementations

In order to demonstrate the superiority of the block implementation of the BFBDF over its predictor-corrector (PC) implementation, we have solved Example1 for k
=
1,2
,
3,4 and the results are displayed in Table 8. We constructed Fractional Explicit Adams Methods and used them as predictors for the FBDF. In Table 8, we observe that as the step length decreases the PC mode gives inaccurate approximations, while the BFBDF still retains its high accuracy. With these observations, we conclude that the BFBDF can be used as a general purpose method for the integration of fractional differential equations.Table 8
Comparison of exact errors for Example1 (α
=
0.75).
N
k
=
1
k
=
2
k
=
3
k
=
4
Block
PC
Block
PC
Block
PC
Block
PC
10
9.434
e
-
02
8.054
e
-
02
3.107
e
-
15
1.611
e
-
12
2.442
e
-
15
6.158
e
-
16
2.663
e
-
12
5.480
e
-
10
20
1.209
e
-
01
1.143
e
-
01
1.887
e
-
15
4.733
e
-
10
2.054
e
-
15
1.409
e
-
15
2.988
e
-
12
4.590
e
-
07
40
1.344
e
-
01
1.306
e
-
01
1.110
e
-
15
7.388
e
-
09
2.866
e
-
15
5.560
e
-
14
3.077
e
-
12
6.973
e
-
05
80
1.412
e
-
01
1.389
e
-
01
3.776
e
-
15
3.203
e
-
06
1.379
e
-
14
4.492
e
-
13
7.365
e
-
12
5.914
e
-
03


## 6. Conclusion

We have constructed a family of fractional linear multistep methods via the interpolation and collocation techniques. The methods developed are implemented as block methods for the numerical approximation of FIVPs. The stability properties of the methods are discussed and numerical examples are given to show that the methods are accurate and efficient, even when applied to large systems arising from the semidiscretization of one-dimensional fractional heat-like partial differential equations.



---

*Source: 650425-2015-09-29.xml*