# Two Hybrid Methods for Solving Two-Dimensional Linear Time-Fractional Partial Differential Equations

Authors: B. A. Jacobs; C. Harley
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A computationally efficient hybridization of the Laplace transform with two spatial discretization techniques is investigated for numerical solutions of time-fractional linear partial differential equations in two space variables. The Chebyshev collocation method is compared with the standard finite difference spatial discretization and the absolute error is obtained for several test problems. Accurate numerical solutions are achieved in the Chebyshev collocation method subject to both Dirichlet and Neumann boundary conditions. The solution obtained by these hybrid methods allows for the evaluation at any point in time without the need for time-marching to a particular point in time.


---

## Body

## 1. Introduction

In recent years, fractional derivatives and fractional partial differential equations (FPDEs) have received great attention both in analysis and application (see [1–3] and references therein). In spite of this, there has been very little work done on solving FPDEs on a bounded domain. Agrawal [4] makes use of the Laplace transform and the finite sine transform to obtain an analytic solution to the fractional diffusion-wave equation on a bounded domain. Many authors have applied He’s variational iterative method (VIM) [5] to FPDEs with great success. However, like the differential transform method [6] (DTM) and Adomian decomposition method [7] (ADM), VIM assumes the FPDE to lie on an infinite domain.A great deal of work on hybrid techniques has been presented where the Laplace transform is coupled with Legendre wavelets, ADM, VIM [8–12]. Yin et al. [13] also couple the VIM with Legendre wavelets for the solution of nonlinear PDEs. Song et al. [14] compare the results obtained by the fractional VIM and ADM concluding that the methodologies are similar; however, the VIM does not require the calculation of Adomian polynomials. Atangana and Kılıçman [15] couple the Homotopy perturbation method (HPM) with the Sumudu transform, an integral transform similar to the Laplace transform, to obtain solutions to certain nonlinear heat-like FPDEs. Bhrawy et al. have contributed a considerable effort to the field, focusing on solutions to fractional differential equations via the shifted Legendre tau method and spectral methods for FPDEs [16–20]. Recently, there has been a large movement toward new techniques for solving FPDEs (see [21–28]). To the best of the authors’ knowledge, these transform techniques are unable to enforce Dirichlet or Neumann boundary conditions on a bounded domain, and as such we investigate a new methodology to attempt to circumvent this.The application of diffusion equations to images is abundant [29–34]. However, the application of a time-fractional partial differential equation has not yet been thoroughly examined. Preserving the spatial topography of an image is imperative to maintaining what is deemed to be useful information in that image. It is with this application in mind that we focus on the diffusion equation in two dimensions, where the introduction of an advection term into the model would propagate information. Similarly, considering the time-fractional diffusion equation with 1
<
α
<
2, we obtain the diffusion-wave equation, which again has propagational properties. We therefore restrict our choice of α to be on [
0,1
]. The introduction of a fractional order derivative raises the question of how to discretize or transform the derivative to produce a form that is amenable to existing techniques.The Grünwald-Letnikov discretization has been used for numerical schemes for fractional partial differential equations [35–37] and is given by
(1)
D
t
α
u
(
t
)
=
lim
⁡
Δ
t
→
0
1
Δ
t
α
∑
k
=
0
j
‍
ω
k
α
u
(
t
-
k
Δ
t
)
,
where
(2)
ω
k
α
=
(
-
1
)
k
(
α
k
)
,
and M
Δ
t
=
t. Computationally, this discretization becomes extremely expensive for long time simulations as each subsequent step in time is dependant onevery time step that has preceded it.The use of Laplace transform allows one to circumvent the problems that arise in the time-domain discretization. However, using the Laplace transform for the fractional order derivative presents the problem of inverting the transform to find a solution. Analytic inversion of the transform is infeasible and hence the numerical scheme for evaluating the Bromwich integral presented by Weideman and Trefethen in [38] is put to extensive use. The Bromwich integral is of the form
(3)
u
(
t
)
=
1
2
π
i
∫
σ
-
i
∞
σ
+
i
∞
‍
e
s
t
U
(
s
)
d
s
,
σ
>
σ
0
,
where σ
0 is the convergence abscissa. Using the parabolic conformal map
(4)
s
=
μ
(
i
w
+
1
)
2
,
-
∞
<
w
<
∞
,
equation (3) becomes
(5)
u
(
t
)
=
1
2
π
i
∫
-
∞
∞
‍
e
s
(
w
)
t
U
(
s
(
w
)
)
s
′
(
w
)
d
w
which can then be approximated simply by the trapezoidal rule, or any other quadrature technique, as
(6)
u
h
;
N
t
(
t
)
=
h
2
π
i
∑
k
=
-
N
t
N
t
‍
e
s
(
w
k
)
t
U
(
s
(
w
k
)
)
s
′
(
w
k
)
,
w
k
=
k
h
.
On the parabolic contour the exponential factor in (3) forces a rapid decay in the integrand making it amenable to quadrature. All that is left is to choose the parameters h and μ optimally which is done by asymptotically balancing the truncation error and discretization errors in each of the half-planes. The methodology described by Weideman and Trefethen achieves near optimal results (see [38] and figures therein) and the interested reader is directed there for a thorough description of the implementation of the method. In this work, we make use of the parabolic contour due to the ease of use and the hyperbolic contour only exhibits a slight improvement in performance over the parabolic contour.We present here an extension to the work conducted by Jacobs and Harley in [39]. This paper extends the aforementioned work to the general form of a time-fractional parabolic partial differential equation as well as including two types of boundary conditions, Dirichlet and Neumann.The following section introduces some preliminary definitions followed by a description of the methods employed, including the different cases for boundary conditions in Section3. Section 5 presents the results for comparison based on three fundamentally different examples of linear FPDEs. A discussion of the results and their relationship to work beyond this research is presented in Section 6 as well as some concluding remarks.


## 2. Preliminaries

In this work we employ Caputo’s definition of a fractional derivative over the Riemann-Louiville derivative due to the fact that the Caputo derivative makes use of the physical boundary conditions, whereas the Riemann-Louiville derivative requires fractional order boundary conditions.Definition 1.
The Riemann-Louiville integral of orderα
>
0 of a function u
(
t
) is
(7)
J
α
u
(
t
)
=
1
Γ
(
α
)
∫
0
t
‍
(
t
-
τ
)
α
-
1
u
(
τ
)
d
τ
,
x
>
0
.Definition 2.
The fractional derivative ofu
(
t
) according to the Caputo definition with m
-
1
<
α
≤
m, m
∈
ℕ, is
(8)
∂
α
u
(
t
)
∂
t
α
=
J
m
-
α
D
m
u
(
t
)
=
D
*
α
u
(
t
)
:
=
1
Γ
(
m
-
α
)
∫
0
t
‍
(
t
-
τ
)
m
-
α
-
1
u
(
m
)
(
τ
)
d
τ
.Ifα
=
1, the Caputo fractional derivative reduces to the ordinary first order derivative. Podlubny [2] illustrates the pleasing property of the Laplace Transform of a Caputo derivative, as can be seen in (9). In our case, where 0
<
α
<
1, we have
(9)
ℒ
{
∂
α
u
(
x
,
y
,
t
)
∂
t
α
}
=
s
α
U
(
x
,
y
)
-
s
α
-
1
u
(
x
,
y
,
0
)
.
This property allows one to treat fractional order derivatives algebraically.Definition 3.
The generalized Mittag-Leffler function of the argumentz is
(10)
E
α
,
β
(
z
)
=
∑
k
=
0
∞
‍
z
k
Γ
(
k
α
+
β
)
.


## 3. Methods: Semidiscrete Hybrid Transform Method

This section introduces the methodologies used for a two-dimensional FPDE, where the one-dimensional case is a simple reduction of the methods presented here.Consider the time-fractional differential equation of the form(11)
∂
α
u
∂
t
α
=
L
(
u
,
u
x
,
u
x
x
,
u
y
,
u
y
y
)
,
(
x
,
y
)
∈
Ω
⊂
ℛ
2
with
(12)
u
(
x
,
y
,
0
)
=
f
(
x
,
y
)
,
where L is a linear function of its arguments, Ω
=
[
-
1,1
]
×
[
-
1,1
] to satisfy the domain required by the Chebyshev polynomials, and f
(
x
,
y
) is a functional representation of our image data or a multivariable function. The boundary conditions may be taken to be Dirichlet or Neumann and will be discussed later. We may now apply a Laplace transform to (11) to obtain
(13)
s
α
U
(
x
,
y
)
-
s
α
-
1
f
(
x
,
y
)
=
ℒ
{
L
(
u
,
u
x
,
u
x
x
,
u
y
,
u
y
y
)
}
,
(14)
s
α
U
(
x
,
y
)
-
s
α
-
1
f
(
x
,
y
)
=
L
(
U
,
U
x
,
U
x
x
,
U
y
,
U
y
y
)
,
where
(15)
ℒ
{
u
(
x
,
y
,
t
)
}
=
U
(
x
,
y
,
s
)
.
Boundary conditions may be in the form of Dirichlet conditions as follows:
(16)
u
(
-
1
,
y
,
t
)
=
a
,
u
(
1
,
y
,
t
)
=
b
,
u
(
x
,
-
1
,
t
)
=
c
,
u
(
x
,
1
,
t
)
=
d
,
and hence,
(17)
U
(
-
1
,
y
)
=
ℒ
{
a
}
,
U
(
1
,
y
)
=
ℒ
{
b
}
,
(18)
U
(
x
,
-
1
)
=
ℒ
{
c
}
,
U
(
x
,
1
)
=
ℒ
{
d
}
.
Alternatively, Neumann boundary conditions give
(19)
u
x
(
-
1
,
y
,
t
)
=
a
,
u
x
(
1
,
y
,
t
)
=
b
,
u
y
(
x
,
-
1
,
t
)
=
c
,
u
y
(
x
,
1
,
t
)
=
d
,
with
(20)
U
x
(
-
1
,
y
)
=
ℒ
{
a
}
,
U
y
(
1
,
y
)
=
ℒ
{
b
}
,
U
y
(
x
,
-
1
)
=
ℒ
{
c
}
,
U
y
(
x
,
1
)
=
ℒ
{
d
}
.
The parameters a, b, c, and d are potentially the functions of the temporal variable and one of the spatial variables; that is, a
=
a
(
y
,
t
). Without loss of generality, we however assume that a, b, c, and d are constant.The spatial components of this model are discretized in two ways: by Chebyshev collocation and by finite differences.


### 3.1. Chebyshev Collocation

Chebyshev polynomials form a basis on[
-
1,1
] and hence we dictate the domain of our PDE to be Ω. We note here, however, that any domain in ℝ
2 can be trivially deformed to match Ω. We discretize our spatial domain using Chebyshev-Gauss-Lobatto points:
(21)
x
i
=
cos
⁡
(
i
π
N
x
)
,
i
=
0,1
,
…
,
N
x
,
y
j
=
cos
⁡
(
j
π
N
y
)
,
j
=
0,1
,
…
,
N
y
.
Note here that x
0
=
1, x
N
x
=
-
1, y
0
=
1, and y
N
y
=
-
1 indicating that the domain is in essence reversed and one must use caution when imposing the boundary conditions.Given that our input function or image has been mapped toΩ, we may assume that N
x
=
N
y; that is, we have equal number of collocation points in each spatial direction. We now define a differentiation matrix D
(
1
)
=
d
k
l:
(22)
d
k
l
=
{
c
k
(
-
1
)
k
+
l
c
l
(
x
k
-
x
l
)
,
k
≠
l
,
-
x
k
2
(
1
-
x
k
2
)
,
k
=
l
,
1
6
(
2
N
x
2
+
1
)
,
k
=
l
=
0
,
-
1
6
(
2
N
x
2
+
1
)
,
k
=
l
=
N
x
,
where
c
k
=
{
2
,
k
=
0
,
N
x
,
1
,
k
=
1
,
…
,
N
x
-
1
.
Bayliss et al. [40] describe a method for minimizing the round-off errors incurred in the calculations of the differentiation matrix. Since we write D
(
2
)
=
D
(
1
)
·
D
(
1
), we implement the method, described in [40], in order to minimize propagation of round-off errors for the second derivative in space.The derivative matrices in thex direction are
(23)
D
^
x
(
1
)
=
D
(
1
)
,
D
^
x
(
2
)
=
D
(
2
)
,
where D
(
1
) is the Chebyshev differentiation matrix of size (
N
x
+
1
)
×
(
N
y
+
1
).Because we have assumedN
x
=
N
y, we derive the pleasing property that D
^
y
(
1
)
=
(
D
^
x
(
1
)
)
T and D
^
y
(
2
)
=
(
D
^
x
(
2
)
)
T.Writing the discretization of (14) in matrix form yields
(24)
D
^
x
(
2
)
U
+
U
D
^
y
(
2
)
-
s
α
U
=
-
F
,
where
(25)
F
i
j
=
s
α
-
1
f
(
x
i
,
y
j
)
.
By expanding (24) in summation notation, we have
(26)
∑
k
=
0
N
x
‍
d
i
k
(
2
)
U
(
x
k
,
y
j
)
+
∑
k
=
0
N
y
‍
d
k
j
(
2
)
U
(
x
i
,
y
k
)
-
s
α
U
(
x
i
,
y
j
)
=
-
s
α
-
1
f
(
x
i
,
y
j
)
,
for i
=
0,1
,
…
,
N
x, j
=
0,1
,
…
,
N
y. By extracting the first and last terms in sums, we obtain
(27)
d
i
0
(
2
)
U
(
x
0
,
y
j
)
+
d
i
N
x
(
2
)
U
(
x
N
x
,
y
j
)
+
d
0
j
(
2
)
U
(
x
i
,
y
0
)
+
d
N
y
j
(
2
)
U
(
x
i
,
y
N
y
)
+
∑
k
=
1
N
x
-
1
‍
d
i
k
(
2
)
U
(
x
k
,
y
j
)
+
∑
k
=
1
N
y
-
1
‍
d
k
j
(
2
)
U
(
x
i
,
y
k
)
-
s
α
U
(
x
i
,
y
j
)
=
-
s
α
-
1
f
(
x
i
,
y
j
)
,
for i
=
1
,
…
,
N
x
-
1, j
=
1
,
…
,
N
y
-
1. We use the form of (27) to impose the boundary conditions.The solutionU
~
=
{
U
(
x
1
,
y
1
)
,
U
(
x
1
,
y
2
)
,
…
,
U
(
x
N
x
-
1, y
N
y
-
1
)
}, which is the matrix of unknown interior points of U, can be obtained by solving the system
(28)
(
A
~
-
s
α
I
)
U
~
+
U
~
B
~
=
-
F
~
,
where A
~ is the matrix of interior points of D
^
x
(
2
) and B
~ is the matrix of interior points of D
^
y
(
2
), so that A
~ and B
~ match the dimensions of U
~. Also
(29)
F
~
i
j
=
s
α
-
1
f
(
x
i
,
y
j
)
+
d
i
0
(
2
)
U
(
x
0
,
y
j
)
+
d
i
N
x
(
2
)
U
(
x
N
x
,
y
j
)
+
d
0
j
(
2
)
U
(
x
i
,
y
0
)
+
d
N
y
j
(
2
)
U
(
x
i
,
y
N
y
)
,
for i
=
1
,
…
,
N
x
-
1, j
=
1
,
…
,
N
y
-
1. By using Kronecker tensor products, denoted by ⊗, and a lexicographic reordering, or reshaping, of U
~ and -
F
~ we may write this as
(30)
(
(
A
~
-
s
α
I
)
⊗
I
+
I
⊗
B
~
T
)
U
~
=
-
F
~
,
which is a linear system that is readily solved by  LinearSolve  in  Mathematica 9.


#### 3.1.1. Dirichlet Boundary Conditions

Dirichlet boundary conditions can be imposed directly by substituting (16) into (27) and collecting all the known terms in F
~.


#### 3.1.2. Neumann Boundary Conditions

Neumann boundary conditions given by (17) are discretized as
(31)
∂
U
∂
x
(
x
N
x
=
-
1
,
y
j
)
≈
∑
k
=
0
N
x
‍
d
N
x
k
U
(
x
k
,
y
j
)
=
ℒ
{
a
}
,
∂
U
∂
x
(
x
0
=
1
,
y
j
)
≈
∑
k
=
0
N
x
‍
d
0
k
U
(
x
k
,
y
j
)
=
ℒ
{
b
}
.
Similarly for (18), we have
(32)
∂
U
∂
y
(
x
i
,
y
N
y
=
-
1
)
≈
∑
k
=
0
N
y
‍
d
k
N
y
U
(
x
i
,
y
k
)
=
ℒ
{
c
}
,
∂
U
∂
y
(
x
i
,
y
0
=
1
)
≈
∑
k
=
0
N
y
‍
d
k
0
U
(
x
i
,
y
k
)
=
ℒ
{
d
}
.
By extracting the first and last terms in the sum, the discretizations can be written as
(33)
(
d
N
x
0
d
N
x
N
x
d
00
d
0
N
x
)
(
U
(
x
0
,
y
j
)
U
(
x
N
x
,
y
j
)
)
=
(
ℒ
{
a
}
-
∑
k
=
1
N
x
-
1
‍
d
N
x
k
U
(
x
k
,
y
j
)
ℒ
{
b
}
-
∑
k
=
1
N
x
-
1
‍
d
0
k
U
(
x
k
,
y
j
)
)
,
(
d
0
N
y
d
N
y
N
y
d
00
d
N
y
0
)
(
U
(
x
i
,
y
0
)
U
(
x
i
,
y
N
y
)
)
=
(
ℒ
{
c
}
-
∑
k
=
1
N
y
-
1
‍
d
k
N
y
U
(
x
i
,
y
k
)
ℒ
{
d
}
-
∑
k
=
1
N
y
-
1
‍
d
k
0
U
(
x
i
,
y
k
)
)
.
The solutions to these linear systems are then substituted into (27).


### 3.2. Finite Difference Discretization

Below we make use of the following finite difference formulae:(34)
∂
U
∂
x
≈
1
2
Δ
x
(
U
i
+
1
,
j
-
U
i
-
1
,
j
)
,
∂
U
∂
y
≈
1
2
Δ
y
(
U
i
,
j
+
1
-
U
i
,
j
-
1
)
,
∂
2
U
∂
x
2
≈
1
Δ
x
2
(
U
i
+
1
,
j
-
2
U
i
,
j
+
U
i
-
1
,
j
)
,
∂
2
U
∂
y
2
≈
1
Δ
y
2
(
U
i
,
j
+
1
-
2
U
i
,
j
+
U
i
,
-
1
j
)
,
where Δ
x
=
2
/
N
x and Δ
y
=
2
/
N
y.


#### 3.2.1. Dirichlet Boundary Conditions

Discretizing (14) using a standard central-difference scheme and writing in matrix notation, we deduce
(35)
(
C
~
-
s
α
I
)
U
~
+
U
~
D
~
=
-
G
~
,
where C
~
and D
~ are tridiagonal matrices corresponding to the finite difference differential matrix with dimension (
N
x
-
1
,
N
y
-
1
) and
(36)
G
~
i
j
=
s
α
-
1
f
(
x
i
,
y
j
)
,
i
=
1
,
…
,
N
x
-
1
,
j
=
1
,
…
,
N
y
-
1
.
We write this as a linear system to be solved, again using LinearSolve in Mathematica 9, as follows:
(37)
(
(
C
~
-
s
α
I
)
⊗
I
+
I
⊗
D
~
T
)
U
~
=
-
G
~
.
The boundary conditions, (16), can be enforced directly onto the matrix U, where the interior points of U are U
~.


#### 3.2.2. Neumann Boundary Conditions

By discretizing the boundary conditions (17) and (18) using a standard central-difference scheme we derive
(38)
U
-
1
,
j
=
U
1
,
j
-
2
Δ
x
ℒ
{
a
}
,
U
N
x
+
1
,
j
=
U
N
x
-
1
,
j
+
2
Δ
x
ℒ
{
b
}
,
U
i
,
-
1
=
U
i
,
1
-
2
Δ
y
ℒ
{
c
}
,
U
i
,
N
y
+
1
=
U
i
,
N
y
-
1
+
2
Δ
y
ℒ
{
d
}
.
Including the above conditions in the matrices C and D each with dimension (
N
x
+
1
,
N
y
+
1
) we can write the entire system as
(39)
(
C
-
s
α
I
)
U
+
U
D
=
-
G
,
where
(40)
G
i
j
=
s
α
-
1
f
(
x
i
,
y
j
)
,
i
=
0
,
…
,
N
x
,
j
=
0
,
…
,
N
y
.
This may be solved as a linear system by writing
(41)
(
(
C
-
s
α
I
)
⊗
I
+
I
⊗
D
T
)
U
=
-
G
.


## 3.1. Chebyshev Collocation

Chebyshev polynomials form a basis on[
-
1,1
] and hence we dictate the domain of our PDE to be Ω. We note here, however, that any domain in ℝ
2 can be trivially deformed to match Ω. We discretize our spatial domain using Chebyshev-Gauss-Lobatto points:
(21)
x
i
=
cos
⁡
(
i
π
N
x
)
,
i
=
0,1
,
…
,
N
x
,
y
j
=
cos
⁡
(
j
π
N
y
)
,
j
=
0,1
,
…
,
N
y
.
Note here that x
0
=
1, x
N
x
=
-
1, y
0
=
1, and y
N
y
=
-
1 indicating that the domain is in essence reversed and one must use caution when imposing the boundary conditions.Given that our input function or image has been mapped toΩ, we may assume that N
x
=
N
y; that is, we have equal number of collocation points in each spatial direction. We now define a differentiation matrix D
(
1
)
=
d
k
l:
(22)
d
k
l
=
{
c
k
(
-
1
)
k
+
l
c
l
(
x
k
-
x
l
)
,
k
≠
l
,
-
x
k
2
(
1
-
x
k
2
)
,
k
=
l
,
1
6
(
2
N
x
2
+
1
)
,
k
=
l
=
0
,
-
1
6
(
2
N
x
2
+
1
)
,
k
=
l
=
N
x
,
where
c
k
=
{
2
,
k
=
0
,
N
x
,
1
,
k
=
1
,
…
,
N
x
-
1
.
Bayliss et al. [40] describe a method for minimizing the round-off errors incurred in the calculations of the differentiation matrix. Since we write D
(
2
)
=
D
(
1
)
·
D
(
1
), we implement the method, described in [40], in order to minimize propagation of round-off errors for the second derivative in space.The derivative matrices in thex direction are
(23)
D
^
x
(
1
)
=
D
(
1
)
,
D
^
x
(
2
)
=
D
(
2
)
,
where D
(
1
) is the Chebyshev differentiation matrix of size (
N
x
+
1
)
×
(
N
y
+
1
).Because we have assumedN
x
=
N
y, we derive the pleasing property that D
^
y
(
1
)
=
(
D
^
x
(
1
)
)
T and D
^
y
(
2
)
=
(
D
^
x
(
2
)
)
T.Writing the discretization of (14) in matrix form yields
(24)
D
^
x
(
2
)
U
+
U
D
^
y
(
2
)
-
s
α
U
=
-
F
,
where
(25)
F
i
j
=
s
α
-
1
f
(
x
i
,
y
j
)
.
By expanding (24) in summation notation, we have
(26)
∑
k
=
0
N
x
‍
d
i
k
(
2
)
U
(
x
k
,
y
j
)
+
∑
k
=
0
N
y
‍
d
k
j
(
2
)
U
(
x
i
,
y
k
)
-
s
α
U
(
x
i
,
y
j
)
=
-
s
α
-
1
f
(
x
i
,
y
j
)
,
for i
=
0,1
,
…
,
N
x, j
=
0,1
,
…
,
N
y. By extracting the first and last terms in sums, we obtain
(27)
d
i
0
(
2
)
U
(
x
0
,
y
j
)
+
d
i
N
x
(
2
)
U
(
x
N
x
,
y
j
)
+
d
0
j
(
2
)
U
(
x
i
,
y
0
)
+
d
N
y
j
(
2
)
U
(
x
i
,
y
N
y
)
+
∑
k
=
1
N
x
-
1
‍
d
i
k
(
2
)
U
(
x
k
,
y
j
)
+
∑
k
=
1
N
y
-
1
‍
d
k
j
(
2
)
U
(
x
i
,
y
k
)
-
s
α
U
(
x
i
,
y
j
)
=
-
s
α
-
1
f
(
x
i
,
y
j
)
,
for i
=
1
,
…
,
N
x
-
1, j
=
1
,
…
,
N
y
-
1. We use the form of (27) to impose the boundary conditions.The solutionU
~
=
{
U
(
x
1
,
y
1
)
,
U
(
x
1
,
y
2
)
,
…
,
U
(
x
N
x
-
1, y
N
y
-
1
)
}, which is the matrix of unknown interior points of U, can be obtained by solving the system
(28)
(
A
~
-
s
α
I
)
U
~
+
U
~
B
~
=
-
F
~
,
where A
~ is the matrix of interior points of D
^
x
(
2
) and B
~ is the matrix of interior points of D
^
y
(
2
), so that A
~ and B
~ match the dimensions of U
~. Also
(29)
F
~
i
j
=
s
α
-
1
f
(
x
i
,
y
j
)
+
d
i
0
(
2
)
U
(
x
0
,
y
j
)
+
d
i
N
x
(
2
)
U
(
x
N
x
,
y
j
)
+
d
0
j
(
2
)
U
(
x
i
,
y
0
)
+
d
N
y
j
(
2
)
U
(
x
i
,
y
N
y
)
,
for i
=
1
,
…
,
N
x
-
1, j
=
1
,
…
,
N
y
-
1. By using Kronecker tensor products, denoted by ⊗, and a lexicographic reordering, or reshaping, of U
~ and -
F
~ we may write this as
(30)
(
(
A
~
-
s
α
I
)
⊗
I
+
I
⊗
B
~
T
)
U
~
=
-
F
~
,
which is a linear system that is readily solved by  LinearSolve  in  Mathematica 9.


### 3.1.1. Dirichlet Boundary Conditions

Dirichlet boundary conditions can be imposed directly by substituting (16) into (27) and collecting all the known terms in F
~.


### 3.1.2. Neumann Boundary Conditions

Neumann boundary conditions given by (17) are discretized as
(31)
∂
U
∂
x
(
x
N
x
=
-
1
,
y
j
)
≈
∑
k
=
0
N
x
‍
d
N
x
k
U
(
x
k
,
y
j
)
=
ℒ
{
a
}
,
∂
U
∂
x
(
x
0
=
1
,
y
j
)
≈
∑
k
=
0
N
x
‍
d
0
k
U
(
x
k
,
y
j
)
=
ℒ
{
b
}
.
Similarly for (18), we have
(32)
∂
U
∂
y
(
x
i
,
y
N
y
=
-
1
)
≈
∑
k
=
0
N
y
‍
d
k
N
y
U
(
x
i
,
y
k
)
=
ℒ
{
c
}
,
∂
U
∂
y
(
x
i
,
y
0
=
1
)
≈
∑
k
=
0
N
y
‍
d
k
0
U
(
x
i
,
y
k
)
=
ℒ
{
d
}
.
By extracting the first and last terms in the sum, the discretizations can be written as
(33)
(
d
N
x
0
d
N
x
N
x
d
00
d
0
N
x
)
(
U
(
x
0
,
y
j
)
U
(
x
N
x
,
y
j
)
)
=
(
ℒ
{
a
}
-
∑
k
=
1
N
x
-
1
‍
d
N
x
k
U
(
x
k
,
y
j
)
ℒ
{
b
}
-
∑
k
=
1
N
x
-
1
‍
d
0
k
U
(
x
k
,
y
j
)
)
,
(
d
0
N
y
d
N
y
N
y
d
00
d
N
y
0
)
(
U
(
x
i
,
y
0
)
U
(
x
i
,
y
N
y
)
)
=
(
ℒ
{
c
}
-
∑
k
=
1
N
y
-
1
‍
d
k
N
y
U
(
x
i
,
y
k
)
ℒ
{
d
}
-
∑
k
=
1
N
y
-
1
‍
d
k
0
U
(
x
i
,
y
k
)
)
.
The solutions to these linear systems are then substituted into (27).


## 3.1.1. Dirichlet Boundary Conditions

Dirichlet boundary conditions can be imposed directly by substituting (16) into (27) and collecting all the known terms in F
~.


## 3.1.2. Neumann Boundary Conditions

Neumann boundary conditions given by (17) are discretized as
(31)
∂
U
∂
x
(
x
N
x
=
-
1
,
y
j
)
≈
∑
k
=
0
N
x
‍
d
N
x
k
U
(
x
k
,
y
j
)
=
ℒ
{
a
}
,
∂
U
∂
x
(
x
0
=
1
,
y
j
)
≈
∑
k
=
0
N
x
‍
d
0
k
U
(
x
k
,
y
j
)
=
ℒ
{
b
}
.
Similarly for (18), we have
(32)
∂
U
∂
y
(
x
i
,
y
N
y
=
-
1
)
≈
∑
k
=
0
N
y
‍
d
k
N
y
U
(
x
i
,
y
k
)
=
ℒ
{
c
}
,
∂
U
∂
y
(
x
i
,
y
0
=
1
)
≈
∑
k
=
0
N
y
‍
d
k
0
U
(
x
i
,
y
k
)
=
ℒ
{
d
}
.
By extracting the first and last terms in the sum, the discretizations can be written as
(33)
(
d
N
x
0
d
N
x
N
x
d
00
d
0
N
x
)
(
U
(
x
0
,
y
j
)
U
(
x
N
x
,
y
j
)
)
=
(
ℒ
{
a
}
-
∑
k
=
1
N
x
-
1
‍
d
N
x
k
U
(
x
k
,
y
j
)
ℒ
{
b
}
-
∑
k
=
1
N
x
-
1
‍
d
0
k
U
(
x
k
,
y
j
)
)
,
(
d
0
N
y
d
N
y
N
y
d
00
d
N
y
0
)
(
U
(
x
i
,
y
0
)
U
(
x
i
,
y
N
y
)
)
=
(
ℒ
{
c
}
-
∑
k
=
1
N
y
-
1
‍
d
k
N
y
U
(
x
i
,
y
k
)
ℒ
{
d
}
-
∑
k
=
1
N
y
-
1
‍
d
k
0
U
(
x
i
,
y
k
)
)
.
The solutions to these linear systems are then substituted into (27).


## 3.2. Finite Difference Discretization

Below we make use of the following finite difference formulae:(34)
∂
U
∂
x
≈
1
2
Δ
x
(
U
i
+
1
,
j
-
U
i
-
1
,
j
)
,
∂
U
∂
y
≈
1
2
Δ
y
(
U
i
,
j
+
1
-
U
i
,
j
-
1
)
,
∂
2
U
∂
x
2
≈
1
Δ
x
2
(
U
i
+
1
,
j
-
2
U
i
,
j
+
U
i
-
1
,
j
)
,
∂
2
U
∂
y
2
≈
1
Δ
y
2
(
U
i
,
j
+
1
-
2
U
i
,
j
+
U
i
,
-
1
j
)
,
where Δ
x
=
2
/
N
x and Δ
y
=
2
/
N
y.


### 3.2.1. Dirichlet Boundary Conditions

Discretizing (14) using a standard central-difference scheme and writing in matrix notation, we deduce
(35)
(
C
~
-
s
α
I
)
U
~
+
U
~
D
~
=
-
G
~
,
where C
~
and D
~ are tridiagonal matrices corresponding to the finite difference differential matrix with dimension (
N
x
-
1
,
N
y
-
1
) and
(36)
G
~
i
j
=
s
α
-
1
f
(
x
i
,
y
j
)
,
i
=
1
,
…
,
N
x
-
1
,
j
=
1
,
…
,
N
y
-
1
.
We write this as a linear system to be solved, again using LinearSolve in Mathematica 9, as follows:
(37)
(
(
C
~
-
s
α
I
)
⊗
I
+
I
⊗
D
~
T
)
U
~
=
-
G
~
.
The boundary conditions, (16), can be enforced directly onto the matrix U, where the interior points of U are U
~.


### 3.2.2. Neumann Boundary Conditions

By discretizing the boundary conditions (17) and (18) using a standard central-difference scheme we derive
(38)
U
-
1
,
j
=
U
1
,
j
-
2
Δ
x
ℒ
{
a
}
,
U
N
x
+
1
,
j
=
U
N
x
-
1
,
j
+
2
Δ
x
ℒ
{
b
}
,
U
i
,
-
1
=
U
i
,
1
-
2
Δ
y
ℒ
{
c
}
,
U
i
,
N
y
+
1
=
U
i
,
N
y
-
1
+
2
Δ
y
ℒ
{
d
}
.
Including the above conditions in the matrices C and D each with dimension (
N
x
+
1
,
N
y
+
1
) we can write the entire system as
(39)
(
C
-
s
α
I
)
U
+
U
D
=
-
G
,
where
(40)
G
i
j
=
s
α
-
1
f
(
x
i
,
y
j
)
,
i
=
0
,
…
,
N
x
,
j
=
0
,
…
,
N
y
.
This may be solved as a linear system by writing
(41)
(
(
C
-
s
α
I
)
⊗
I
+
I
⊗
D
T
)
U
=
-
G
.


## 3.2.1. Dirichlet Boundary Conditions

Discretizing (14) using a standard central-difference scheme and writing in matrix notation, we deduce
(35)
(
C
~
-
s
α
I
)
U
~
+
U
~
D
~
=
-
G
~
,
where C
~
and D
~ are tridiagonal matrices corresponding to the finite difference differential matrix with dimension (
N
x
-
1
,
N
y
-
1
) and
(36)
G
~
i
j
=
s
α
-
1
f
(
x
i
,
y
j
)
,
i
=
1
,
…
,
N
x
-
1
,
j
=
1
,
…
,
N
y
-
1
.
We write this as a linear system to be solved, again using LinearSolve in Mathematica 9, as follows:
(37)
(
(
C
~
-
s
α
I
)
⊗
I
+
I
⊗
D
~
T
)
U
~
=
-
G
~
.
The boundary conditions, (16), can be enforced directly onto the matrix U, where the interior points of U are U
~.


## 3.2.2. Neumann Boundary Conditions

By discretizing the boundary conditions (17) and (18) using a standard central-difference scheme we derive
(38)
U
-
1
,
j
=
U
1
,
j
-
2
Δ
x
ℒ
{
a
}
,
U
N
x
+
1
,
j
=
U
N
x
-
1
,
j
+
2
Δ
x
ℒ
{
b
}
,
U
i
,
-
1
=
U
i
,
1
-
2
Δ
y
ℒ
{
c
}
,
U
i
,
N
y
+
1
=
U
i
,
N
y
-
1
+
2
Δ
y
ℒ
{
d
}
.
Including the above conditions in the matrices C and D each with dimension (
N
x
+
1
,
N
y
+
1
) we can write the entire system as
(39)
(
C
-
s
α
I
)
U
+
U
D
=
-
G
,
where
(40)
G
i
j
=
s
α
-
1
f
(
x
i
,
y
j
)
,
i
=
0
,
…
,
N
x
,
j
=
0
,
…
,
N
y
.
This may be solved as a linear system by writing
(41)
(
(
C
-
s
α
I
)
⊗
I
+
I
⊗
D
T
)
U
=
-
G
.


## 4. Analysis

### 4.1. Solvability

We have reduced all four cases to a system of linear systems, presented in (30), (37), and (41). As mentioned in the previous sections, we implemented theMathematica solver LinearSolve, which analyzes the matrix structure and adaptively selects the most appropriate method of solution. We have compared a number of standard solution methods, such as SOR, without finding an algorithm that is faster thanLinearSolve. A system M
x
=
b has a unique solution, x
=
M
-
1
b, provided the matrix M has an inverse M
-
1 that exists. In the following analysis we denote the length in one dimension of the respective matrix M by L, since this length is scheme dependent.Proposition 4.
IfM is an irreducible diagonally dominant matrix for which |
m
i
i
|
>
∑
k
=
1
k
≠
i
L
|
m
i
k
| for at least one i, then M is invertible [41].The proof of the above proposition may be found in [41]. All that is left is to show that our matrix M is always irreducible and satisfies Proposition 4 for some i. A useful characterization of an irreducible matrix is as follows:given a system
M
x
=
b
, the matrix
M
is irreducible if a change in any component of
b
will cause a change in the solution
x
.Proof.
LetM
x
=
b and consider M
x
~
=
b
+
ϵ; then, M
(
x
-
x
~
)
=
ϵ. Assuming x
-
x
~
=
0
, then M
(
x
-
x
~
)
=
0
, a contradiction. So x
≠
x
~ and any change in b will result in a change in x. Hence M is irreducible.


#### 4.1.1. Finite Difference Scheme

The structure of the matrix on the right-hand side of (37) and (41) reduces our diagonal dominance criterion
(42)
|
m
i
i
|
>
∑
k
=
1
k
≠
i
L
‍
|
m
i
k
|
to
(43)
|
4
+
s
α
|
>
2
.
In approximating the Bromwich integral by the trapezoidal rule, we truncate the limits of integration from (
-
∞
,
∞
) to [
-
N
t
,
N
t
]. The contour path of the integral is defined to be a parabola, denoted as s above, in the complex plane with a minimum proportional to the truncation parameter N
t and inversely proportional to the final time we are integrating to, t
1. This means that provided the ratio N
t
/
t
1 is sufficiently large the parabola will traverse the complex space avoiding the lower bound required by the condition (43). Hence, the finite difference scheme is solvable provided one chooses parameters for the evaluation of the Bromwich integral that satisfy the above condition. An example of the above condition is illustrated in Figure 1 where the quarter-circle of radius 2 represents the lower bound required by the above condition.Figure 1
Integration contour in the real-imaginary plane nonintersecting with the lower bound forα
=
0.75.


#### 4.1.2. Chebyshev Collocation

In the case of Chebyshev collocation the solvability criterion is a little more difficult to satisfy, and hence we are only able to derive a necessary criterion that our choice of inversion parameters must satisfy. Once again the structure of (30) dictates that the condition
(44)
|
m
i
i
|
>
∑
k
=
1
k
≠
i
L
‍
|
m
i
k
|
reduces to
(45)
|
2
d
i
i
(
2
)
-
s
α
|
>
∑
k
=
1
k
≠
i
L
‍
2
|
d
i
k
|
.
The above condition is easy to verify in practice; however, due to the dependence of the derivative matrix on the resolution parameters N
x and N
y, it is difficult to write a general condition in closed form. Appropriate choices of parameters in the inversion resolution ensure that the scheme is solvable as described in the previous section. Given that the derivative matrix for Chebyshev collocation is full, the right-hand side of the solvability condition (45) produces a lower bound with a much larger magnitude than in the finite difference case; moreover, the lower bound grows with N
x. Examples of contours are given for N
x
=
5 in Figure 2 and for N
x
=
20 in Figure 3 illustrating the adherence of a well-constructed contour to the solvability condition (45).Figure 2
Integration contour in the real-imaginary plane nonintersecting with the lower bound forα
=
0.75 and N
x
=
5.Figure 3
Integration contour in the real-imaginary plane nonintersecting with the lower bound forα
=
0.75 and N
x
=
20.


### 4.2. Accuracy

#### 4.2.1. Finite Difference Scheme

Theoretically, the accuracy of the method is well known to be𝒪
(
Δ
x
2
)
+
𝒪
(
Δ
y
2
) [41]. The errors obtained in Section 5 concur with theoretical bounds.


#### 4.2.2. Chebyshev Collocation

The work of Breuer and Everson [42] and Don and Solomonoff [43] both present arguments on the accuracy of Chebyshev collocation. In practice, the accuracy of the method is measured by numerically differentiating a function and comparing the numerical derivative with the analytic result. The results obtained in the current work are consistent with the results presented in [42, 43] obtaining very good errors for small values of N
x and N
y. The errors in Chebyshev collocation tend to increase rather drastically for very large values of N
x contradicting the typical rule of thumb. The aforementioned work explains that while the truncation error decreases as resolution increases, the round-off errors accumulate dramatically and dominate.


## 4.1. Solvability

We have reduced all four cases to a system of linear systems, presented in (30), (37), and (41). As mentioned in the previous sections, we implemented theMathematica solver LinearSolve, which analyzes the matrix structure and adaptively selects the most appropriate method of solution. We have compared a number of standard solution methods, such as SOR, without finding an algorithm that is faster thanLinearSolve. A system M
x
=
b has a unique solution, x
=
M
-
1
b, provided the matrix M has an inverse M
-
1 that exists. In the following analysis we denote the length in one dimension of the respective matrix M by L, since this length is scheme dependent.Proposition 4.
IfM is an irreducible diagonally dominant matrix for which |
m
i
i
|
>
∑
k
=
1
k
≠
i
L
|
m
i
k
| for at least one i, then M is invertible [41].The proof of the above proposition may be found in [41]. All that is left is to show that our matrix M is always irreducible and satisfies Proposition 4 for some i. A useful characterization of an irreducible matrix is as follows:given a system
M
x
=
b
, the matrix
M
is irreducible if a change in any component of
b
will cause a change in the solution
x
.Proof.
LetM
x
=
b and consider M
x
~
=
b
+
ϵ; then, M
(
x
-
x
~
)
=
ϵ. Assuming x
-
x
~
=
0
, then M
(
x
-
x
~
)
=
0
, a contradiction. So x
≠
x
~ and any change in b will result in a change in x. Hence M is irreducible.


### 4.1.1. Finite Difference Scheme

The structure of the matrix on the right-hand side of (37) and (41) reduces our diagonal dominance criterion
(42)
|
m
i
i
|
>
∑
k
=
1
k
≠
i
L
‍
|
m
i
k
|
to
(43)
|
4
+
s
α
|
>
2
.
In approximating the Bromwich integral by the trapezoidal rule, we truncate the limits of integration from (
-
∞
,
∞
) to [
-
N
t
,
N
t
]. The contour path of the integral is defined to be a parabola, denoted as s above, in the complex plane with a minimum proportional to the truncation parameter N
t and inversely proportional to the final time we are integrating to, t
1. This means that provided the ratio N
t
/
t
1 is sufficiently large the parabola will traverse the complex space avoiding the lower bound required by the condition (43). Hence, the finite difference scheme is solvable provided one chooses parameters for the evaluation of the Bromwich integral that satisfy the above condition. An example of the above condition is illustrated in Figure 1 where the quarter-circle of radius 2 represents the lower bound required by the above condition.Figure 1
Integration contour in the real-imaginary plane nonintersecting with the lower bound forα
=
0.75.


### 4.1.2. Chebyshev Collocation

In the case of Chebyshev collocation the solvability criterion is a little more difficult to satisfy, and hence we are only able to derive a necessary criterion that our choice of inversion parameters must satisfy. Once again the structure of (30) dictates that the condition
(44)
|
m
i
i
|
>
∑
k
=
1
k
≠
i
L
‍
|
m
i
k
|
reduces to
(45)
|
2
d
i
i
(
2
)
-
s
α
|
>
∑
k
=
1
k
≠
i
L
‍
2
|
d
i
k
|
.
The above condition is easy to verify in practice; however, due to the dependence of the derivative matrix on the resolution parameters N
x and N
y, it is difficult to write a general condition in closed form. Appropriate choices of parameters in the inversion resolution ensure that the scheme is solvable as described in the previous section. Given that the derivative matrix for Chebyshev collocation is full, the right-hand side of the solvability condition (45) produces a lower bound with a much larger magnitude than in the finite difference case; moreover, the lower bound grows with N
x. Examples of contours are given for N
x
=
5 in Figure 2 and for N
x
=
20 in Figure 3 illustrating the adherence of a well-constructed contour to the solvability condition (45).Figure 2
Integration contour in the real-imaginary plane nonintersecting with the lower bound forα
=
0.75 and N
x
=
5.Figure 3
Integration contour in the real-imaginary plane nonintersecting with the lower bound forα
=
0.75 and N
x
=
20.


## 4.1.1. Finite Difference Scheme

The structure of the matrix on the right-hand side of (37) and (41) reduces our diagonal dominance criterion
(42)
|
m
i
i
|
>
∑
k
=
1
k
≠
i
L
‍
|
m
i
k
|
to
(43)
|
4
+
s
α
|
>
2
.
In approximating the Bromwich integral by the trapezoidal rule, we truncate the limits of integration from (
-
∞
,
∞
) to [
-
N
t
,
N
t
]. The contour path of the integral is defined to be a parabola, denoted as s above, in the complex plane with a minimum proportional to the truncation parameter N
t and inversely proportional to the final time we are integrating to, t
1. This means that provided the ratio N
t
/
t
1 is sufficiently large the parabola will traverse the complex space avoiding the lower bound required by the condition (43). Hence, the finite difference scheme is solvable provided one chooses parameters for the evaluation of the Bromwich integral that satisfy the above condition. An example of the above condition is illustrated in Figure 1 where the quarter-circle of radius 2 represents the lower bound required by the above condition.Figure 1
Integration contour in the real-imaginary plane nonintersecting with the lower bound forα
=
0.75.


## 4.1.2. Chebyshev Collocation

In the case of Chebyshev collocation the solvability criterion is a little more difficult to satisfy, and hence we are only able to derive a necessary criterion that our choice of inversion parameters must satisfy. Once again the structure of (30) dictates that the condition
(44)
|
m
i
i
|
>
∑
k
=
1
k
≠
i
L
‍
|
m
i
k
|
reduces to
(45)
|
2
d
i
i
(
2
)
-
s
α
|
>
∑
k
=
1
k
≠
i
L
‍
2
|
d
i
k
|
.
The above condition is easy to verify in practice; however, due to the dependence of the derivative matrix on the resolution parameters N
x and N
y, it is difficult to write a general condition in closed form. Appropriate choices of parameters in the inversion resolution ensure that the scheme is solvable as described in the previous section. Given that the derivative matrix for Chebyshev collocation is full, the right-hand side of the solvability condition (45) produces a lower bound with a much larger magnitude than in the finite difference case; moreover, the lower bound grows with N
x. Examples of contours are given for N
x
=
5 in Figure 2 and for N
x
=
20 in Figure 3 illustrating the adherence of a well-constructed contour to the solvability condition (45).Figure 2
Integration contour in the real-imaginary plane nonintersecting with the lower bound forα
=
0.75 and N
x
=
5.Figure 3
Integration contour in the real-imaginary plane nonintersecting with the lower bound forα
=
0.75 and N
x
=
20.


## 4.2. Accuracy

### 4.2.1. Finite Difference Scheme

Theoretically, the accuracy of the method is well known to be𝒪
(
Δ
x
2
)
+
𝒪
(
Δ
y
2
) [41]. The errors obtained in Section 5 concur with theoretical bounds.


### 4.2.2. Chebyshev Collocation

The work of Breuer and Everson [42] and Don and Solomonoff [43] both present arguments on the accuracy of Chebyshev collocation. In practice, the accuracy of the method is measured by numerically differentiating a function and comparing the numerical derivative with the analytic result. The results obtained in the current work are consistent with the results presented in [42, 43] obtaining very good errors for small values of N
x and N
y. The errors in Chebyshev collocation tend to increase rather drastically for very large values of N
x contradicting the typical rule of thumb. The aforementioned work explains that while the truncation error decreases as resolution increases, the round-off errors accumulate dramatically and dominate.


## 4.2.1. Finite Difference Scheme

Theoretically, the accuracy of the method is well known to be𝒪
(
Δ
x
2
)
+
𝒪
(
Δ
y
2
) [41]. The errors obtained in Section 5 concur with theoretical bounds.


## 4.2.2. Chebyshev Collocation

The work of Breuer and Everson [42] and Don and Solomonoff [43] both present arguments on the accuracy of Chebyshev collocation. In practice, the accuracy of the method is measured by numerically differentiating a function and comparing the numerical derivative with the analytic result. The results obtained in the current work are consistent with the results presented in [42, 43] obtaining very good errors for small values of N
x and N
y. The errors in Chebyshev collocation tend to increase rather drastically for very large values of N
x contradicting the typical rule of thumb. The aforementioned work explains that while the truncation error decreases as resolution increases, the round-off errors accumulate dramatically and dominate.


## 5. Results

Example 5 (one-dimensional time-fractional diffusion equation with homogenous Neumann boundary conditions).
We consider first the time-fractional diffusion equation in one dimension:(46)
∂
α
u
∂
t
α
=
∂
2
u
∂
x
2
,
subject to the boundary conditions
(47)
u
x
(
-
1
,
t
)
=
0
,
u
x
(
1
,
t
)
=
0
and initial condition
(48)
u
(
x
,
0
)
=
cos
⁡
(
π
2
(
x
+
1
)
)
.
The results obtained by finite differences and Chebyshev collocation were compared to the exact solution given by Kazemi and Erjaee [44] as
(49)
u
(
x
,
t
)
=
cos
⁡
(
π
2
(
x
+
1
)
)
(
1
-
t
α
E
α
,
1
+
α
(
-
t
α
)
)
,
where E
α
,
β
(
z
) is the generalized Mittag-Leffler function of the argument z. We select α
=
0.8 in line with [44]; however, experimental results indicate that the methods are robust for virtually any value of 0
<
α
≤
1. We note here that the domain was originally [
0
,
π
] and hence a linear transformation in the spatial variables is required to map the domain to [
-
1,1
], which is in accordance with the domain of the Chebyshev polynomials. The errors in the finite difference and Chebyshev schemes are tabulated in Table 1.
Figures4 and 5 illustrate the diminishing error incurred in the process of inverting the Laplace transform. Inversion of the Laplace transform, even in the analytic case, can lead to a singularity at t
=
0. This arises in the numerical inversion and results in a relatively large error in the present schemes near t
=
0. However, this error diminishes rapidly and our schemes obtain an accurate solution.Table 1
Maximum absolute error in the presented method's solution of the problem described by Example5 at time t
=
0
.
5.
N
x
Error in finite difference
Error in Chebyshev collocation
5
0.009927020087821425
0.0009403962041890646
10
0.002471078174960284
2.5800430680789077
×
10
-
8
15
0.0010985435173019864
9.547918011776346
×
10
-
15
20
0.0006179865325516287
8.992806499463768
×
10
-
15Figure 4
Log plot illustrating the diminishing average error in the Chebyshev collocation scheme with time for Example5.Figure 5
Log plot illustrating the diminishing average error in the finite difference scheme with time for Example5.Example 6 (diffusion-advection equation with Dirichlet boundary conditions).
This example considers the time-fractional diffusion-advection equation in one dimension:(50)
∂
α
u
∂
t
α
=
∂
2
u
∂
x
2
+
∂
u
∂
x
,
subject to
(51)
u
(
-
1
,
t
)
=
e
-
1
,
u
(
1
,
t
)
=
e
-
1
,
u
(
x
,
0
)
=
e
-
x
2
.
To the authors’ knowledge, no exact solution exists for the time-fractional diffusion-advection equation. Comparing the present methods withNDSolve in Mathematica 9yields satisfactory results for α
=
1, but no solution can be found for fractional α. We instead compare our solutions in the Laplace domain, where we obtain an exact solution to the transformed equation usingDSolve in Mathematica 9. This allows one to compare the performance of the present methods for various values for α. The errors obtained, for α
=
0.7, are presented in Table 2. Given that these errors are valid in the transform domain we note that the numerical error of 𝒪
(
8.1
2
-
N
t
) is incurred upon inversion of the Laplace transform, as presented by Weideman and Trefethen in [38] and where N
t is typically 50.Table 2
Maximum absolute error in the presented method's solution of the problem described by Example6 in the Laplace domain at s
=
50.
N
x
Error in finite difference
Error in Chebyshev collocation
5
0.0002846591036554467
0.000031320115677323235
10
0.0003424383347704913
5.981351275385904
×
10
-
9
15
0.00028950142395197005
5.763522471780025
×
10
-
12
20
0.0002476113730374256
5.538704044016907
×
10
-
12Example 7 (two-dimensional time-fractional diffusion equation with homogenous Dirichlet boundary conditions).
We now consider the two-dimensional time-fractional diffusion equation(52)
∂
α
u
∂
t
α
=
∂
2
u
∂
x
2
+
∂
2
u
∂
y
2
,
with boundary conditions
(53)
u
(
-
1
,
y
,
t
)
=
0
,
u
(
1
,
y
,
t
)
=
0
,
u
(
x
,
-
1
,
t
)
=
0
,
u
(
x
,
1
,
t
)
=
0
and initial condition
(54)
u
(
x
,
y
,
0
)
=
sin
(
π
(
x
+
1
)
)
sin
(
π
(
y
+
1
)
)
,
so that the boundary conditions are consistent with the initial condition. The parameter α is taken to be 0.8. Momani [45] gives an exact solution as
(55)
u
(
x
,
y
,
t
)
=
sin
(
π
(
x
+
1
)
)
sin
(
π
(
y
+
1
)
)
E
α
(
-
2
t
α
)
,
where E
α
(
z
)
=
E
α
,
1
(
z
) is the Mittag-Leffler function of order α. The efficacy of these methods for fractional order derivatives is illustrated in Table 3.Table 3
Maximum absolute error in the presented method's solution of the problem described by Example7 at time t
=
0
.
5.
N
x
Error in finite difference
Error in Chebyshev collocation
5
0.038120869340055985
0.0064580077147820825
10
0.009341572386976915
8.103172060291985
×
10
-
7
15
0.004521845750043552
3.937461467984349
×
10
-
12
20
0.002567981542938025
6.517009154549669
×
10
-
13


## 6. Concluding Remarks

The results above strongly advocate the use of Chebyshev collocation as a spatial discretization method given the rapid error reduction with increasing spatial resolution. These hybrid methods present a robust way in which one can solve linear time-fractional partial differential equations on a bounded domain with Neumann or Dirichlet boundary conditions, particularly given discrete initial data.Chebyshev collocation presents extremely small errors when compared to the exact solution. We use numerical experiments as substantiation of the method for applying discrete initial conditions where an exact solution may not exist. Figure6 illustrates the decreasing error with increasing number of collocation points for the examples presented in the previous section.Figure 6
Log plot illustrating the increasing accuracy with increasing collocation points for all examples.The efficiency of the Laplace transform within the context of this hybridized method over a time-marching scheme is threefold. First, we are able to treat the fractional order derivative algebraically. The error incurred in the temporal dimension is only attributed to the evaluation of the Bromwich integral and, furthermore, this error drops off rapidly with increased resolution as illustrated in [38]. Finally the solution obtained is a function of time so that we may evaluate our solution at any time rather than needing to march to that time. The Grünwald-Letnikov discretization presented in (1) is an example of a time-marching scheme. The computational time required for a long time solution via the Grünwald-Letnikov discretization is enormous, due to the fractionality being dependant on every time step that precedes the current time. Moreover, every time step incurs a truncation error, so that the further the solution marches the greater the error is; contrastingly, the present method’s error diminishes as time evolves. As a counter-point, if one were seeking a solution after a very short time, then a time-marching scheme may be better suited.The method described and implemented above is a direct one and is therefore free of any iterative scheme. Hence the convergence of the scheme is difficult to speak of. Figure6 illustrates how the approximate solution obtained converges to the exact solution with increasing number of collocation points for the present examples. However, due to the well-known phenomenon of the Chebyshev collocation method exhibiting large round-off errors for large number of collocation points, N
x
>
100, due to finite-precision error accumulation [42, 43], the collocation hybrid scheme is not consistent with any given problem for N
x
→
∞ or Δ
x
→
0.This research has presented a numerical experimental comparison between the standard finite difference method and the Chebyshev collocation method as a means of spatial discretization when hybridized with the Laplace transform. These methods enjoy the benefits of an exact transform in temporal variable and furthermore allow one to easily and efficiently deal with a fractional order derivative, at the cost of numerically inverting the Laplace transform.The goal of these methods is to apply a fractional order diffusion equation to an image on a bounded two-dimensional domain. The use of a discretization is therefore unavoidable given that the initial condition may in fact be discrete.The solution to the discretized equations is found by writing a two-dimensional system of sizeN
x
×
N
x and N
y
×
N
y as a one-dimensional system of size N
x
N
y
×
N
x
N
y. While this is more computationally expensive, it does exhibit an elegance in construction. An alternative approach would be to implement an alternating-direction implicit (ADI) type scheme [46], where each dimension is acted on in turn rather than at once.Due to the Laplace transform being a linear operator, this method is not suitable for nonlinear problems, nor is it applicable to FPDEs with both fractional spatial derivatives and fractional temporal derivatives.We have shown a hybrid method combining the Laplace transform and a spatial discretization which can be extremely effective at solving linear FPDEs on a two-dimensional bounded domain with Dirichlet or Neumann boundary conditions.



---

*Source: 757204-2014-04-08.xml*