# Solving Fractional Diffusion Equation via the Collocation Method Based on Fractional Legendre Functions

Authors: Muhammed Syam; Mohammed Al-Refai
Journal: Journal of Computational Methods in Physics
 (2014)

---

## Abstract

A formulation of the fractional Legendre functions is constructed to solve the generalized time-fractional diffusion equation. The fractional derivative is described in the Caputo sense. The method is based on the collection Legendre and path following methods. Analysis for the presented method is given and numerical results are presented.


---

## Body

## 1. Introduction

We consider the generalized time-fractional diffusion equation of the form(1)
D
t
α
u
(
x
,
t
)
=
a
(
x
,
t
)
D
x
2
u
(
x
,
t
)
+
f
(
x
,
t
)
,
m
m
m
m
m
l
m
m
m
x
∈
(
-
1,1
)
,
t
∈
(
0
,
T
)
,
with initial and boundary conditions
(2)
u
(
-
1
,
t
)
=
h
1
(
t
)
,
u
(
1
,
t
)
=
h
2
(
t
)
,
u
(
x
,
0
)
=
g
(
x
)
,
where a, f
∈
C
1
(
[
-
1,1
]
×
[
0
,
T
]
), g
∈
C
[
-
1,1
], h
1, h
2
∈
C
[
0
,
T
], T
>
0, and 0
<
α
≤
1. For α
=
1, the fractional diffusion equation is reduced to a conventional diffusion-reaction equation which is well studied, so we focus on 0
<
α
<
1. Some existence and uniqueness results of Problem (1)-(2) were established in [1]. In recent years, great interests were devoted to the analytical and numerical treatments of fractional differential equations (FDEs). Usually, FDEs appear as generalizations to existing models with integer derivative and they also present new models for some physical problems [2, 3]. In general, FDEs do not possess exact solutions in closed forms, and, therefore, numerical methods such as the variational iteration (VIM) [4, 5], the homotopy analysis method (HAM) [6, 7], and the Adomian decomposition method (ADM) [8, 9] have been implemented for several types of FDEs. Also, the maximum principle and the method of lower and upper solutions have been extended to deal with FDEs and obtain analytical and numerical results [10, 11]. The Tau method, the pseudospectral method, and the wavelet method based on the Legendre polynomials have been implemented for several types of FDEs [12–14]. Kazem et al. [12] have constructed the Legendre functions of fractional order and discussed some of their properties. The resulting Legendre function operational and product matrices, together with the Tau method, have been implemented to solve linear and nonlinear fractional differential equations. The effectiveness of the approach has been examined through several examples. In [13], a fractional diffusion equation is considered, where the fractional derivative of order 1
<
α
≤
2 refers to the spatial variable x. The Legendre pseudospectral method is implemented to solve the problem, where the solution is expanded with regular Legendre polynomials. As a result, a system of linear equation has been obtained and integrated using the finite difference method. However, in solving fractional differential equations of order α using series expansions, it is common and more efficient to expand the solution with fractional functions of the form ∑
i
=
0
N
‍
c
i
x
i
α. Rawashdeh [14] has implemented the Legendre wavelets method for integrodifferential equations with fractional order.The Legendre collocation method has been implemented for wide classes of differential equations and the effectiveness of the method is illustrated [15]. In the recent work, we intend to apply the collocation method based on the shifted fractional Legendre functions to integrate the Problem (1)-(2). To the best of our knowledge, the method has not been developed to integrate fractional diffusion equations of the form (1)-(2). We organize this paper as follows. In Section 2, we present basic definitions and results of fractional derivative. In Section 3, we present the numerical technique for solving Problem (1)-(2). In Section 4, we present some numerical results to illustrate the efficiency of the presented method. Finally we conclude with some comments in Section 5.


## 2. Preliminaries

In this section, we present the definition and some preliminary results of the Caputo fractional derivative, as well as the definition of the fractional-order Legendre functions and their properties.Definition 1.
A real functionf
(
t
), t
>
0, is said to be in the space C
μ, μ
∈
R, if there exists a real number p
>
μ, such that f
(
t
)
=
t
p
f
1
(
t
), where f
1
(
t
)
∈
C
[
0
,
∞
), and it is said to be in the space C
μ
m if f
(
m
)
∈
C
μ, m
∈
N.Definition 2.
The left Riemann-Liouville fractional integral of orderδ
≥
0, of a function f
∈
C
μ, μ
≥
-
1, is defined by
(3)
I
δ
f
(
t
)
=
{
1
Γ
(
δ
)
∫
0
t
‍
(
t
-
s
)
δ
-
1
f
(
s
)
d
s
,
δ
>
0
,
f
(
t
)
,
δ
=
0
.Definition 3.
Forδ
>
0, m
-
1
<
δ
≤
m, m
∈
N, t
>
0, and f
∈
C
-
1
m, the left Caputo fractional derivative is defined by
(4)
D
δ
f
(
t
)
=
{
1
Γ
(
m
-
δ
)
∫
0
t
‍
(
t
-
s
)
m
-
1
-
δ
f
(
m
)
(
s
)
d
s
,
m
-
1
<
δ
<
m
,
f
(
m
)
(
t
)
,
δ
=
m
,
where Γ is the well-known Gamma function.The Caputo derivative defined in (4) is related to the Riemann-Liouville fractional integral, I
δ, of order δ
∈
R
+, by D
δ
f
(
t
)
=
I
m
-
δ
f
(
m
)
(
t
). The Caputo fractional derivative satisfies the following properties for f
∈
C
μ, μ
≥
-
1, and α
≥
0 (see [16]):(1)
D
α
I
α
f
(
t
)
=
f
(
t
),
(2)
I
α
D
α
f
(
t
)
=
f
(
t
)
-
∑
k
=
0
n
-
1
‍
f
(
k
)
(
0
)
(
t
k
/
k
!
),
(3)
D
α
c
=
0, where c is constant,
(4)
D
α
t
γ
=
{

0
,
γ
<
α
,
γ
∈
{
0,1
,
2
,
…
}
(
Γ
(
γ
+
1
)
/
Γ
(
γ
-
α
+
1
)
)
t
γ
-
α
,
otherwise

},
(5)
D
α
(
∑
k
=
0
m
‍
c
i
f
i
(
t
)
)
=
∑
k
=
0
m
‍
c
i
D
α
f
i
(
t
), where c
1
,
c
2
,
…
,
c
m are constants.The basic concept of this paper is the Legendre polynomials. For this reason, we study some of their properties.Definition 4.
The Legendre polynomials{
L
k
(
x
)
:
k
=
0,1
,
2
,
…
} are the eigenfunctions of the Sturm-Liouville problem:
(5)
(
(
1
-
x
2
)
L
k
′
(
x
)
)
′
+
k
(
k
+
1
)
L
k
(
x
)
=
0
,
x
∈
[
-
1,1
]
.Among the properties of the Legendre polynomials we list the following [19]:(1)
∫
-
1
1
‍
L
i
(
x
)
L
j
(
x
)
d
x
=
(
2
/
(
2
i
+
1
)
)
δ
i
j, where δ
i
j
=
{

0
,
i
≠
j
1
,
i
=
j

},
(2)
L
i
+
1
(
x
)
=
(
(
2
i
+
1
)
/
(
i
+
1
)
)
x
L
i
(
x
)
-
(
i
/
(
i
+
1
)
)
L
i
-
1
(
x
), for i
⩾
1,
(3)
L
i
(
±
1
)
=
(
±
1
)
i.In order to use these polynomials on the interval[
0,1
], we define the shifted Legendre polynomials by S
i
(
x
)
=
L
i
(
2
x
-
1
). Using the change of variable z
=
2
x
-
1, S
i
(
x
) has the following properties:(1)
∫
0
1
‍
S
i
(
x
)
S
j
(
x
)
d
x
=
(
1
/
(
2
i
+
1
)
)
δ
i
j,
(2)
S
i
+
1
(
x
)
=
(
(
2
i
+
1
)
/
(
i
+
1
)
)
(
2
x
-
1
)
S
i
(
x
)
-
(
i
/
(
i
+
1
)
)
S
i
-
1
(
x
), for i
⩾
1,
(3)
S
i
(
0
)
=
(
-
1
)
i and S
i
(
1
)
=
1.The analytic closed form of the shifted Legendre polynomials of degreei is given by
(6)
S
i
(
x
)
=
∑
k
=
0
i
‍
(
-
1
)
i
+
k
(
i
+
k
)
!
(
i
-
k
)
!
x
k
(
k
!
)
2
.One of the common and efficient methods for solving fractional differential equations of orderα
>
0 is using series expansion of the form ∑
k
=
0
i
‍
c
k
t
α
k. For this reason, we define the fractional-order Legendre function by F
i
α
(
t
)
=
S
i
(
t
α
). Using the properties of the shifted Legendre polynomials, it is easy to verify that [20](1)
(
(
t
-
t
1
+
α
)
F
i
α
′
(
t
)
)
′
+
i
(
i
+
1
)
α
2
t
α
-
1
F
i
α
(
t
)
=
0, t
∈
(
0,1
),
(2)
F
i
α
(
t
)
=
(
(
2
i
+
1
)
/
(
i
+
1
)
)
(
2
t
α
-
1
)
F
i
α
(
t
)
-
(
i
/
(
i
+
1
)
)
F
i
-
1
α
(
t
), for i
⩾
1,
(3)
F
0
α
(
t
)
=
1 and F
1
α
(
t
)
=
2
t
α
-
1,
(4)
F
i
α
(
0
)
=
(
-
1
)
i and F
i
α
(
1
)
=
1.In addition,{
F
i
α
(
t
)
:
i
=
0,1
,
2
,
…
} are orthogonal functions with respect to the weight function w
(
t
)
=
t
α
-
1 on (
0,1
) with
(7)
∫
0
1
‍
F
i
α
(
t
)
F
j
α
(
t
)
w
(
t
)
d
t
=
1
(
2
i
+
1
)
α
δ
i
j
.The closed form ofF
i
α
(
t
) is given by
(8)
F
i
α
(
t
)
=
∑
k
=
0
i
‍
(
-
1
)
i
+
k
(
i
+
k
)
!
(
i
-
k
)
!
t
α
k
(
k
!
)
2
.Using properties (4) and (5) of the Caputo fractional derivative, we have(9)
D
α
F
i
α
(
t
)
=
∑
k
=
1
i
‍
(
-
1
)
i
+
k
(
i
+
k
)
!
(
i
-
k
)
!
(
k
!
)
2
Γ
(
k
α
+
1
)
Γ
(
(
k
-
1
)
α
+
1
)
t
(
k
-
1
)
α
.The following result is important, since it facilitates applying the collection method.Theorem 5.
Letu
∈
C
1
[
0,1
] and let u
′′
(
t
) be a piecewise continuous function on [
0,1
]. Then,(i)
u
(
t
) can be represented by infinite series expansion as u
(
t
)
=
∑
k
=
0
∞
‍
u
k
F
k
α
(
t
), where
(10)
u
k
=
(
2
i
+
1
)
α
∫
0
1
‍
u
(
t
)
F
k
α
(
t
)
w
(
t
)
d
t
.
(ii)
∑
k
=
0
∞
‍
u
k
(
α
)
F
k
α
(
t
) converges uniformly on [
0,1
] to D
α
u
(
t
), where u
k
(
α
)
=
∑
j
=
k
+
1
∞
‍
a
j
k
u
j and a
j
k
=
(
2
k
+
1
)
α
∫
0
1
‍
D
α
F
j
α
(
t
)
F
k
α
(
t
)
w
(
t
)
d
t, for k
=
0,1
,
2
,
… and j
=
k
+
1
,
k
+
2
,
… .Proof.
(i) Sinceu
∈
C
1
[
-
1,1
] and u
′′
(
x
) is a piecewise continuous function on [
-
1,1
],  ∑
k
=
0
∞
v
k
L
k
(
x
) converges uniformly to u
(
x
) on [
-
1,1
], where {
v
k
} can be computed by the orthogonality relation of the Legendre polynomials; see [11]. Since r
:
[
0,1
]
→
[
-
1,1
], defined by r
(
t
)
=
2
t
α
-
1, is a bijection continuous function, the infinite series ∑
k
=
0
∞
u
k
F
k
α
(
t
) converges uniformly to u
(
t
) on [
0,1
], where the value of u
k follows from the orthogonality relation of {
F
i
α
(
t
)
:
i
=
0,1
,
2
,
…
} with respect to the weight function w
(
t
)
=
t
α
-
1 on [
0,1
], which completes the proof.
(ii) LetP
n
(
t
)
=
∑
k
=
0
n
u
k
F
k
α
(
t
) for n
=
0,1
,
2
,
…. From Part (i), P
n
(
t
) converges uniformly to u
(
t
) on [
0,1
]. Since u
∈
C
1
[
0,1
] and u
′′
(
t
) is a piecewise continuous function on [
0,1
],
(11)
d
d
t
(
Lim
n
→
∞
P
n
(
t
)
)
=
Lim
n
→
∞
(
d
d
t
P
n
(
t
)
)
,
and (
d
/
d
t
)
P
n
(
t
) converges uniformly to (
d
/
d
t
)
u
(
t
) on [
0,1
]. Thus, ∫
0
x
(
P
n
′
(
t
)
/
(
x
-
t
)
α
)
d
t converges uniformly to ∫
0
x
(
u
′
(
t
)
/
(
x
-
t
)
α
)
d
t on [
0,1
] which gives the result of the second part. The value of a
j
k follows from the orthogonality relation of {
F
i
α
(
t
)
:
i
=
0,1
,
2
,
…
} with respect to the weight function w
(
t
)
=
t
α
-
1 on [
0,1
].


## 3. Collocation Method

In this section, we use the fractional-order Legendre collocation method to discretize Problem (1)-(2). For simplicity, we assume that a
(
x
,
t
)
=
a
(
x
) and T
=
1. If T
≠
1, we use the change of variable τ
=
t
/
T to make the t-domain (0,1). Expand the solution u
(
x
,
t
) in terms of the fractional-order Legendre function as follows:
(12)
U
N
(
x
,
t
)
=
∑
k
=
0
N
+
2
‍
u
k
(
t
)
L
k
(
x
)
.
Thus,
(13)
D
x
2
U
N
(
x
,
t
)
=
∑
k
=
0
N
‍
u
k
(
2
)
(
t
)
L
k
(
x
)
,
where u
k
(
2
)
=
(
k
+
(
1
/
2
)
)
∑
j
=
k
+
2
,
j
+
k
even
[
j
(
j
+
1
)
-
k
(
k
+
1
)
]
u
j for k
=
0,1
,
…
,
N. Therefore, for U
N
(
x
,
t
), the residual is given by
(14)
R
(
U
N
(
x
,
t
)
)
≈
D
t
α
U
N
(
x
,
t
)
-
a
(
x
,
t
)
D
x
2
U
N
(
x
,
t
)
-
f
(
x
,
t
)
.Orthogonalize the residual with respect to the Dirac delta function as follows:(15)
〈
R
(
U
N
(
x
,
t
)
)
,
δ
(
x
-
x
j
)
〉
=
∫
0
1
R
(
U
N
(
x
,
t
)
)
δ
(
x
-
x
j
)
d
x
=
0
,
for
j
=
0
:
N
,
where x
j are the collocation points. We choose the collocation points to be the roots of L
N
+
2
′. Therefore, (15) leads to the elementwise equation:
(16)
D
t
α
U
N
(
x
j
,
t
)
-
a
(
x
j
,
t
)
D
x
′′
U
N
(
x
j
,
t
)
+
f
(
x
j
,
t
)
=
0
,
h
h
l
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
for
j
=
0
:
N
,
or
(17)
∑
k
=
0
N
+
2
‍
[
D
t
α
u
k
(
t
)
]
L
k
(
x
j
)
-
a
(
x
j
,
t
)
∑
k
=
0
N
‍
u
k
(
2
)
(
t
)
L
k
(
x
j
)
-
f
(
x
j
,
t
)
=
0
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
for
j
=
0
:
N
.
Let
(18)
U
(
t
)
=
[
u
0
(
t
)
u
1
(
t
)
⋮
u
N
+
2
(
t
)
]
,
U
(
2
)
(
t
)
=
[
u
0
(
2
)
(
t
)
u
1
(
2
)
(
t
)
⋮
u
N
(
2
)
(
t
)
]
,
F
(
t
)
=
[
f
(
x
0
,
t
)
f
(
x
1
,
t
)
⋮
f
(
x
N
,
t
)
]
.Thus, we can rewrite (17) in the matrix form as
(19)
A
1
D
t
α
U
(
t
)
-
A
A
2
U
(
2
)
(
t
)
-
F
(
t
)
=
0
,
where
(20)
A
1
=
[
L
0
(
x
0
)
L
1
(
x
0
)
⋯
L
N
+
2
(
x
0
)
L
0
(
x
1
)
L
1
(
x
1
)
⋯
L
N
+
2
(
x
1
)
⋮
⋮
⋱
⋮
L
0
(
x
N
)
L
1
(
x
N
)
⋯
L
N
+
2
(
x
N
)
]
,
A
=
[
a
(
x
0
)
0
⋯
0
0
a
(
x
1
)
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
a
(
x
N
)
]
,
A
2
=
[
L
0
(
x
0
)
L
1
(
x
0
)
⋯
L
N
(
x
0
)
L
0
(
x
1
)
L
1
(
x
1
)
⋯
L
N
(
x
1
)
⋮
⋮
⋱
⋮
L
0
(
x
N
)
L
1
(
x
N
)
⋯
L
N
(
x
N
)
]
.
Since u
k
(
2
)
=
(
k
+
(
1
/
2
)
)
∑
j
=
k
+
2
,
j
+
k
even
[
j
(
j
+
1
)
-
k
(
k
+
1
)
]
u
j for k
=
0,1
,
…
,
N, it is easy to see that U
(
2
)
(
t
)
=
A
3
U
(
t
), where A
3 is (
N
+
1
)
×
(
N
+
3
) matrix. Therefore, System (19) becomes
(21)
A
1
D
t
α
U
(
t
)
-
A
A
2
A
3
U
(
t
)
-
F
(
t
)
=
0
.Now, we study the boundary conditions on the variablex. From (12), one can see that
(22)
U
N
(
±
1
,
t
)
=
∑
k
=
0
N
+
2
‍
u
k
(
t
)
L
k
(
±
1
)
=
∑
k
=
0
N
+
1
‍
(
±
1
)
k
u
k
(
t
)
,
which implies that
(23)
A
4
U
(
t
)
=
[
h
1
(
t
)
h
2
(
t
)
]
,
where A
4
=
[
-
1
1
-
1
1
⋯
(
-
1
)
N
+
1
1
1
1
1
⋯
1
]
2
×
(
N
+
2
). From Systems (21) and (23), we obtain the following fractional system:
(24)
B
1
D
t
α
U
(
t
)
-
B
2
U
(
t
)
=
R
(
t
)
,
where
(25)
B
1
=
[
A
1
0
2
×
(
N
+
2
)
]
,
B
2
=
[
A
A
2
A
3
A
4
]
,
R
(
t
)
=
[
F
(
t
)
h
1
(
t
)
h
2
(
t
)
]
.From (2) and (12), we see that
(26)
g
(
x
)
=
U
N
(
x
,
0
)
=
∑
k
=
0
N
+
2
‍
u
k
(
0
)
L
k
(
x
)
.Using the orthogonality property of the Legendre polynomials, we get(27)
u
k
(
0
)
=
∫
-
1
1
g
(
x
)
L
k
(
x
)
d
x
∫
-
1
1
L
k
2
(
x
)
d
x
=
2
k
+
1
2
∫
-
1
1
g
(
x
)
L
k
(
x
)
d
x
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
l
h
h
h
h
h
h
h
l
for
k
=
0
:
N
+
2
.Therefore, our initial fractional value problem becomes(28)
B
1
D
t
α
U
(
t
)
-
B
2
U
(
t
)
=
R
(
t
)
,
(29)
U
(
0
)
=
B
,
where
(30)
B
=
[
1
2
∫
-
1
1
g
(
x
)
L
0
(
x
)
d
x
3
2
∫
-
1
1
g
(
x
)
L
1
(
x
)
d
x
⋮
2
N
+
5
2
∫
-
1
1
g
(
x
)
L
N
+
2
(
x
)
d
x
]
.To solve the System (28)-(29), we use the fractional-order Legendre collocation method. Approximate the solution U
(
t
) in terms of the fractional-order Legendre function as follows:
(31)
V
M
(
t
)
=
∑
k
=
0
M
+
1
‍
U
k
F
k
α
(
t
)
.
Thus,
(32)
D
t
α
V
M
(
x
,
t
)
=
∑
k
=
0
M
‍
U
k
(
α
)
(
t
)
F
k
α
(
t
)
,
where U
k
(
α
)
=
∑
j
=
k
+
1
M
+
1
a
j
k
U
j, and a
j
k
=
(
2
k
+
1
)
α
∫
0
1
D
α
F
j
α
(
t
)
F
k
α
(
t
)
w
(
t
)
d
t, for k
=
0,1
,
2
,
…
,
M and j
=
k
+
1
,
k
+
2
,
…
,
M
+
1. Therefore, for V
M
(
t
), the residual is given by
(33)
R
(
V
M
(
t
)
)
≈
B
1
D
t
α
V
M
(
t
)
-
B
2
V
M
(
t
)
-
R
(
t
)
.Orthogonalize the residual with respect to the Dirac delta function as follows:(34)
〈
R
(
V
M
(
t
)
)
,
δ
(
t
-
t
j
)
〉
=
∫
0
1
R
(
V
M
(
t
)
)
δ
(
t
-
t
j
)
d
t
=
0
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
for
j
=
0
:
M
,
where t
j are the collocation points. We choose the collocation points to be
(35)
{
(
t
0
+
1
2
)
-
α
,
(
t
1
+
1
2
)
-
α
,
…
,
(
t
M
+
1
2
)
-
α
}
,
where t
0
,
t
1
,
…
,
t
M are the roots of L
M
+
1
′
(
t
). Therefore, (34) leads to the elementwise equation:
(36)
B
1
D
t
α
V
M
(
t
j
)
-
B
2
V
M
(
t
j
)
-
R
(
t
j
)
=
0
,
for
j
=
0
:
M
or
(37)
B
1
∑
k
=
0
M
‍
U
k
(
α
)
(
t
)
F
k
α
(
t
)
-
B
2
∑
k
=
0
M
+
1
‍
U
k
F
k
α
(
t
)
-
R
(
t
j
)
=
0
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
for
j
=
0
:
M
.
Let
(38)
V
=
[
U
0
U
1
⋮
U
M
~
+
1
]
,
V
(
α
)
=
[
U
0
(
α
)
U
1
(
α
)
⋮
U
M
(
α
)
]
,
Λ
=
[
R
(
t
0
)
R
(
t
1
)
⋮
R
(
t
M
)
]
.Thus, we can rewrite (37) in the matrix form as
(39)
C
3
C
1
V
(
α
)
-
C
4
C
2
V
-
Λ
=
0
,
where(40)
C
1
=
[
B
1
0
⋯
0
0
B
1
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
B
1
]
(
N
+
3
)
(
M
+
1
)
×
(
N
+
3
)
(
M
+
1
)
,
C
2
=
[
B
2
0
⋯
0
0
B
2
⋯
0
⋮
⋮
⋱
⋮
0
0
⋯
B
2
]
(
N
+
3
)
(
M
+
2
)
×
(
N
+
3
)
(
M
+
2
)
,
C
3
=
[
F
0
α
(
t
0
)
I
N
+
3
F
1
α
(
t
0
)
I
N
+
3
⋯
F
M
α
(
t
0
)
I
N
+
3
F
0
α
(
t
1
)
I
N
+
3
F
1
α
(
t
1
)
I
N
+
3
⋯
F
M
α
(
t
1
)
I
N
+
3
⋮
⋮
⋱
⋮
F
0
α
(
t
M
)
I
N
+
3
F
1
α
(
t
M
)
I
N
+
3
⋯
F
M
α
(
t
M
)
I
N
+
3
]
(
N
+
3
)
(
M
+
1
)
×
(
N
+
3
)
(
M
+
1
)
,
C
4
=
[
F
0
α
(
t
0
)
I
N
+
3
F
1
α
(
t
0
)
I
N
+
3
⋯
F
M
+
1
α
(
t
0
)
I
N
+
3
F
0
α
(
t
1
)
I
N
+
3
F
1
α
(
t
1
)
I
N
+
3
⋯
F
M
+
1
α
(
t
1
)
I
N
+
3
⋮
⋮
⋱
⋮
F
0
α
(
t
M
)
I
N
+
3
F
1
α
(
t
M
)
I
N
+
3
⋯
F
M
+
1
α
(
t
M
)
I
N
+
3
]
(
N
+
3
)
(
M
+
1
)
×
(
N
+
3
)
(
M
+
2
)
.SinceU
k
(
α
)
=
∑
j
=
k
+
1
M
+
1
a
j
k
U
j, it is easy to see that V
(
α
)
=
C
5
V, where(41)
A
3
=
[
0
I
N
+
3
a
10
I
N
+
3
a
20
I
N
+
3
a
30
I
N
+
3
⋯
a
M
0
I
N
+
3
a
(
M
+
1
)
0
I
N
+
3
0
I
N
+
3
0
I
N
+
3
a
21
I
N
+
3
a
31
I
N
+
3
⋯
a
M
1
I
N
+
3
a
(
M
+
1
)
1
I
N
+
3
0
I
N
+
3
0
I
N
+
3
0
I
N
+
3
a
32
I
N
+
3
⋯
a
M
2
I
N
+
3
a
(
M
+
1
)
2
I
N
+
3
⋮
⋮
⋮
⋮
⋱
⋮
⋮
0
I
N
+
3
0
I
N
+
3
0
I
N
+
3
0
I
N
+
3
⋯
a
M
M
I
N
+
3
a
(
M
+
1
)
M
I
N
+
3
]
,and a
j
k
=
(
2
k
+
1
)
α
∫
0
1
D
α
F
j
α
(
t
)
F
k
α
(
t
)
w
(
t
)
d
t for k
=
0,1
,
2
,
…
,
M and j
=
k
+
1
:
M
+
1. Therefore, System (19) becomes
(42)
C
3
C
1
C
5
V
-
C
4
C
2
V
-
Λ
=
0
.Now, we study the initial condition on the variablex. From (29), one can see that
(43)
B
=
U
(
0
)
=
∑
k
=
0
M
+
1
‍
U
k
F
k
β
(
0
)
=
∑
k
=
0
M
+
1
‍
(
-
1
)
k
U
k
,
which implies that
(44)
C
6
V
=
B
,
where C
6
=
[
-
I
N
+
3
I
N
+
3
-
I
N
+
3
I
N
+
3
⋯
(
-
1
)
N
+
1
I
N
+
3
]
(
N
+
3
)
×
(
N
+
3
)
(
M
+
2
). From Systems (21) and (23), we obtain the following linear system:
(45)
Ω
V
=
Ψ
,
where
(46)
Ω
=
[
C
3
C
1
C
5
-
C
4
C
2
C
6
]
,
Ψ
=
[
Λ
B
]
.
Finally, we use Mathematica to solve the above linear system.


## 4. Numerical Results

In this section, we implement the proposed numerical technique for four examples.Example 1.
Consider the fractional diffusion equation:(47)
D
t
α
u
(
x
,
t
)
=
D
x
2
u
(
x
,
t
)
,
x
∈
(
-
1,1
)
,
t
∈
(
0,1
)
,
0
<
α
<
1
,
u
(
±
1
,
t
)
=
±
E
α
(
-
t
α
)
sin
⁡
(
1
)
,
0
<
t
<
1
,
u
(
x
,
0
)
=
sin
⁡
x
,
-
1
<
x
<
1
,
where E
α
(
t
) is the Mittag-Leffler function. Since D
t
α
E
α
(
-
t
α
)
=
-
E
α
(
-
t
α
) for 0
<
α
<
1, it is easy to see that the exact solution is
(48)
u
(
x
,
t
)
=
E
α
(
-
t
α
)
sin
⁡
x
.
The exact solution for α
=
1 is
(49)
u
(
x
,
t
)
=
e
-
t
sin
⁡
x
.
The approximate solutions generated by the proposed method are presented in Figure 1, for different values of α and N
=
M
=
10. Figure 2 depicts the exact solution (red) and the approximate solution (green) for α
=
1, and N
=
M
=
10. Define the error
(50)
Er
=
max
⁡
⁡
{
|
u
(
x
i
,
t
j
)
-
u
app
(
x
i
,
t
j
)
|
h
h
h
l
h
h
h
h
h
:
x
i
=
-
1
+
0.01
(
i
-
1
)
,
t
j
=
0.01
(
j
-
1
)
,
h
h
h
h
h
l
h
h
h
i
:
1
:
201
,
j
=
1
:
101
|
u
(
x
i
,
t
j
)
-
u
app
(
x
i
,
t
j
)
|
}
,
where u
app
(
x
,
t
) is the approximate solution generated by the proposed method for N
=
M
=
10. Table 1 presents the error for different values of α.Table 1
Error for Example1.
α
Er
0.5
4.4 ∗ 10−9
0.9
8.9 ∗ 10−10
0.99
6.1 ∗ 10−10
0.9999
5.2 ∗ 10−12The approximate solutions of Example1.
(a)
(b)
(c)
(d)Figure 2
Exact solution (red) and approximate solution (green) of Example1.Example 2.
Consider the fractional diffusion equation:(51)
D
t
α
u
(
x
,
t
)
=
D
x
2
u
(
x
,
t
)
-
2
t
3
+
6
Γ
(
4
-
α
)
t
3
-
α
x
2
,
h
h
h
h
h
h
h
h
h
h
l
x
∈
(
-
1,1
)
,
t
∈
(
0,1
)
,
0
<
α
<
1
,
u
(
±
1
,
t
)
=
t
3
,
0
<
t
<
1
,
u
(
x
,
0
)
=
0
,
-
1
<
x
<
1
,
where u
(
x
,
t
)
=
t
3
x
2 is the exact solution. The approximate solutions generated by the proposed method are presented in Figure 3 for different values of α and N
=
M
=
6. Figure 4 depicts the exact solution (red) and the approximate solution (green) for α
=
1 and N
=
M
=
6.
Table2 presents the error for different values of α.Table 2
Error for Example2.
α
Er
0.5
7.2 ∗ 10−10
0.9
3.3 ∗ 10−11
0.99
4.1 ∗ 10−12
0.9999
7.2 ∗ 10−13The approximate solutions of Example2.
(a)
(b)
(c)
(d)Figure 4
Exact solution (red) and approximate solution (green) for Example2.Example 3.
Consider the fractional diffusion equation presented in [17]:
(52)
D
t
α
u
(
x
,
t
)
=
D
x
2
u
(
x
,
t
)
+
2
x
(
1
-
x
)
t
2
-
α
Γ
(
3
-
α
)
+
2
(
t
2
+
1
)
,
h
l
h
h
h
h
h
h
h
h
h
h
h
h
h
h
x
∈
(
0,1
)
,
t
∈
(
0,1
)
,
0
<
α
<
1
,
u
(
0
,
t
)
=
u
(
1
,
t
)
=
0
,
0
<
t
<
1
,
u
(
x
,
0
)
=
x
(
x
-
1
)
,
-
1
<
x
<
1
.
The exact solution is
(53)
u
(
x
,
t
)
=
x
(
x
-
1
)
(
t
2
+
1
)
.
To apply the proposed method, we shall do the following change of variable X
=
2
x
-
1. In this case, the x-domain becomes [
-
1,1
]. The approximate solutions generated by the proposed method and the exact solution are presented in Figure 5 for α
=
0.92 and α
=
0.98 at t
=
0.01 and N
=
M
=
10. Table 3 presents a comparison between the error in our results and the ones obtained by the finite difference method (FDM) [17] for α
=
0.92,0.98 and t
=
0.01.Table 3
Errors in FDM [17] and the proposed method at t
=
0.01.
α
=
0
.
92
α
=
0
.
98
x
Error in FDM
Er
Error in FDM
Er
0.1
1.7 ∗ 10−6
3.1 ∗ 10−11
1.4 ∗ 10−6
2.9 ∗ 10−11
0.2
3.1 ∗ 10−6
4.2 ∗ 10−11
2.5 ∗ 10−6
2.9 ∗ 10−11
0.3
4.0 ∗ 10−6
4.7 ∗ 10−11
3.3 ∗ 10−6
3.1 ∗ 10−11
0.4
4.6 ∗ 10−6
4.9 ∗ 10−11
3.8 ∗ 10−6
3.3 ∗ 10−11
0.5
4.8 ∗ 10−6
5.2 ∗ 10−11
4.0 ∗ 10−6
3.6 ∗ 10−11
0.6
4.6 ∗ 10−6
3.9 ∗ 10−11
3.8 ∗ 10−6
3.3 ∗ 10−11
0.7
4.0 ∗ 10−6
3.2 ∗ 10−11
3.3 ∗ 10−6
3.1 ∗ 10−11
0.8
3.1 ∗ 10−6
2.8 ∗ 10−11
2.5 ∗ 10−6
2.2 ∗ 10−11
0.9
1.7 ∗ 10−6
2.1 ∗ 10−11
1.4 ∗ 10−6
1.9 ∗ 10−11Exact and approximate solutions of Example3.
(a)
(b)Example 4.
Consider the fractional diffusion equation presented in [18]:
(54)
D
t
α
u
(
x
,
t
)
=
D
x
2
u
(
x
,
t
)
+
sin
(
π
x
)
(
π
2
t
2
+
2
t
2
-
α
Γ
(
3
-
α
)
)
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
l
x
∈
(
0,1
)
,
t
∈
(
0,1
)
,
0
<
α
<
1
,
u
(
0
,
t
)
=
u
(
1
,
t
)
=
0
,
0
<
t
<
1
,
u
(
x
,
0
)
=
0
,
-
1
<
x
<
1
.
The exact solution is
(55)
u
(
x
,
t
)
=
sin
(
π
x
)
t
2
.
To apply the proposed method, we will do the following change of variable X
=
2
x
-
1. In this case, the x-domain becomes [
-
1,1
]. To make a comparison with the results of [18], assume that e
I, e
II, and e
III are the errors in [18] using uniform mesh, quasiuniform mesh, and nonuniform mesh for T
=
1. Let e
pro be the error in the proposed method for T
=
1 and N
=
M
=
10. Results are presented in Table 4.Table 4
Results of the method in [18] and the proposed method.
α
Number of mesh points in [18]
e
I
e
II
e
III
Er
0.25
10
3.5
*
1
0
-
4
2.4
*
1
0
-
4
1.9
*
1
0
-
4
2.3
*
1
0
-
11
20
1.1
*
1
0
-
4
7.2
*
1
0
-
5
5.0
*
1
0
-
5
40
3.5
*
1
0
-
5
2.2
*
1
0
-
5
1.3
*
1
0
-
5
80
1.1
*
1
0
-
5
6.5
*
1
0
-
6
3.5
*
1
0
-
6
0.5
10
1.3
*
1
0
-
3
8.9
*
1
0
-
4
6.1
*
1
0
-
4
4.1
*
1
0
-
12
20
4.8
*
1
0
-
4
3.1
*
1
0
-
4
1.9
*
1
0
-
4
40
1.7
*
1
0
-
4
1.1
*
1
0
-
4
5.8
*
1
0
-
5
80
6.1
*
1
0
-
5
3.8
*
1
0
-
5
1.9
*
1
0
-
5
0.75
10
3.8
*
1
0
-
3
2.7
*
1
0
-
3
1.9
*
1
0
-
3
6.7
*
1
0
-
11
20
1.6
*
1
0
-
3
1.1
*
1
0
-
3
7.4
*
1
0
-
4
40
6.9
*
1
0
-
4
4.7
*
1
0
-
4
2.9
*
1
0
-
4
80
2.9
*
1
0
-
4
2.0
*
1
0
-
4
1.2
*
1
0
-
4


## 5. Conclusion

In this paper, we use series expansion based on the shifted fractional Legendre functions to solve fractional diffusions equations of Caputo’s type. We write the coefficients of the fractional derivative in terms of the shifted fractional Legendre functions as indicated in Theorem5 and give explicit relationship between them. Then, we use the collocation method to compute these coefficients. To the best of our knowledge, the method has not been developed to integrate fractional diffusion equations of the form (1)-(2). We test the proposed technique for several examples and present four of them in this paper. These examples show the efficiency and the accuracy of the proposed method, where in few terms we achieved accuracy up to 1
0
-
10. In Examples 3 and 4, we compare our results with the ones obtained by FDM in [17, 18]. Both examples show that the proposed method works more efficiently and accurately than the methods in [17, 18].



---

*Source: 381074-2014-07-24.xml*