# A New Legendre Collocation Method for Solving a Two-Dimensional Fractional Diffusion Equation

Authors: A. H. Bhrawy
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A new spectral shifted Legendre Gauss-Lobatto collocation (SL-GL-C) method is
developed and analyzed to solve a class of two-dimensional initial-boundary fractional diffusion
equations with variable coefficients. The method depends basically on the fact that an expansion
in a series of shifted Legendre polynomialsP
L
,
n
(
x
)
P
L
,
m
(
y
), for the function and its space-fractional
derivatives occurring in the partial fractional differential equation (PFDE), is assumed; the expansion
coefficients are then determined by reducing the PFDE with its boundary and initial conditions
to a system of ordinary differential equations (SODEs) for these coefficients. This system may be
solved numerically by using the fourth-order implicit Runge-Kutta (IRK) method. This method,
in contrast to common finite-difference and finite-element methods, has the exponential rate of
convergence for the two spatial discretizations. Numerical examples are presented in the form of
tables and graphs to make comparisons with the results obtained by other methods and with the
exact solutions more easier.


---

## Body

## 1. Introduction

In recent years there has been a high level of interest of employing spectral methods for numerically solving many types of integral and differential equations, due to their ease of applying them for finite and infinite domains [1–8]. The speed of convergence is one of the great advantages of spectral method. Besides, spectral methods have exponential rates of convergence; they also have high level of accuracy. From the overview of spectral approximation to differential equations, the spectral methods have been divided to four types, namely, collocation [9–11], tau [12, 13], Galerkin [14, 15], and Petrov Galerkin [16, 17] methods. The main idea of all versions of spectral methods is to express the spectral solution of the problem as a finite sum of certain basis functions (combination of orthogonal polynomials or functions) and then to choose the coefficients in order to minimize the difference between the exact and numerical solutions as possible. The spectral collocation method is a specific type of spectral methods, which is more applicable and widely used to solve almost types of differential equations [18–21].Fractional differential equations [22–31] model many phenomena in several fields such as biology, viscoelasticity, finance, fluid mechanics, physics, chemistry, and engineering. Most fractional differential equations do not have exact solutions, so approximation and numerical techniques must be used. Here, we focus on the application of SL-GL-C scheme for the numerical solution of the time dependent fractional diffusion equation in two dimensions. This problem models phenomena exhibiting [32–34] anomalous diffusion that cannot be modeled accurately by the classical two-dimensional diffusion equations with variable coefficients. For the previous numerical methods used for time dependent fractional diffusion equation, see, for example, the papers in [35–40].The solutionu of this equation is approximated as u
N
,
M which can be expressed as a finite expansion of shifted Legendre polynomials for the space variables and then evaluate the spatial fractional partial derivatives of the finite expansion of u
N
,
M at the shifted Legendre-Gauss-Lobatto (SL-GL) quadrature points. Substituting these approximations in the underlined PFDE provides a SODEs in time. This system may be solved by IRK scheme. This scheme is one of the suitable schemes for solving SODEs of first order. Since the IRK scheme has excellent stability properties, then it allows us to reduce computational costs. Indeed, this is the first work for employing SL-GL-C scheme to solve two-dimensional fractional diffusion equations.This paper is organized as follows. We present few relevant properties of fractional integration in the Riemann-Liouville and shifted Legendre polynomials in the coming section. Section3 presents the numerical solution of the two-dimensional fractional diffusion equations with initial-boundary conditions. Section 4 is devoted to solve two test examples. Finally, some concluding remarks are given in the last section.


## 2. Preliminaries and Notation

### 2.1. The Fractional Integration in the Riemann-Liouville Sense

There are several definitions of the fractional integration of orderν
>
0 and not necessarily equivalent to each other; see [41, 42]. The most used definition is due to Riemann-Liouville, which is defined as(1)
J
ν
f
(
x
)
=
1
Γ
(
ν
)
∫
0
x
‍
(
x
-
t
)
ν
-
1
f
(
t
)
d
t
,
ν
>
0
,
x
>
0
,
J
0
f
(
x
)
=
f
(
x
)
.
One of the basic properties of the operator J
ν is
(2)
J
ν
x
β
=
Γ
(
β
+
1
)
Γ
(
β
+
1
+
ν
)
x
β
+
ν
.
The Riemann-Liouville fractional derivative of order ν will be denoted by D
ν. The next relation defines the Riemann-Liouville fractional derivative of order ν as
(3)
D
ν
f
(
x
)
=
d
m
d
x
m
(
J
m
-
ν
f
(
x
)
)
,
where m
-
1
<
ν
≤
m
,
m
∈
N, and m is the smallest integer greater than ν.Lemma 1.
Ifm
-
1
<
ν
≤
m
,
m
∈
N, then
(4)
D
ν
J
ν
f
(
x
)
=
f
(
x
)
,
J
ν
D
ν
f
(
x
)
=
f
(
x
)
-
∑
i
=
0
m
-
1
‍
f
(
i
)
(
0
+
)
x
i
i
!
,
x
>
0
.


### 2.2. Properties of Shifted Legendre Polynomials

The well-known Legendre polynomialsP
i
(
x
) are defined on the interval (
-
1,1
). Some properties about the standard Legendre polynomials will be recalled in this section. The Legendre polynomials P
k
(
x
)
(
k
=
0,1
,
…
,
) satisfy the following Rodrigue's formula:
(5)
P
k
(
x
)
=
(
-
1
)
k
2
k
k
!
D
k
(
(
1
-
x
2
)
k
)
.
We recall also that P
k
(
x
) is a polynomial of degree k, and therefore P
k
(
q
)
(
x
) (the qth derivative of P
k
(
x
)) is
(6)
P
k
(
q
)
(
x
)
=
∑
i
=
0
(
k
+
i
=
even
)
k
-
q
‍
C
q
(
k
,
i
)
P
i
(
x
)
,
where
(7)
C
q
(
k
,
i
)
=
2
q
-
1
(
2
i
+
1
)
Γ
(
q
+
k
-
i
2
)
Γ
(
q
+
k
+
i
+
1
2
)
×
(
Γ
(
q
)
Γ
(
2
-
q
+
k
-
i
2
)
Γ
(
3
-
q
+
k
+
i
2
)
)
-
1
.
The Legendre polynomials satisfy the following relations:
(8)
P
0
(
x
)
=
1
,
P
1
(
x
)
=
x
,
P
k
+
2
(
x
)
=
2
k
+
3
k
+
2
x
P
k
+
1
(
x
)
-
k
+
1
k
+
2
P
k
(
x
)
,
and the orthogonality relation
(9)
(
P
k
(
x
)
,
P
l
(
x
)
)
=
∫
-
1
1
P
k
(
x
)
P
l
(
x
)
w
(
x
)
d
x
=
h
k
δ
l
k
‍
,
where w
(
x
)
≡
1, h
k
=
2
/
(
2
k
+
1
), and δ
l
k is the Kronecker delta function. The Legendre-Gauss-Lobatto quadrature is used to evaluate the previous integral accurately. For any ϕ
∈
S
2
N
-
1
[
-
1,1
], we get
(10)
∫
-
1
1
ϕ
(
x
)
d
x
‍
=
∑
j
=
0
N
‍
ϖ
N
,
j
ϕ
(
x
N
,
j
)
.
We introduce the following discrete inner product:
(11)
(
u
,
v
)
N
=
∑
j
=
0
N
‍
u
(
x
N
,
j
)
v
(
x
N
,
j
)
ϖ
N
,
j
,
where x
N
,
j (0
≤
j
≤
N) and ϖ
N
,
j (0
≤
j
≤
N) are used as usual nodes and the corresponding Christoffel numbers in the interval [
-
1,1
], respectively. Let the shifted Legendre polynomials P
i
(
(
2
x
/
L
)
-
1
) be denoted by P
L
,
i
(
x
); x
∈
(
0
,
L
). Then P
L
,
i
(
x
) can be obtained with the aid of the following recurrence formula:
(12)
(
i
+
1
)
P
L
,
i
+
1
(
x
)
=
(
2
i
+
1
)
(
2
x
L
-
1
)
P
L
,
i
(
x
)
-
i
P
L
,
i
-
1
(
x
)
,
i
=
1,2
,
…
.
The analytic form of the shifted Legendre polynomials P
L
,
i
(
x
) of degree i is given by
(13)
P
L
,
i
(
x
)
=
∑
k
=
0
i
‍
(
-
1
)
i
+
k
(
i
+
k
)
!
(
i
-
k
)
!
(
k
!
)
2
L
k
x
k
.
The Riemann-Liouville fractional integration of shifted Legendre polynomials P
L
,
i
(
x
) of degree i is given by
(14)
J
ν
P
L
,
i
(
x
)
=
∑
k
=
0
i
(
-
1
)
i
+
k
(
i
+
k
)
!
(
i
-
k
)
!
(
k
!
)
2
L
k
J
ν
x
k
‍
=
∑
k
=
0
i
(
-
1
)
i
+
k
(
i
+
k
)
!
k
!
(
i
-
k
)
!
(
k
!
)
2
L
k
Γ
(
k
+
ν
+
1
)
x
k
+
ν
,
‍
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
i
=
0,1
,
…
,
N
,
D
ν
P
L
,
i
(
x
)
=
∑
k
=
0
i
‍
(
-
1
)
i
+
k
(
i
+
k
)
!
(
i
-
k
)
!
(
k
!
)
2
L
k
J
ν
x
k
=
∑
k
=
0
i
‍
(
-
1
)
i
+
k
(
i
+
k
)
!
Γ
(
k
+
1
)
(
i
-
k
)
!
(
k
!
)
2
L
k
Γ
(
k
+
1
-
ν
)
x
k
-
ν
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
i
=
0,1
,
…
,
N
,
where P
L
,
i
(
0
)
=
(
-
1
)
i. The orthogonality condition is
(15)
∫
0
L
P
L
,
j
(
x
)
P
L
,
k
(
x
)
w
L
(
x
)
d
x
‍
=
h
L
,
k
δ
j
k
,
where w
L
(
x
)
≡
1  and  h
L
,
k
=
L
/
(
2
k
+
1
).The functionu
(
x
), square integrable in (
0
,
L
), may be expressed in terms of shifted Legendre polynomials as
(16)
u
(
x
)
=
∑
j
=
0
∞
‍
c
j
P
L
,
j
(
x
)
,
where the coefficients c
j are given by
(17)
c
j
=
1
h
L
,
j
∫
0
L
u
(
x
)
P
L
,
j
(
x
)
d
x
‍
,
j
=
0,1
,
2
,
…
.
In practice, only the first (
N
+
1
)-terms shifted Legendre polynomials are considered. Hence u
(
x
) can be expressed in the form
(18)
u
N
(
x
)
≃
∑
j
=
0
N
‍
c
j
P
L
,
j
(
x
)
.


## 2.1. The Fractional Integration in the Riemann-Liouville Sense

There are several definitions of the fractional integration of orderν
>
0 and not necessarily equivalent to each other; see [41, 42]. The most used definition is due to Riemann-Liouville, which is defined as(1)
J
ν
f
(
x
)
=
1
Γ
(
ν
)
∫
0
x
‍
(
x
-
t
)
ν
-
1
f
(
t
)
d
t
,
ν
>
0
,
x
>
0
,
J
0
f
(
x
)
=
f
(
x
)
.
One of the basic properties of the operator J
ν is
(2)
J
ν
x
β
=
Γ
(
β
+
1
)
Γ
(
β
+
1
+
ν
)
x
β
+
ν
.
The Riemann-Liouville fractional derivative of order ν will be denoted by D
ν. The next relation defines the Riemann-Liouville fractional derivative of order ν as
(3)
D
ν
f
(
x
)
=
d
m
d
x
m
(
J
m
-
ν
f
(
x
)
)
,
where m
-
1
<
ν
≤
m
,
m
∈
N, and m is the smallest integer greater than ν.Lemma 1.
Ifm
-
1
<
ν
≤
m
,
m
∈
N, then
(4)
D
ν
J
ν
f
(
x
)
=
f
(
x
)
,
J
ν
D
ν
f
(
x
)
=
f
(
x
)
-
∑
i
=
0
m
-
1
‍
f
(
i
)
(
0
+
)
x
i
i
!
,
x
>
0
.


## 2.2. Properties of Shifted Legendre Polynomials

The well-known Legendre polynomialsP
i
(
x
) are defined on the interval (
-
1,1
). Some properties about the standard Legendre polynomials will be recalled in this section. The Legendre polynomials P
k
(
x
)
(
k
=
0,1
,
…
,
) satisfy the following Rodrigue's formula:
(5)
P
k
(
x
)
=
(
-
1
)
k
2
k
k
!
D
k
(
(
1
-
x
2
)
k
)
.
We recall also that P
k
(
x
) is a polynomial of degree k, and therefore P
k
(
q
)
(
x
) (the qth derivative of P
k
(
x
)) is
(6)
P
k
(
q
)
(
x
)
=
∑
i
=
0
(
k
+
i
=
even
)
k
-
q
‍
C
q
(
k
,
i
)
P
i
(
x
)
,
where
(7)
C
q
(
k
,
i
)
=
2
q
-
1
(
2
i
+
1
)
Γ
(
q
+
k
-
i
2
)
Γ
(
q
+
k
+
i
+
1
2
)
×
(
Γ
(
q
)
Γ
(
2
-
q
+
k
-
i
2
)
Γ
(
3
-
q
+
k
+
i
2
)
)
-
1
.
The Legendre polynomials satisfy the following relations:
(8)
P
0
(
x
)
=
1
,
P
1
(
x
)
=
x
,
P
k
+
2
(
x
)
=
2
k
+
3
k
+
2
x
P
k
+
1
(
x
)
-
k
+
1
k
+
2
P
k
(
x
)
,
and the orthogonality relation
(9)
(
P
k
(
x
)
,
P
l
(
x
)
)
=
∫
-
1
1
P
k
(
x
)
P
l
(
x
)
w
(
x
)
d
x
=
h
k
δ
l
k
‍
,
where w
(
x
)
≡
1, h
k
=
2
/
(
2
k
+
1
), and δ
l
k is the Kronecker delta function. The Legendre-Gauss-Lobatto quadrature is used to evaluate the previous integral accurately. For any ϕ
∈
S
2
N
-
1
[
-
1,1
], we get
(10)
∫
-
1
1
ϕ
(
x
)
d
x
‍
=
∑
j
=
0
N
‍
ϖ
N
,
j
ϕ
(
x
N
,
j
)
.
We introduce the following discrete inner product:
(11)
(
u
,
v
)
N
=
∑
j
=
0
N
‍
u
(
x
N
,
j
)
v
(
x
N
,
j
)
ϖ
N
,
j
,
where x
N
,
j (0
≤
j
≤
N) and ϖ
N
,
j (0
≤
j
≤
N) are used as usual nodes and the corresponding Christoffel numbers in the interval [
-
1,1
], respectively. Let the shifted Legendre polynomials P
i
(
(
2
x
/
L
)
-
1
) be denoted by P
L
,
i
(
x
); x
∈
(
0
,
L
). Then P
L
,
i
(
x
) can be obtained with the aid of the following recurrence formula:
(12)
(
i
+
1
)
P
L
,
i
+
1
(
x
)
=
(
2
i
+
1
)
(
2
x
L
-
1
)
P
L
,
i
(
x
)
-
i
P
L
,
i
-
1
(
x
)
,
i
=
1,2
,
…
.
The analytic form of the shifted Legendre polynomials P
L
,
i
(
x
) of degree i is given by
(13)
P
L
,
i
(
x
)
=
∑
k
=
0
i
‍
(
-
1
)
i
+
k
(
i
+
k
)
!
(
i
-
k
)
!
(
k
!
)
2
L
k
x
k
.
The Riemann-Liouville fractional integration of shifted Legendre polynomials P
L
,
i
(
x
) of degree i is given by
(14)
J
ν
P
L
,
i
(
x
)
=
∑
k
=
0
i
(
-
1
)
i
+
k
(
i
+
k
)
!
(
i
-
k
)
!
(
k
!
)
2
L
k
J
ν
x
k
‍
=
∑
k
=
0
i
(
-
1
)
i
+
k
(
i
+
k
)
!
k
!
(
i
-
k
)
!
(
k
!
)
2
L
k
Γ
(
k
+
ν
+
1
)
x
k
+
ν
,
‍
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
i
=
0,1
,
…
,
N
,
D
ν
P
L
,
i
(
x
)
=
∑
k
=
0
i
‍
(
-
1
)
i
+
k
(
i
+
k
)
!
(
i
-
k
)
!
(
k
!
)
2
L
k
J
ν
x
k
=
∑
k
=
0
i
‍
(
-
1
)
i
+
k
(
i
+
k
)
!
Γ
(
k
+
1
)
(
i
-
k
)
!
(
k
!
)
2
L
k
Γ
(
k
+
1
-
ν
)
x
k
-
ν
,
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
i
=
0,1
,
…
,
N
,
where P
L
,
i
(
0
)
=
(
-
1
)
i. The orthogonality condition is
(15)
∫
0
L
P
L
,
j
(
x
)
P
L
,
k
(
x
)
w
L
(
x
)
d
x
‍
=
h
L
,
k
δ
j
k
,
where w
L
(
x
)
≡
1  and  h
L
,
k
=
L
/
(
2
k
+
1
).The functionu
(
x
), square integrable in (
0
,
L
), may be expressed in terms of shifted Legendre polynomials as
(16)
u
(
x
)
=
∑
j
=
0
∞
‍
c
j
P
L
,
j
(
x
)
,
where the coefficients c
j are given by
(17)
c
j
=
1
h
L
,
j
∫
0
L
u
(
x
)
P
L
,
j
(
x
)
d
x
‍
,
j
=
0,1
,
2
,
…
.
In practice, only the first (
N
+
1
)-terms shifted Legendre polynomials are considered. Hence u
(
x
) can be expressed in the form
(18)
u
N
(
x
)
≃
∑
j
=
0
N
‍
c
j
P
L
,
j
(
x
)
.


## 3. Shifted Legendre Spectral Collocation Method

In this section, we derive the SL-GL-C method for solving two-dimensional fractional diffusion equation of the form(19)
∂
v
(
x
,
y
,
t
)
∂
t
=
g
1
(
x
,
y
)
∂
ν
1
v
(
x
,
y
,
t
)
∂
x
ν
1
+
g
2
(
x
,
y
)
∂
ν
2
v
(
x
,
y
,
t
)
∂
y
ν
2
+
f
(
x
,
y
,
t
)
,
(
x
,
y
,
t
)
∈
[
0
,
L
1
]
×
[
0
,
L
2
]
×
(
0
,
T
]
,
which is subject to the initial condition
(20)
v
(
x
,
y
,
0
)
=
g
3
(
x
,
y
)
,
(
x
,
y
)
∈
[
0
,
L
1
]
×
[
0
,
L
2
]
and boundary conditions
(21)
v
(
0
,
y
,
t
)
=
g
4
(
y
,
t
)
,
v
(
L
1
,
y
,
t
)
=
g
5
(
y
,
t
)
,
hhhhhhhhhhhhhhhhhh
(
y
,
t
)
∈
[
0
,
L
2
]
×
[
0
,
T
]
,
v
(
x
,
0
,
t
)
=
g
6
(
x
,
t
)
,
v
(
x
,
L
2
,
t
)
=
g
7
(
x
,
t
)
hhhhhhhhhhhhhhhhhh
(
x
,
t
)
∈
[
0
,
L
1
]
×
[
0
,
T
]
.Denote byx
N
,
j
(
x
L
,
N
,
j
)
,
0
⩽
j
⩽
N, and ϖ
N
,
j
(
ϖ
L
,
N
,
j
)
,
0
≤
i
≤
N, the nodes and Christoffel numbers of the standard (shifted) Legendre-Gauss-Lobatto quadrature on the intervals (
-
1,1
)
,
(
0
,
L
), respectively. Then one can clearly deduce that
(22)
x
L
,
N
,
j
=
L
2
(
x
N
,
j
+
1
)
,
0
≤
j
≤
N
,
ϖ
L
,
N
,
j
=
(
L
2
)
α
+
β
+
1
ϖ
N
,
j
,
0
≤
j
≤
N
,
and if S
N
(
0
,
L
) denotes the set of all polynomials of degree at most N, then it follows for any ϕ
∈
S
2
N
-
1
(
0
,
L
) that
(23)
∫
0
L
‍
w
L
(
x
)
ϕ
(
x
)
d
x
=
(
L
2
)
∫
-
1
1
‍
ϕ
(
L
2
(
x
+
1
)
)
d
x
=
(
L
2
)
∑
j
=
0
N
‍
ϖ
N
,
j
ϕ
(
L
2
(
x
N
,
j
+
1
)
)
=
∑
j
=
0
N
‍
ϖ
L
,
N
,
j
ϕ
(
x
L
,
N
,
j
)
.Here, the SL-GL-C method will be employed to transform the previous two-dimensional fractional diffusion equation into SODEs. Let us expand the dependent variables in the form(24)
v
(
x
,
y
,
t
)
=
∑
i
=
0
N
‍
∑
j
=
0
M
‍
a
i
,
j
(
t
)
P
L
1
,
i
(
x
)
P
L
2
,
j
(
y
)
.
We observe from (15) that the coefficients a
i
,
j
(
t
) can be approximated by
(25)
a
i
,
j
(
t
)
=
1
h
L
1
,
i
h
L
2
,
j
×
∑
l
=
0
N
‍
∑
k
=
0
M
‍
(
P
L
2
,
j
(
y
L
2
,
M
,
k
)
ϖ
M
,
k
P
L
1
,
i
(
x
L
1
,
N
,
l
)
ϖ
N
,
l
)
∑
l
=
0
N
‍
h
h
∑
k
=
0
M
‍
×
v
(
x
L
1
,
N
,
l
,
y
L
2
,
M
,
k
,
t
)
.
The expression (24) can be rewritten as
(26)
v
(
x
,
y
,
t
)
=
∑
i
=
0
N
‍
∑
j
=
0
M
‍
∑
l
=
0
N
‍
∑
k
=
0
M
(
P
L
2
,
j
(
y
L
2
,
M
,
k
)
ϖ
M
,
k
P
L
1
,
i
(
x
L
1
,
N
,
l
)
ϖ
N
,
l
)
h
L
1
,
i
h
L
2
,
j
‍
∑
l
=
0
N
‍
∑
k
=
0
M
‍
∑
l
=
0
N
h
‍
∑
k
=
0
M
‍
×
P
L
1
,
i
(
x
)
P
L
2
,
j
(
y
)
v
l
,
k
(
t
)
,
where v
(
x
L
1
,
N
,
l
,
y
L
2
,
M
,
k
,
t
)
=
v
l
,
k
(
t
). In what follows, the approximation of the fractional spatial partial derivative ∂
ν
1
v
(
x
,
y
,
t
)
/
∂
x
ν
1 with respect to x for v
(
x
,
y
,
t
) can be computed as
(27)
∂
ν
1
v
(
x
,
y
,
t
)
∂
x
ν
1
=
∑
i
=
0
N
‍
∑
j
=
0
M
‍
∑
l
=
0
N
‍
∑
k
=
0
M
‍
(
P
L
2
,
j
(
y
L
2
,
M
,
k
)
ϖ
M
,
k
P
L
1
,
i
(
x
L
1
,
N
,
l
)
ϖ
N
,
l
)
h
L
1
,
i
h
L
2
,
j
∑
l
=
0
N
‍
∑
k
=
0
M
‍
∑
l
=
0
N
‍
h
∑
k
=
0
M
‍
×
D
(
ν
1
)
(
P
L
1
,
i
(
x
)
)
P
L
2
,
j
(
y
)
v
l
,
k
(
t
)
.Also, the fractional spatial partial derivative of fractional order∂
ν
2
v
(
x
,
y
,
t
)
/
∂
y
ν
2 with respect to y of the approximate solution is
(28)
∂
ν
2
v
(
x
,
y
,
t
)
∂
y
ν
2
=
∑
l
=
0
N
‍
∑
k
=
0
M
‍
∑
i
=
0
N
‍
∑
j
=
0
M
(
P
L
2
,
j
(
y
L
2
,
M
,
k
)
ϖ
M
,
k
P
L
1
,
i
(
x
L
1
,
N
,
l
)
ϖ
N
,
l
)
h
L
1
,
i
h
L
2
,
j
‍
∑
l
=
0
N
‍
∑
k
=
0
M
‍
∑
l
=
0
N
‍
h
∑
k
=
0
M
‍
×
P
L
1
,
i
(
x
)
D
(
ν
2
)
(
P
L
2
,
j
(
y
)
)
v
l
,
k
(
t
)
.
The previous fractional derivatives can be computed at the SL-GL interpolation nodes as
(29)
(
∂
ν
1
v
(
x
,
y
,
t
)
∂
x
ν
1
)
y
=
y
L
2
,
M
,
m
x
=
x
L
1
,
N
,
n
=
∑
l
=
0
N
‍
∑
k
=
0
M
‍
ρ
l
,
k
n
,
m
v
l
,
k
(
t
)
,
hhhhl
n
=
0,1
,
…
,
N
,
m
=
0,1
,
…
,
M
,
(
∂
ν
2
v
(
x
,
y
,
t
)
∂
y
ν
2
)
y
=
y
L
2
,
M
,
m
x
=
x
L
1
,
N
,
n
=
∑
l
=
0
N
‍
∑
k
=
0
M
‍
λ
l
,
k
n
,
m
v
l
,
k
(
t
)
,
hhhhh
n
=
0,1
,
…
,
N
,
m
=
0,1
,
…
,
M
,
where
(30)
ρ
l
,
k
n
,
m
=
∑
i
=
0
N
‍
∑
j
=
0
M
‍
(
P
L
2
,
j
(
y
L
2
,
M
,
k
)
ϖ
M
,
k
P
L
1
,
i
(
x
L
1
,
N
,
l
)
ϖ
N
,
l
)
h
i
h
j
∑
i
=
0
N
‍
.
∑
j
=
0
M
‍
×
D
(
ν
1
)
P
L
1
,
i
(
x
L
1
,
N
,
n
)
P
L
2
,
j
(
y
L
2
,
M
,
m
)
,
λ
l
,
k
n
,
m
=
∑
i
=
0
N
‍
∑
j
=
0
M
‍
(
P
L
2
,
j
(
y
L
2
,
M
,
k
)
ϖ
M
,
k
P
L
1
,
i
(
x
L
1
,
N
,
l
)
ϖ
N
,
l
)
h
i
h
j
∑
i
=
0
N
‍
.
∑
j
=
0
M
‍
×
D
(
ν
2
)
P
L
2
,
j
(
y
L
2
,
M
,
m
)
P
L
1
,
i
(
x
L
1
,
N
,
n
)
.
In the proposed SL-GL-C method the residual of (19) is set to zero at (
N
-
1
)
×
(
M
-
1
) of SL-GL points. Moreover, the values of v
0
,
k
(
t
), v
N
,
k
(
t
), v
l
,
0
(
t
), and v
l
,
M
(
t
) are given by
(31)
v
0
,
k
(
t
)
=
g
4
(
y
L
2
,
M
,
k
,
t
)
,
v
N
,
k
(
t
)
=
g
5
(
y
L
2
,
M
,
k
,
t
)
,
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
l
h
h
k
=
0
,
…
,
M
,
v
l
,
0
(
t
)
=
g
6
(
x
L
1
,
N
,
l
,
t
)
,
v
l
,
M
(
t
)
=
g
7
(
x
L
1
,
N
,
l
,
t
)
,
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
l
=
0
,
…
,
N
.
Therefore, adapting (26)–(31) enables one to write (19)–(21) as a (
N
-
1
)
×
(
M
-
1
) SODEs:
(32)
v
˙
n
,
m
(
t
)
=
g
1
n
,
m
∑
l
=
0
N
‍
∑
k
=
0
M
‍
ρ
l
,
k
n
,
m
v
l
,
k
(
t
)
+
g
2
n
,
m
∑
l
=
0
N
‍
∑
k
=
0
M
‍
λ
l
,
k
n
,
m
v
l
,
k
(
t
)
+
f
n
,
m
(
t
)
,
n
=
1
,
…
,
N
-
1
,
m
=
1
,
…
,
M
-
1
,
which is subject to the initial conditions:
(33)
v
n
,
m
(
0
)
=
g
3
n
,
m
,
n
=
1
,
…
,
N
-
1
,
m
=
1
,
…
,
M
-
1
,
where
(34)
g
1
(
x
L
1
,
N
,
n
,
y
L
2
,
M
,
m
)
=
g
1
n
,
m
,
g
2
(
x
L
1
,
N
,
n
,
y
L
2
,
M
,
m
)
=
g
2
n
,
m
,
g
3
(
x
L
1
,
N
,
n
,
y
L
2
,
M
,
m
)
=
g
3
n
,
m
,
f
(
x
L
1
,
N
,
n
,
y
L
2
,
M
,
m
,
t
)
=
f
n
,
m
(
t
)
.
Finally, we can arrange (32) and (33) to their matrix formulation, which can be solved by IRK scheme:
(35)
(
v
˙
1,1
(
t
)
⋯
v
˙
1
,
M
-
1
(
t
)
v
˙
2,1
(
t
)
⋮
v
˙
2
,
M
-
1
(
t
)
⋯
⋱
⋯
⋯
⋱
⋯
⋯
⋱
⋯
v
˙
N
-
2,1
(
t
)
⋮
v
˙
N
-
2
,
M
-
1
(
t
)
v
˙
N
-
1,1
(
t
)
⋯
v
˙
N
-
1
,
M
-
1
(
t
)
)
=
(
H
1,1
(
t
,
v
1
,
…
,
v
N
-
1
)
⋯
H
1
,
M
-
1
(
t
,
v
1
,
…
,
v
N
-
1
)
H
2,1
(
t
,
v
1
,
…
,
v
N
-
1
)
⋮
H
2
,
M
-
1
(
t
,
v
1
,
…
,
v
N
-
1
)
⋯
⋱
⋯
⋯
⋱
⋯
⋯
⋱
⋯
H
N
-
2,1
(
t
,
v
1
,
…
,
v
N
-
1
)
⋮
H
N
-
2
,
M
-
1
(
t
,
v
1
,
…
,
v
N
-
1
)
H
N
-
1,1
(
t
,
v
1
,
…
,
v
N
-
1
)
⋯
H
N
-
1
,
M
-
1
(
t
,
v
1
,
…
,
v
N
-
1
)
)
,
(
v
1,1
(
0
)
⋯
v
1
,
M
-
1
(
0
)
v
2,1
(
0
)
⋮
v
2
,
M
-
1
(
0
)
⋯
⋱
⋯
⋯
⋱
⋯
⋯
⋱
⋯
v
N
-
2,1
(
0
)
⋮
v
N
-
2
,
M
-
1
(
0
)
v
N
-
1,1
(
0
)
⋯
v
N
-
1
,
M
-
1
(
0
)
)
=
(
g
3
1,1
⋯
g
3
1
,
M
-
1
g
3
2,1
⋮
g
3
2
,
M
-
1
⋯
⋱
⋯
⋯
⋱
⋯
⋯
⋱
⋯
g
3
N
-
2,1
⋮
g
3
N
-
2
,
M
-
1
g
3
N
-
1,1
⋯
g
3
N
-
1
,
M
-
1
)
,
t
∈
[
-
τ
,
0
]
,
where
(36)
H
n
,
m
(
t
,
v
1,1
,
…
,
v
N
-
1
,
M
-
1
)
=
g
1
n
,
m
∑
l
=
0
N
‍
∑
k
=
0
M
‍
ρ
l
,
k
n
,
m
v
l
,
k
(
t
)
+
g
2
n
,
m
∑
l
=
0
N
‍
∑
k
=
0
M
‍
λ
l
,
k
n
,
m
v
l
,
k
(
t
)
+
f
n
,
m
(
t
)
,
n
=
1
,
…
,
N
-
1
,
m
=
1
,
…
,
M
-
1
.The Runge-Kutta method can be expressed as one of powerful numerical integrations tools used for initial value SODEs of first order. The IRK method presents a subclass of the well-known family of Runge-Kutta methods and has many applications in the numerical solution of systems of ordinary differential equations.


## 4. Numerical Results

In this section, we present two numerical examples to show the accuracy, robustness, and applicability of the proposed method. We compare the results obtained by our method with those obtained in [43].The difference between the measured value of the approximate solution and the exact solution (absolute error) is given by(37)
E
(
x
,
y
,
t
)
=
|
u
(
x
,
y
,
t
)
-
u
N
(
x
,
y
,
t
)
|
,
where x, y, t, u
(
x
,
y
,
t
), and u
N
(
x
,
y
,
t
) are the space vectors, time, exact, and numerical solutions, respectively. Moreover, the maximum absolute error (MAEs) is given by
(38)
M
E
=
Max
⁡
{
E
(
x
,
y
,
t
)
:
over
all
domains
}
.
Also, we define the norm infinity as
(39)
M
E
t
=
Max
⁡
{
E
(
x
,
y
,
t
k
)
:
over
all
space
domains
}
.Example 1.
Consider the two-dimensional problem [43] is
(40)
∂
u
∂
t
(
x
,
y
,
t
)
=
g
1
(
x
,
t
)
∂
1.8
u
(
x
,
y
,
t
)
∂
x
1.8
+
g
2
(
x
,
t
)
∂
1.6
u
(
x
,
y
,
t
)
∂
x
1.6
+
f
(
x
,
y
,
t
)
,
(
x
,
y
,
t
)
∈
[
0,1
]
×
[
0,1
]
×
[
0
,
T
]
,with the diffusion coefficients(41)
g
1
(
x
,
t
)
=
Γ
(
2.2
)
x
2.8
y
6
,
g
2
(
x
,
t
)
=
2
x
y
2.6
Γ
(
4.6
)
,
and consider the forcing function is
(42)
f
(
x
,
y
,
t
)
=
-
(
1
+
2
x
y
)
e
-
t
x
3
y
3.6
,
with the initial condition
(43)
u
(
x
,
y
,
0
)
=
x
3
y
3.6
and Dirichlet boundary conditions
(44)
u
(
0
,
y
,
t
)
=
u
(
x
,
0
,
t
)
=
0
,
u
(
1
,
y
,
t
)
=
e
-
t
y
3.6
,
u
(
x
,
1
,
t
)
=
e
-
t
x
3
∀
t
≥
0
.
The exact solution of this two-dimensional fractional diffusion equation is given by
(45)
u
(
x
,
y
,
t
)
=
e
-
t
x
3
y
3.6
.From the first look on Table1, we see high accurate results based on MAEs at different choices of N and M, compared with the results given in [43]. Moreover, results listed in Table 2 confirm that the proposed method is very accurate.Table 1
Maximum absolute errors for problem (40) for T
=
1.
N
=
M
SL-GL-C method
[43]
M
E
t
M
E
Δ
x
=
Δ
y
=
Δ
t
M
E
4
3.41785
×
1
0
-
3
9.11427
×
1
0
-
3
—
—
6
2.13706
×
1
0
-
5
6.83861
×
1
0
-
5
1
5
1.05564
×
1
0
-
3
8
1.32866
×
1
0
-
6
4.85909
×
1
0
-
6
1
10
2.61694
×
1
0
-
4
10
1.74023
×
1
0
-
7
7.07142
×
1
0
-
7
1
20
6.57067
×
1
0
-
5
12
3.43827
×
1
0
-
8
1.52415
×
1
0
-
7
1
40
1.55866
×
1
0
-
5Table 2
Maximum absolute errors for problem (40) for T
=
1 and N
=
M
=
12.
x
y
t
E
(
x
,
y
,
t
)
x
y
t
E
(
x
,
y
,
t
)
x
y
t
E
(
x
,
y
,
t
)
0.1
0.1
1
2.74515
×
1
0
-
11
0.1
0.5
1
9.17923
×
1
0
-
12
0.1
1
1
4.41521
×
1
0
-
11
0.2
2.13241
×
1
0
-
10
0.2
7.57278
×
1
0
-
11
0.2
3.58585
×
1
0
-
10
0.3
7.22081
×
1
0
-
10
0.3
2.61927
×
1
0
-
10
0.3
1.20622
×
1
0
-
9
0.4
1.71192
×
1
0
-
9
0.4
6.43406
×
1
0
-
10
0.4
2.93613
×
1
0
-
9
0.5
3.35365
×
1
0
-
9
0.5
1.29893
×
1
0
-
9
0.5
5.8445
×
1
0
-
9
0.6
5.80344
×
1
0
-
9
0.6
2.33485
×
1
0
-
9
0.6
1.02404
×
1
0
-
8
0.7
9.20721
×
1
0
-
9
0.7
3.82976
×
1
0
-
9
0.7
1.65193
×
1
0
-
8
0.8
1.39884
×
1
0
-
8
0.8
5.84963
×
1
0
-
9
0.8
2.50569
×
1
0
-
8
0.9
2.06847
×
1
0
-
8
0.9
8.5766
×
1
0
-
9
0.9
3.64012
×
1
0
-
8
1
2.75218
×
1
0
-
8
1
1.26484
×
1
0
-
8
1
5.60815
×
1
0
-
8In Figure1, we display the numerical solution of problem (40) computed by the present method at N
=
M
=
12, while we plot the space graph of absolute error of problem (40) computed by the present method at N
=
M
=
12 in Figure 2. As shown in Figures 3 and 4, we see the agreement of numerical and exact solutions curves.Figure 1
Space graph of numerical solution of problem (40).Figure 2
Space graph of absolute error of problem (40).Figure 3
x-direction curves of numerical and exact solutions of problem (40).Figure 4
y-direction curves of numerical and exact solutions of problem (40).Example 2.
Consider the two-dimensional fractional diffusion equation is(46)
∂
u
∂
t
(
x
,
y
,
t
)
=
g
1
(
x
,
t
)
∂
ν
1
u
(
x
,
y
,
t
)
∂
x
ν
1
+
g
2
(
x
,
t
)
∂
ν
2
u
(
x
,
y
,
t
)
∂
y
ν
2
+
f
(
x
,
y
,
t
)
,
(
x
,
y
,
t
)
∈
[
0,1
]
×
[
0,1
]
×
[
0
,
T
]
,with the diffusion coefficients(47)
g
1
(
x
,
y
)
=
(
3
-
2
x
)
Γ
(
3
-
ν
1
)
2
,
g
2
(
x
,
y
)
=
(
4
-
y
)
Γ
(
4
-
ν
2
)
6
,
and consider the forcing function is
(48)
f
(
x
,
y
,
t
)
=
e
-
t
(
x
2
(
-
y
3
/
2
+
y
-
4
)
y
3
/
2
+
x
(
2
x
-
3
)
y
3
)
,
with the initial condition
(49)
u
(
x
,
y
,
0
)
=
x
2
y
3
and the Dirichlet boundary conditions
(50)
u
(
0
,
y
,
t
)
=
u
(
x
,
0
,
t
)
=
0
,
u
(
1
,
y
,
t
)
=
e
-
t
y
3
,
u
(
x
,
1
,
t
)
=
e
-
t
x
2
∀
t
≥
0
.
The exact solution of this two-dimensional fractional diffusion equation is given by
(51)
u
(
x
,
y
,
t
)
=
e
-
t
x
2
y
3
.More accurate results for maximum absolute and absolute errors of (46) have been obtained in Tables 3 and 4, respectively, with values of parameters listed in their captions, the approximate solution u
N
,
M at ν
1
=
ν
2
=
1.5 and N
=
M
=
10 is plotted in Figure 5. Moreover, we plot the absolute error for this problem at ν
1
=
ν
2
=
1.5 and N
=
M
=
10 in Figure 6. Finally, the x- and y-directions curves of absolute error at ν
1
=
ν
2
=
1.5 and N
=
M
=
10 are displayed in Figures 7 and 8.Table 3
Maximum absolute errors for problem (46) for T
=
1 and ν
1
=
ν
2
=
1.5.
N
=
M
M
E
t
M
E
6
2.57399
×
1
0
-
8
1.02853
×
1
0
-
7
8
1.51431
×
1
0
-
8
1.21098
×
1
0
-
7
10
8.24341
×
1
0
-
14
4.13336
×
1
0
-
13Table 4
Maximum absolute errors for problem (46) for T
=
1 and ν
1
=
ν
2
=
1.5.
x
y
t
E
(
x
,
y
,
t
)
x
y
t
E
(
x
,
y
,
t
)
x
y
t
E
(
x
,
y
,
t
)
0.1
0.1
1
4.06522
×
1
0
-
16
0.1
0.5
1
7.05989
×
1
0
-
15
0.1
1
1
2.35055
×
1
0
-
15
0.2
3.12066
×
1
0
-
15
0.2
7.445
×
1
0
-
15
0.2
1
6.63792
×
1
0
-
15
0.3
4.22817
×
1
0
-
15
0.3
1.73386
×
1
0
-
15
0.3
1
6.92588
×
1
0
-
15
0.4
8.7365
×
1
0
-
15
0.4
4.25354
×
1
0
-
15
0.4
1
1.70948
×
1
0
-
14
0.5
1.40339
×
1
0
-
14
0.5
1.56351
×
1
0
-
14
0.5
1
2.83974
×
1
0
-
14
0.6
1.81075
×
1
0
-
14
0.6
1.14631
×
1
0
-
14
0.6
1
3.44343
×
1
0
-
14
0.7
2.69246
×
1
0
-
14
0.7
1.8454
×
1
0
-
14
0.7
1
5.44946
×
1
0
-
14
0.8
3.42964
×
1
0
-
14
0.8
3.14054
×
1
0
-
14
0.8
1
6.78416
×
1
0
-
14
0.9
4.3951
×
1
0
-
14
0.9
2.78319
×
1
0
-
14
0.9
1
9.03513
×
1
0
-
14
1
3.1379
×
1
0
-
15
1
1.20095
×
1
0
-
14
1
1
1.0103
×
1
0
-
14Figure 5
Space graph of numerical solution of problem (46) at ν
1
=
ν
2
=
1.5 and N
=
M
=
10.Figure 6
Space graph of absolute error of problem (46) at ν
1
=
ν
2
=
1.5 and N
=
M
=
10.Figure 7
x-direction curve of absolute error of problem (46) at ν
1
=
ν
2
=
1.5 and N
=
M
=
10.Figure 8
y-direction curves of absolute error of problem (46) at ν
1
=
ν
2
=
1.5 and N
=
M
=
10.


## 5. Concluding Remarks and Future Work

The SL-GL-C method was investigated successfully in spatial discretizations to get accurate approximate solution of the two-dimensional fractional diffusion equations with variable coefficients. This problem was transformed into SODEs in time variable which greatly simplifies the problem. The IRK scheme was then applied to the resulting system. From the numerical experiments, the obtained results demonstrated the effectiveness and the high accuracy of SL-GL-C method for solving the mentioned problem.The technique can be extended to more sophisticated problems in two- and three-dimensional spaces. In principle, this method may be extended to related problems in mathematical physics. It is possible to use other orthogonal polynomials, say, Chebyshev polynomials or Jacobi polynomials, to solve the mentioned problem in this paper. Furthermore, the proposed spectral method might be developed by considering the Legendre pseudospectral approximation in both temporal and spatial discretizations. We should note that, as a numerical method, we are restricted to solving problem over a finite domain. Also, the pseudospectral approximation might be employed based on generalized Laguerre or modified generalized Laguerre polynomials [44] to solve similar problems in a semi-infinite spatial intervals.



---

*Source: 636191-2014-05-04.xml*