# A Computational Study of an Implicit Local Discontinuous Galerkin Method for Time-Fractional Diffusion Equations

Authors: Leilei Wei; Xindong Zhang
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We propose, analyze, and test a fully discrete local discontinuous Galerkin (LDG) finite element method for a time-fractional diffusion equation. The proposed method is based on a finite difference scheme in time and local discontinuous Galerkin methods in space. By choosing the numerical fluxes carefully, we prove that our scheme is unconditionally stable and convergent. Finally, numerical examples are performed to illustrate the effectiveness and the accuracy of the method.


---

## Body

## 1. Introduction

Fractional calculus which is considered as the generalization of the integer order calculus attracts much attention recently its of their numerous applications in physics and engineering. They provide an excellent instrument for the description of memory and hereditary properties of various materials and processes. This is the main advantage of fractional derivatives in comparison with classic integral-order models, in which such effects are, in fact, neglected. Interest of some scholars has been shown in research on the problems involving the fractional order partial differential equations (PDEs) [1–19]. Machado et al. [20] introduced the recent history of fractional calculus; as for the detailed theory and applications of fractional integrals and derivatives, we can refer to [21, 22] and the references therein. Due to their numerous applications in the areas of physics and engineering, solving such equations and numerical schemes for fractional differential equations has been stimulated.Fractional equations arise in continuous-time random walks, modeling of anomalous diffusive and subdiffusive systems, unification of diffusion and wave propagation phenomenon, and simplification of the results. There are only a few numerical works in the literature to solve fractional diffusion equations. Liu et al. [23] use a first-order finite difference scheme in both time and space directions for this equation, where some stability conditions are derived. In [14] Lin and Xu examine a practical finite difference/Legendre spectral method to solve the initial-boundary value time-fractional diffusion problem on a finite domain. In [11], Jiang and Ma use high-order finite element methods to solve the equation and prove an optimal convergence rate.In this paper, we consider the following time-fractional diffusion equation:(1)
∂
α
u
(
x
,
t
)
∂
t
α
-
∂
2
u
(
x
,
t
)
∂
x
2
=
f
(
x
,
t
)
,
(
x
,
t
)
∈
[
a
,
b
]
×
[
0
,
T
]
,
u
(
x
,
0
)
=
u
0
(
x
)
,
x
∈
[
a
,
b
]
,
where 0
<
α
<
1 is the order of the time-fractional derivatives. f and u
0 are given smooth functions. We do not pay attention to boundary condition in this paper; hence, the solution is considered to be either periodic or compactly supported.We define(
∂
α
u
(
x
,
t
)
)
/
(
∂
t
α
) as the Caputo fractional derivatives of order α [24, 25],
(2)
∂
α
u
(
x
,
t
)
∂
t
α
=
1
Γ
(
1
-
α
)
∫
0
t
‍
∂
u
(
x
,
s
)
∂
s
d
s
(
t
-
s
)
α
,
0
<
α
<
1
;
here, Γ
(
·
) is the Gamma function.In the present paper, we propose a fully discrete local discontinuous Galerkin (LDG) finite element method for solving the time-fractional diffusion equation. Our fully discrete scheme is based on a finite difference scheme in time and local discontinuous Galerkin methods in space. By choosing the numerical fluxes carefully, we prove that our scheme is unconditionally stable and gives an error estimate.What remains of this paper is organized as follows. We begin by introducing some basic notations and mathematical preliminaries which are required for establishing our results. In Section3, we discuss the LDG scheme for the fractional equation (1), and we prove that the fully discrete scheme is unconditionally stable and convergent. Numerical experiments to illustrate the accuracy and capability of the method are given in Section 4. Finally, in Section 5, concluding remarks are provided.


## 2. Notations and Auxiliary Results

In this section, we introduce notations and definitions to be used later in the paper and also present some auxiliary results.Given a spatial grida
=
x
1
/
2
<
x
3
/
2
<
⋯
<
x
N
+
(
1
/
2
)
=
b, define the mesh I
j
=
[
x
j
-
(
1
/
2
)
,
x
j
+
(
1
/
2
)
], for j
=
1
,
…
N and the cell lengths Δ
x
j
=
x
j
+
(
1
/
2
)
-
x
j
-
(
1
/
2
)
,
1
≤
j
≤
N, and h
=
max
⁡
1
≤
j
≤
N
Δ
x
j.We denote byu
j
+
(
1
/
2
)
+ and u
j
+
(
1
/
2
)
- the values of u at x
j
+
1
/
2, from the right cell I
j
+
1 and from the left cell I
j. [
u
]
j
+
(
1
/
2
) is used to denote u
j
+
(
1
/
2
)
+
-
u
j
+
(
1
/
2
)
-, that is, the jump of u at cell interfaces.We define the piecewise-polynomial spaceV
h
k as the space of polynomials of the degree up to k in each cell I
j; that is,
(3)
V
h
k
=
{
v
:
v
∈
P
k
(
I
j
)
,
x
∈
I
j
,
j
=
1,2
,
…
N
}
.For error estimates, we will be using two projections in one dimension[
a
,
b
], denoted by P; that is, for each j,
(4)
∫
I
j
‍
(
P
ω
(
x
)
-
ω
(
x
)
)
v
(
x
)
=
0
,
∀
v
∈
P
k
(
I
j
)
,
and special projection P
±; that is, for each j,
(5)
∫
I
j
‍
(
P
+
ω
(
x
)
-
ω
(
x
)
)
v
(
x
)
=
0
,
∀
v
∈
P
k
-
1
(
I
j
)
,
P
+
ω
(
x
j
-
(
1
/
2
)
+
)
=
ω
(
x
j
-
(
1
/
2
)
)
,
(6)
∫
I
j
‍
(
P
-
ω
(
x
)
-
ω
(
x
)
)
v
(
x
)
=
0
,
∀
v
∈
P
k
-
1
(
I
j
)
,
P
-
ω
(
x
j
+
(
1
/
2
)
-
)
=
ω
(
x
j
+
(
1
/
2
)
)
.For the two projections, the following inequality holds [26–28]:
(7)
∥
ω
e
∥
+
h
∥
ω
e
∥
∞
+
h
1
/
2
∥
ω
e
∥
τ
h
≤
C
0
h
k
+
1
,
where ω
e
=
P
ω
-
ω or ω
e
=
P
±
ω
-
ω. The positive constant C
0, solely depending on ω, is independent of h. τ
h denotes the set of boundary points of all elements I
j.In the present paper, we useC
i
(
i
=
0,1
,
2
) to denote a positive constant which may have a different value in each occurrence. The usual notation of norms in Sobolev spaces will be used. Let the scalar inner product on L
2
(
D
) be denoted by (
·
,
·
)
D and the associated norm by ∥
·
∥
D. If D
=
Ω, we drop D.


## 3. Fully Discrete LDG Scheme

LetΔ
t
=
T
/
M be the time meshsize, let M be a positive integer, let t
n
=
n
Δ
t
,
n
=
0,1
,
…
,
M be mesh point. First, we estimate the time-fractional derivatives (
∂
α
u
(
x
,
t
)
)
/
(
∂
t
α
) at t
n as follows [11, 14]:
(8)
∂
α
u
(
x
,
t
n
)
∂
t
α
=
(
Δ
t
)
1
-
α
Γ
(
2
-
α
)
×
∑
i
=
0
n
-
1
‍
b
i
u
(
x
,
t
n
-
i
)
-
u
(
x
,
t
n
-
i
-
1
)
Δ
t
+
γ
n
(
x
)
,
where b
i
=
(
i
+
1
)
1
-
α
-
i
1
-
α, γ
n
(
x
)
≤
C
1
(
Δ
t
)
2
-
α, C
1 is dependent on u
,
T
,
α.We rewrite (1) as a first-order system:
(9)
p
=
u
x
,
∂
α
u
(
x
,
t
)
∂
t
α
-
p
x
=
f
,
and we give the weak form of (9) at t
n as follows:
(10)
∫
Ω
‍
u
(
x
,
t
n
)
v
d
x
-
β
λ
1
(
∫
Ω
‍
u
(
x
,
t
n
)
v
x
d
x
-
∑
j
=
1
N
(
(
u
(
x
,
t
n
)
v
-
)
j
+
(
1
/
2
)
‍
-
(
u
(
x
,
t
n
)
v
+
)
j
-
(
1
/
2
)
)
(
∫
Ω
‍
u
(
x
,
t
n
)
v
x
d
x
)
+
β
λ
2
(
∫
Ω
‍
p
(
x
,
t
n
)
v
x
d
x
-
∑
j
=
1
N
‍
(
(
p
(
x
,
t
n
)
v
-
)
j
+
(
1
/
2
)
-
(
p
(
x
,
t
n
)
v
+
)
j
-
(
1
/
2
)
)
(
∫
Ω
‍
p
(
x
,
t
n
)
v
x
d
x
)
=
∑
i
=
1
n
-
1
‍
(
b
i
-
1
-
b
i
)
∫
Ω
‍
u
(
x
,
t
n
-
i
)
v
d
x
+
b
n
-
1
∫
Ω
‍
u
(
x
,
t
0
)
v
d
x
-
β
∫
Ω
‍
γ
n
(
x
)
v
d
x
+
β
∫
Ω
‍
f
(
x
,
t
n
)
v
d
x
,
∫
Ω
‍
p
(
x
,
t
n
)
w
d
x
+
∫
Ω
‍
u
(
x
,
t
n
)
w
x
d
x
-
∑
j
=
1
N
‍
(
(
u
(
x
,
t
n
)
w
-
)
j
+
(
1
/
2
)
-
(
u
(
x
,
t
n
)
w
+
)
j
-
(
1
/
2
)
)
=
0
,
where β
=
(
Δ
t
)
α
Γ
(
2
-
α
).Letu
h
n
,
p
h
n
∈
V
h
k be the approximation of u
(
·
,
t
n
)
,
p
(
·
,
t
n
); respectively, f
n
(
x
)
=
f
(
x
,
t
n
). We define a fully discrete local discontinuous Galerkin scheme as follows: find u
h
n
,
p
h
n
∈
V
h
k, such that, for all test functions v
,
w
∈
V
h
k,
(11)
∫
Ω
‍
u
h
n
v
d
x
+
β
(
-
∑
j
=
1
N
‍
(
(
p
h
n
^
v
-
)
j
+
(
1
/
2
)
-
(
p
h
n
^
v
+
)
j
-
(
1
/
2
)
)
)
∫
Ω
‍
p
h
n
v
x
d
x
-
∑
j
=
1
N
‍
(
(
p
h
n
^
v
-
)
j
+
(
1
/
2
)
-
(
p
h
n
^
v
+
)
j
-
(
1
/
2
)
)
)
=
∑
i
=
1
n
-
1
(
b
i
-
1
-
b
i
)
∫
Ω
‍
u
h
n
-
i
v
d
x
+
b
n
-
1
∫
Ω
‍
u
h
0
v
d
x
+
β
∫
Ω
‍
f
n
v
d
x
,
∫
Ω
‍
p
h
n
w
d
x
+
∫
Ω
‍
u
h
n
w
x
d
x
-
∑
j
=
1
N
‍
(
(
u
h
n
^
w
-
)
j
+
(
1
/
2
)
-
(
u
h
n
^
w
+
)
j
-
(
1
/
2
)
)
=
0
.The “hat” terms in (11) in the cell boundary terms from integration by parts are the so-called “numerical fluxes,” which are single valued functions defined on the edges and should be designed based on different guiding principles for different PDEs to ensure stability. It turns out that we can take the simple choices such that
(12)
u
h
n
^
=
(
u
h
n
)
-
,
p
h
n
^
=
(
p
h
n
)
+
.We remark that the choice for the fluxes (12) is not unique. In fact, the crucial part is taking u
h
n
^ and p
h
n
^ from opposite sides. We know the truncation error is γ
n
(
x
) from (8).In order to simplify the notations and without loss of generality, we consider the casef
=
0 in its numerical analysis. Now, we consider the stability for the scheme (11), we have the following result.Theorem 1.
For periodic or compactly supported boundary conditions, the fully-discrete LDG scheme (11) is unconditionally stable, and the numerical solution u
h
n satisfies
(13)
∥
u
h
n
∥
2
+
2
β
∥
p
h
n
∥
2
≤
∥
u
h
0
∥
2
,
n
=
1,2
…
,
M
.Proof.
We will prove Theorem1 by mathematical induction. When n
=
1, scheme (11) is
(14)
∫
Ω
‍
u
h
1
v
d
x
+
β
(
-
∑
j
=
1
N
(
(
p
h
1
^
v
-
)
j
+
(
1
/
2
)
+
(
p
h
1
^
v
+
)
j
-
(
1
/
2
)
)
)
∫
Ω
‍
p
h
1
v
x
d
x
-
∑
j
=
1
N
(
(
p
h
1
^
v
-
)
j
+
(
1
/
2
)
+
(
p
h
1
^
v
+
)
j
-
(
1
/
2
)
)
)
=
∫
Ω
‍
u
h
0
v
d
x
,
∫
Ω
‍
p
h
1
w
d
x
+
∫
Ω
‍
u
h
1
w
x
d
x
-
∑
j
=
1
N
‍
(
(
u
h
1
^
w
-
)
j
+
(
1
/
2
)
+
(
u
h
1
^
w
+
)
j
-
(
1
/
2
)
)
=
0
.
Taking the test functionsv
=
u
h
1
,
w
=
β
p
h
1, we obtain
(15)
∥
u
h
1
∥
Ω
2
+
β
∥
p
h
1
∥
Ω
2
+
∑
j
=
1
N
‍
β
τ
[
u
h
1
]
j
-
(
1
/
2
)
2
+
∑
j
=
1
N
‍
β
(
Ψ
(
p
h
1
,
u
h
1
)
j
+
(
1
/
2
)
-
Ψ
(
p
h
1
,
u
h
1
)
j
-
(
1
/
2
)
+
Θ
(
p
h
1
,
u
h
1
)
j
-
(
1
/
2
)
)
=
∫
Ω
‍
u
h
0
u
h
1
d
x
.
Here,
(16)
Ψ
(
p
h
1
,
u
h
1
)
=
(
p
h
1
)
-
(
u
h
1
)
-
-
p
h
1
^
(
u
h
1
)
-
+
u
h
1
^
(
p
h
1
)
-
,
Θ
(
p
h
1
,
u
h
1
)
=
(
p
h
1
)
-
(
u
h
1
)
-
-
(
p
h
1
)
+
(
u
h
1
)
+
-
p
h
1
~
(
u
h
1
)
-
+
p
h
1
~
(
u
h
1
)
+
+
u
h
1
^
(
p
h
1
)
-
-
u
h
1
^
(
p
h
1
)
+
.
If we take fluxes (12) and after some manual calculation, we can easily obtain Θ
(
p
h
1
,
u
h
1
)
j
-
(
1
/
2
)
=
0.
From the fact that(17)
∫
Ω
‍
u
h
0
u
h
1
d
x
≤
1
2
∥
u
h
1
∥
Ω
2
+
1
2
∥
u
h
0
∥
Ω
2
,
we can get
(18)
∥
u
h
1
∥
Ω
≤
∥
u
h
0
∥
Ω
.
Now, suppose the following inequality holds(19)
∥
u
h
m
∥
Ω
≤
∥
u
h
0
∥
Ω
,
m
=
1,2
,
…
,
P
.
We need to prove ∥
u
h
P
+
1
∥
Ω
≤
∥
u
h
0
∥
Ω. Let n
=
P
+
1 and take the test functions v
=
u
h
P
+
1
,
w
=
β
p
h
P
+
1 in scheme (11); we can obtain
(20)
∥
u
h
P
+
1
∥
Ω
2
+
β
∥
p
h
P
+
1
∥
Ω
2
+
∑
j
=
1
N
‍
β
τ
[
u
h
P
+
1
]
j
-
(
1
/
2
)
2
+
∑
j
=
1
N
‍
β
(
Ψ
(
p
h
P
+
1
,
u
h
P
+
1
)
j
+
(
1
/
2
)
-
Ψ
(
p
h
P
+
1
,
u
h
P
+
1
)
j
-
(
1
/
2
)
+
Θ
(
p
h
P
+
1
,
u
h
P
+
1
)
j
-
(
1
/
2
)
)
=
∑
i
=
1
P
‍
(
b
i
-
1
-
b
i
)
∫
Ω
‍
u
h
P
+
1
-
i
u
h
P
+
1
d
x
+
b
P
∫
Ω
‍
u
h
0
u
h
P
+
1
d
x
≤
∑
i
=
1
P
‍
(
b
i
-
1
-
b
i
)
∥
u
h
P
+
1
-
i
∥
Ω
∥
u
h
P
+
1
∥
Ω
+
b
P
∥
u
h
0
∥
Ω
∥
u
h
P
+
1
∥
Ω
≤
(
∑
i
=
1
P
‍
(
b
i
-
1
-
b
i
)
+
b
P
)
∥
u
h
0
∥
Ω
∥
u
h
P
+
1
∥
Ω
.
Taking fluxes (12), we can easily obtain Θ
j
-
(
1
/
2
)
(
p
h
P
+
1
,
u
h
P
+
1
)
=
0. Then, the last inequality gives
(21)
∥
u
h
P
+
1
∥
Ω
≤
∥
u
h
0
∥
Ω
.
This finishes the proof of the stability result.Theorem 2.
Letu
(
x
,
t
n
) be the exact solution of the problem (1), which is sufficiently smooth with bounded derivatives. Let u
h
n be the numerical solution of the fully discrete LDG scheme (11); then, there hold the following error estimates when 0
<
α
<
1:
(22)
∥
u
(
x
,
t
n
)
-
u
h
n
∥
≤
C
T
α
1
-
α
(
(
Δ
t
)
-
α
h
k
+
1
+
(
Δ
t
)
2
-
α
+
(
Δ
t
)
-
α
/
2
h
k
+
(
1
/
2
)
+
h
k
+
1
)
and when α
→
1:
(23)
∥
u
(
x
,
t
n
)
-
u
h
n
∥
≤
T
C
(
(
Δ
t
)
-
1
h
k
+
1
+
Δ
t
+
(
Δ
t
)
-
1
/
2
h
k
+
(
1
/
2
)
+
h
k
+
1
)
.Proof.
We denote(24)
e
u
n
=
u
(
x
,
t
n
)
-
u
h
n
=
P
-
e
u
n
-
(
P
-
u
(
x
,
t
n
)
-
u
(
x
,
t
n
)
)
,
e
p
n
=
p
(
x
,
t
n
)
-
p
h
n
=
P
+
e
p
n
-
(
P
+
p
(
x
,
t
n
)
-
p
(
x
,
t
n
)
)
.
Subtracting (11) from (10) and with fluxes (12), we can obtain the error equation
(25)
∫
Ω
‍
e
u
n
v
d
x
+
β
(
∫
Ω
‍
e
p
n
v
x
d
x
-
∑
j
=
1
N
‍
(
(
(
e
p
n
)
+
v
-
)
j
+
(
1
/
2
)
+
(
(
e
p
n
)
+
v
+
)
j
-
(
1
/
2
)
)
)
-
∑
i
=
1
n
-
1
‍
(
b
i
-
1
-
b
i
)
∫
Ω
‍
e
u
n
-
i
v
d
x
-
b
n
-
1
∫
Ω
‍
e
u
0
v
d
x
+
β
∫
Ω
‍
γ
n
(
x
)
v
d
x
+
∫
Ω
‍
e
p
n
w
d
x
+
∫
Ω
‍
e
u
n
w
x
d
x
-
∑
j
=
1
N
‍
(
(
(
e
u
n
)
-
w
-
)
j
+
(
1
/
2
)
+
(
(
e
u
n
)
-
w
+
)
j
-
(
1
/
2
)
)
=
0
.
Using (24), the error equation (25) can be written as
(26)
∫
Ω
‍
P
-
e
u
n
v
d
x
+
β
(
∫
Ω
‍
P
+
e
p
n
v
x
d
x
-
∑
j
=
1
N
‍
(
(
(
P
+
e
p
n
)
+
v
-
)
j
+
(
1
/
2
)
+
(
(
P
+
e
p
n
)
+
v
+
)
j
-
(
1
/
2
)
)
)
+
∫
Ω
‍
P
+
e
p
n
w
d
x
+
∫
Ω
‍
P
-
e
u
n
w
x
d
x
-
∑
j
=
1
N
‍
(
(
(
P
-
e
u
n
)
-
w
-
)
j
+
(
1
/
2
)
+
(
(
P
-
e
u
n
)
-
w
+
)
j
-
(
1
/
2
)
)
=
∑
i
=
1
n
-
1
‍
(
b
i
-
1
-
b
i
)
∫
Ω
‍
P
-
e
u
n
-
i
v
d
x
+
b
n
-
1
∫
Ω
‍
P
-
e
u
0
v
d
x
-
β
∫
Ω
‍
γ
n
(
x
)
v
d
x
+
∫
Ω
‍
(
P
-
u
(
x
,
t
n
)
-
u
(
x
,
t
n
)
)
v
d
x
+
β
(
∫
Ω
‍
(
P
+
p
(
x
,
t
n
)
-
p
(
x
,
t
n
)
)
v
x
d
x
-
∑
j
=
1
N
‍
(
(
(
P
+
p
(
x
,
t
n
)
-
p
(
x
,
t
n
)
)
+
v
-
)
j
+
(
1
/
2
)
+
(
(
P
+
p
(
x
,
t
n
)
-
p
(
x
,
t
n
)
)
+
v
+
)
j
-
(
1
/
2
)
)
(
∫
Ω
‍
(
P
+
p
(
x
,
t
n
)
-
p
(
x
,
t
n
)
)
v
x
d
x
)
-
∑
i
=
1
n
-
1
‍
(
b
i
-
1
-
b
i
)
∫
Ω
‍
(
P
-
u
(
x
,
t
n
-
i
)
-
u
(
x
,
t
n
-
i
)
)
v
d
x
-
b
n
-
1
∫
Ω
‍
(
P
-
u
(
x
,
t
0
)
-
u
(
x
,
t
0
)
)
v
d
x
+
∫
Ω
‍
(
P
+
p
(
x
,
t
n
)
-
p
(
x
,
t
n
)
)
w
d
x
+
∫
Ω
‍
(
P
-
u
(
x
,
t
n
)
-
u
(
x
,
t
n
)
)
w
x
d
x
-
∑
j
=
1
N
‍
(
(
(
P
-
u
(
x
,
t
n
)
-
u
(
x
,
t
n
)
)
-
w
-
)
j
+
(
1
/
2
)
+
(
(
P
-
u
(
x
,
t
n
)
-
u
(
x
,
t
n
)
)
-
w
+
)
j
-
(
1
/
2
)
)
.
Taking the test functions v
=
P
-
e
u
n
,
w
=
β
P
+
e
p
n in (26), using the properties (4) and (6), then the following equality holds:
(27)
∫
Ω
‍
(
P
-
e
u
n
)
2
d
x
+
β
∫
Ω
‍
(
P
+
e
p
n
)
2
d
x
=
∑
i
=
1
n
-
1
‍
(
b
i
-
1
-
b
i
)
∫
Ω
‍
P
-
e
u
n
-
i
P
-
e
u
n
d
x
+
b
n
-
1
∫
Ω
‍
P
-
e
u
0
P
-
e
u
n
d
x
-
β
∫
Ω
‍
γ
n
(
x
)
P
-
e
u
n
d
x
+
∫
Ω
‍
(
P
-
u
(
x
,
t
n
)
-
u
(
x
,
t
n
)
)
P
-
e
u
n
d
x
+
β
(
∫
Ω
‍
(
P
+
p
(
x
,
t
n
)
-
p
(
x
,
t
n
)
)
(
P
-
e
u
n
)
x
d
x
-
∑
j
=
1
N
‍
(
(
(
P
+
p
(
x
,
t
n
)
-
p
(
x
,
t
n
)
)
+
(
P
-
e
u
n
)
-
)
j
+
(
1
/
2
)
+
(
(
P
+
p
(
x
,
t
n
)
-
p
(
x
,
t
n
)
)
+
(
P
-
e
u
n
)
+
)
j
-
(
1
/
2
)
)
(
∫
Ω
‍
(
P
+
p
(
x
,
t
n
)
-
p
(
x
,
t
n
)
)
(
P
-
e
u
n
)
x
d
x
)
-
∑
i
=
1
n
-
1
‍
(
b
i
-
1
-
b
i
)
∫
Ω
‍
(
P
-
u
(
x
,
t
n
-
i
)
-
u
(
x
,
t
n
-
i
)
)
P
-
e
u
n
d
x
-
b
n
-
1
∫
Ω
‍
(
P
-
u
(
x
,
t
0
)
-
u
(
x
,
t
0
)
)
P
-
e
u
n
d
x
+
∫
Ω
‍
(
P
+
p
(
x
,
t
n
)
-
p
(
x
,
t
n
)
)
(
β
P
+
e
p
n
)
d
x
+
∫
Ω
‍
(
P
-
u
(
x
,
t
n
)
-
u
(
x
,
t
n
)
)
(
β
P
+
e
p
n
)
x
d
x
-
∑
j
=
1
N
‍
(
(
(
P
-
u
(
x
,
t
n
)
-
u
(
x
,
t
n
)
)
-
(
β
P
+
e
p
n
)
-
)
j
+
(
1
/
2
)
+
(
(
P
-
u
(
x
,
t
n
)
-
u
(
x
,
t
n
)
)
-
(
β
P
+
e
p
n
)
+
)
j
-
(
1
/
2
)
)
=
∑
i
=
1
n
-
1
‍
(
b
i
-
1
-
b
i
)
∫
Ω
‍
P
-
e
u
n
-
i
P
-
e
u
n
d
x
+
b
n
-
1
∫
Ω
‍
P
-
e
u
0
P
-
e
u
n
d
x
-
β
∫
Ω
‍
γ
n
(
x
)
P
-
e
u
n
d
x
+
∫
Ω
‍
(
P
-
u
(
x
,
t
n
)
-
u
(
x
,
t
n
)
)
P
-
e
u
n
d
x
-
∑
i
=
1
n
-
1
‍
(
b
i
-
1
-
b
i
)
∫
Ω
‍
(
P
-
u
(
x
,
t
n
-
i
)
-
u
(
x
,
t
n
-
i
)
)
P
-
e
u
n
d
x
-
b
n
-
1
∫
Ω
‍
(
P
-
u
(
x
,
t
0
)
-
u
(
x
,
t
0
)
)
P
-
e
u
n
d
x
+
∫
Ω
‍
(
P
+
p
(
x
,
t
n
)
-
p
(
x
,
t
n
)
)
(
β
P
+
e
p
n
)
d
x
;
that is,
(28)
∫
Ω
‍
(
P
-
e
u
n
)
2
d
x
+
β
∫
Ω
‍
(
P
+
e
p
n
)
2
d
x
≤
(
∑
i
=
1
n
-
1
‍
(
b
i
-
1
-
b
i
)
∥
P
-
e
u
n
-
i
∥
+
b
n
-
1
∥
P
-
e
u
0
∥
+
β
∥
γ
n
(
x
)
∥
+
∥
P
-
u
(
x
,
t
n
)
-
u
(
x
,
t
n
)
∥
+
∑
i
=
1
n
-
1
‍
(
b
i
-
1
-
b
i
)
∥
P
-
u
(
x
,
t
n
-
i
)
-
u
(
x
,
t
n
-
i
)
∥
+
b
n
-
1
∥
P
-
u
(
x
,
t
0
)
-
u
(
x
,
t
0
)
∥
(
∑
i
=
1
n
-
1
‍
(
b
i
-
1
-
b
i
)
∥
P
-
e
u
n
-
i
∥
)
∥
P
-
e
u
n
∥
+
β
∥
P
+
p
(
x
,
t
n
)
-
p
(
x
,
t
n
)
∥
∥
P
+
e
p
n
∥
.
Based on the fact that a
2
+
b
2
≤
(
a
+
b
)
2, we can obtain
(29)
∥
P
-
e
u
n
∥
≤
∑
i
=
1
n
-
1
‍
(
b
i
-
1
-
b
i
)
∥
P
-
e
u
n
-
i
∥
+
b
n
-
1
∥
P
-
e
u
0
∥
+
β
∥
γ
n
(
x
)
∥
+
∥
P
-
u
(
x
,
t
n
)
-
u
(
x
,
t
n
)
∥
+
∑
i
=
1
n
-
1
‍
(
b
i
-
1
-
b
i
)
∥
P
-
u
(
x
,
t
n
-
i
)
-
u
(
x
,
t
n
-
i
)
∥
+
b
n
-
1
∥
P
-
u
(
x
,
t
0
)
-
u
(
x
,
t
0
)
∥
+
β
∥
P
+
p
(
x
,
t
n
)
-
p
(
x
,
t
n
)
∥
.
For the sake of convenience, we denote
(30)
β
=
O
(
(
Δ
t
)
α
)
=
C
2
(
Δ
t
)
α
.
(1) We start with the following estimate:
(31)
∥
P
-
e
u
n
∥
≤
b
n
-
1
-
1
C
(
h
k
+
1
+
(
Δ
t
)
2
+
(
Δ
t
)
α
/
2
h
k
+
(
1
/
2
)
)
.
When n
=
1, (35) becomes
(32)
∥
P
-
e
u
1
∥
≤
∥
P
-
e
u
0
∥
+
β
∥
γ
1
(
x
)
∥
+
∥
P
-
u
(
x
,
t
1
)
-
u
(
x
,
t
1
)
∥
+
∥
P
-
u
(
x
,
t
0
)
-
u
(
x
,
t
0
)
∥
+
β
∥
P
+
p
(
x
,
t
1
)
-
p
(
x
,
t
1
)
∥
≤
C
1
C
2
(
Δ
t
)
2
+
3
C
0
h
k
+
1
+
C
2
C
0
(
Δ
t
)
α
/
2
h
k
+
1
.
Denoting C
=
max
⁡
{
3
C
0
,
C
1
C
2
,
C
2
C
0
}, then we can obtain
(33)
∥
P
-
e
u
1
∥
≤
b
0
-
1
C
(
(
Δ
t
)
2
+
h
k
+
1
+
(
Δ
t
)
α
/
2
h
k
+
1
)
.
Next, we suppose the following inequality holds:
(34)
∥
P
-
e
u
m
∥
≤
b
m
-
1
-
1
C
(
h
k
+
1
+
(
Δ
t
)
2
+
(
Δ
t
)
α
/
2
h
k
+
1
)
,
m
=
1,2
,
…
,
K
.
Let n
=
K
+
1 in the inequality (35); we deduce
(35)
∥
P
-
e
u
K
+
1
∥
≤
∑
i
=
1
K
‍
(
b
i
-
1
-
b
i
)
∥
P
-
e
u
K
+
1
-
i
∥
+
b
K
∥
P
-
e
u
0
∥
+
β
∥
γ
K
+
1
(
x
)
∥
+
∥
P
-
u
(
x
,
t
K
+
1
)
-
u
(
x
,
t
K
+
1
)
∥
+
∑
i
=
1
K
‍
(
b
i
-
1
-
b
i
)
∥
P
-
u
(
x
,
t
K
+
1
-
i
)
-
u
(
x
,
t
K
+
1
-
i
)
∥
+
b
K
∥
P
-
u
(
x
,
t
0
)
-
u
(
x
,
t
0
)
∥
+
β
∥
P
+
p
(
x
,
t
K
+
1
)
-
p
(
x
,
t
K
+
1
)
∥
≤
∑
i
=
1
K
‍
(
b
i
-
1
-
b
i
)
∥
P
-
e
u
K
+
1
-
i
∥
+
C
1
C
2
(
Δ
t
)
2
+
3
C
0
h
k
+
1
+
C
2
C
0
(
Δ
t
)
α
/
2
h
k
+
1
≤
∑
i
=
1
K
‍
(
b
i
-
1
-
b
i
)
b
K
-
i
-
1
C
(
h
k
+
1
+
(
Δ
t
)
2
+
(
Δ
t
)
α
/
2
h
k
+
1
)
+
C
(
h
k
+
1
+
(
Δ
t
)
2
+
(
Δ
t
)
α
/
2
h
k
+
1
)
.
Notice the fact that(36)
b
i
-
1
-
1
<
b
i
-
1
;
we know
(37)
∥
P
-
e
u
K
+
1
∥
≤
∑
i
=
1
K
‍
(
b
i
-
1
-
b
i
+
b
K
)
b
K
-
1
C
(
h
k
+
1
+
(
Δ
t
)
2
+
(
Δ
t
)
α
/
2
h
k
+
1
)
;
that is,
(38)
∥
P
-
e
u
K
+
1
∥
≤
b
K
-
1
C
(
h
k
+
1
+
(
Δ
t
)
2
+
(
Δ
t
)
α
/
2
h
k
+
1
)
.
Inequality (31) follows.
By some calculations and analyses, we know thatn
-
α
b
n
-
1
-
1 increasingly tends to 1
/
(
1
-
α
). For more details of the proof, we refer to [14]. So we can obtain
(39)
∥
P
-
e
u
n
∥
≤
b
n
-
1
-
1
C
(
h
k
+
1
+
(
Δ
t
)
2
+
(
Δ
t
)
α
/
2
h
k
+
1
)
≤
n
α
n
-
α
b
n
-
1
-
1
C
(
h
k
+
1
+
(
Δ
t
)
2
+
(
Δ
t
)
α
/
2
h
k
+
1
)
≤
C
T
α
1
-
α
(
(
Δ
t
)
-
α
h
k
+
1
+
(
Δ
t
)
2
-
α
+
(
Δ
t
)
-
α
/
2
h
k
+
1
)
.
(2) The above estimate has no meaning whenα
→
1 due to 1
/
(
1
-
α
)
→
∞. So we must reconsider it for the case α
→
1.
We suppose the following estimate holds:(40)
∥
P
-
e
u
n
∥
≤
n
C
(
h
k
+
1
+
(
Δ
t
)
2
+
(
Δ
t
)
α
/
2
h
k
+
1
)
.
By the similar techniques used in (1) and that in [14], we can obtain (40) easily. Here, we omitted the proof to save space. Then, we know that when α
→
1,
(41)
∥
P
-
e
u
n
∥
≤
T
C
(
(
Δ
t
)
-
1
h
k
+
1
+
Δ
t
+
(
Δ
t
)
-
1
/
2
h
k
+
1
)
.
Thus, Theorem2 follows by the triangle inequality and the interpolating property (7).


## 4. Numerical Examples

In this section, we offer some numerical examples to illustrate the accuracy and capability of the method. For this purpose, we calculate the numerical results of the exact solutions (for the cases where exact solutions are available). We mainly focus on the spatial accuracy, so a small time step is used such that errors stemming from the temporal approximation are negligible. With the aid of successive mesh refinements, we have verified that the results shown are numerically convergent.Example 1.
We consider time-fractional equation (1) in Ω
=
[
0,1
]; the corresponding forcing term f
(
x
,
t
) is of the form
(42)
f
(
x
,
t
)
=
2
t
2
-
α
Γ
(
3
-
α
)
sin
(
2
π
x
)
+
4
π
t
2
sin
(
2
π
x
)
;
Then, the exact solution is u
(
x
,
t
)
=
t
2
sin
(
2
π
x
). The space and time step is h
=
1
/
N
,
Δ
t
=
1
/
M, respectively. We check the spatial accuracy by fixing the time step sufficiently small to avoid contamination of the temporal error. From Tables 1, 2, 3, 4, 5, 6, 7, 8, and 9, we can see that the errors in L
2-norm and L
∞-norm attain optimal order of accuracy for piecewise P
k polynomials for α
=
0.1,0.2
,
…
,
0.9. In Table 10, we show the errors in L
1-norm and L
2-norm attains 2
-
α order of accuracy for two values of α
:
0.5 and 0.7.Table 1
Spatial accuracy test for the fractional order equation with the forcing term (42). α
=
0.1, M
=
1
0
4, T
=
1.
N
L
2-Error
Order
L
∞-Error
Order
P
0
5
0.266079725401477
—
0.625909610947696
—
10
0.129418994391760
1.04
0.313875271857693
1.00
15
8.584650747857661
E
-
002
1.01
0.209356339112292
1.00
20
6.427202957081933
E
-
002
1.01
0.157044742089558
1.00
P
1
5
6.743149391082712
E
-
002
—
0.250201723135112
—
10
1.695823938890830
E
-
002
1.99
6.470264790292768
E
-
002
1.95
15
7.544622677866573
E
-
003
2.00
2.865534240293233
E
-
002
2.01
20
4.245325459728799
E
-
003
2.00
1.631373162490857
E
-
002
1.96
P
2
5
6.687590470572136
E
-
003
—
3.176309264202759
E
-
002
—
10
8.508409837025869
E
-
004
2.97
3.972045444261703
E
-
003
3.00
15
2.529255442619073
E
-
004
2.99
1.219537695428569
E
-
003
2.91
20
1.068251766123077
E
-
004
3.00
5.116601796458933
E
-
004
3.02Table 2
Spatial accuracy test for the fractional order equation with the forcing term (42). α
=
0.2, M
=
1
0
4, T
=
1.
N
L
2-Error
Order
L
∞-Error
Order
P
0
5
0.265992689941441
—
0.625698918739395
—
10
0.129409013964343
1.04
0.313850364903468
1.00
15
8.584359829944356
E
-
002
1.01
0.209349038923774
1.00
20
6.427081020663458
E
-
002
1.01
0.157041676426717
1.00
P
1
5
6.742575102958956
E
-
002
—
0.250172479773477
—
10
1.695793298787001
E
-
002
1.99
6.470108658023166
E
-
002
1.95
15
7.544564047272043
E
-
003
2.00
2.865508053806298
E
-
002
2.01
20
4.245307149025943
E
-
003
2.00
1.631368024213253
E
-
002
1.96
P
2
5
6.687158434739263
E
-
003
—
3.176105007478780
E
-
002
—
10
8.508277902535355
E
-
004
2.97
3.971971730332332
E
-
003
3.00
15
2.529238517808841
E
-
004
2.99
1.219537991111880
E
-
003
2.91
20
1.068248546234930
E
-
004
3.00
5.116578768526442
E
-
004
3.02Table 3
Spatial accuracy test for the fractional order equation with the forcing term (42). α
=
0.3, M
=
1
0
4, T
=
1.
N
L
2-Error
Order
L
∞-Error
Order
P
0
5
0.265902352109771
—
0.625479684568295
—
10
0.129398459085195
1.04
0.313823958040159
1.00
15
8.584042638776486
E
-
002
1.01
0.209341059113568
1.00
20
6.426942470350358
E
-
002
1.01
0.157038183985132
1.00
P
1
5
6.741974383634558
E
-
002
—
0.250139776852748
—
10
1.695759921523912
E
-
002
1.99
6.469739429987975
E
-
002
1.95
15
7.544495988113592
E
-
003
2.00
2.865280596382414
E
-
002
2.01
20
4.245284134117610
E
-
003
2.00
1.631154850118743
E
-
002
1.96
P
2
5
6.686707388830991
E
-
003
—
3.175767063389789
E
-
002
—
10
8.508142602066554
E
-
004
2.97
3.971887377547693
E
-
003
3.00
15
2.529242023200877
E
-
004
2.99
1.219110576456323
E
-
003
2.91
20
1.068298199014039
E
-
004
3.00
5.116544492389075
E
-
004
3.02Table 4
Spatial accuracy test for the fractional order equation with the forcing term (42). α
=
0.4, M
=
1
0
4, T
=
1.
N
L
2-Error
Order
L
∞-Error
Order
P
0
5
0.265810951781775
—
0.625257299208153
—
10
0.129388142089229
1.04
0.313798079922729
1.00
15
8.583749784514737
E
-
002
1.01
0.209333672586633
1.00
20
6.426824345720851
E
-
002
1.01
0.157035198723482
1.00
P
1
5
6.741372975727104
E
-
002
—
0.250110840146267
—
10
1.695728922719653
E
-
002
1.99
6.469745500225199
E
-
002
1.95
15
7.544440107132382
E
-
003
2.00
2.865418145532883
E
-
002
2.01
20
4.245268103040704
E
-
003
2.00
1.631320854177420
E
-
002
1.96
P
2
5
6.686254118771560
E
-
003
—
3.175655562951962
E
-
002
—
10
8.508000331773612
E
-
004
2.97
3.971816136897089
E
-
003
3.00
15
2.529200607167737
E
-
004
2.99
1.219463687311600
E
-
003
2.91
20
1.068236428236778
E
-
004
3.00
5.116528775108040
E
-
004
3.02Table 5
Spatial accuracy test for the fractional order equation with the forcing term (42). α
=
0.5, M
=
1
0
4, T
=
1.
N
L
2-Error
Order
L
∞-Error
Order
P
0
5
0.265717546415908
—
0.625029434969908
—
10
0.129377426180298
1.04
0.313771131147625
1.00
15
8.583436861020562
E
-
002
1.01
0.209325759569002
1.00
20
6.426692835664491
E
-
002
1.01
0.157031866794147
1.00
P
1
5
6.740756825656390
E
-
002
—
0.250079212525491
—
10
1.695695972605060
E
-
002
1.99
6.469564272790318
E
-
002
1.95
15
7.544376786237617
E
-
003
2.00
2.865377151494430
E
-
002
2.01
20
4.245248206073015
E
-
003
2.00
1.631302015908276
E
-
002
1.96
P
2
5
6.685790666908454
E
-
003
—
3.175428377723111
E
-
002
—
10
8.507858388639393
E
-
004
2.97
3.971736600682197
E
-
003
3.00
15
2.529181785104898
E
-
004
2.99
1.219436646278105
E
-
003
2.91
20
1.068231569593730
E
-
004
3.00
5.116503423953015
E
-
004
3.02Table 6
Spatial accuracy test for the fractional order equation with the forcing term (42). α
=
0.6, M
=
1
0
4, T
=
1.
N
L
2-Error
Order
L
∞-Error
Order
P
0
5
0.265623584016600
—
0.624799593887340
—
10
0.129366663521458
1.04
0.313743992355115
1.00
15
8.583123035523127
E
-
002
1.01
0.209317802516665
1.00
20
6.426561199139737
E
-
002
1.01
0.157028522738007
1.00
P
1
5
6.740140755904378
E
-
002
—
0.250047569176264
—
10
1.695663105052171
E
-
002
1.99
6.469391604272101
E
-
002
1.95
15
7.544313817414976
E
-
003
2.00
2.865344855707630
E
-
002
2.01
20
4.245228502731458
E
-
003
2.00
1.631292276055152
E
-
002
1.96
P
2
5
6.685327402126262
E
-
003
—
3.175206692006394
E
-
002
—
10
8.507716798808087
E
-
004
2.97
3.971657446348387
E
-
003
3.00
15
2.529163241406188
E
-
004
2.99
1.219428402258027
E
-
003
2.91
20
1.068227212010240
E
-
004
3.00
5.116478542937988
E
-
004
3.02Table 7
Spatial accuracy test for the fractional order equation with the forcing term (42). α
=
0.7, M
=
1
0
4, T
=
1.
N
L
2-Error
Order
L
∞-Error
Order
P
0
5
0.265530364772941
—
0.624570952554277
—
10
0.129356176022042
1.04
0.313717476972879
1.00
15
8.582825921887373
E
-
002
1.01
0.209310249447525
1.00
20
6.426441648813166
E
-
002
1.01
0.157025477918926
1.00
P
1
5
6.739538689689466
E
-
002
—
0.250018407776982
—
10
1.695632222252693
E
-
002
1.99
6.469410377138518
E
-
002
1.95
15
7.544258629091435
E
-
003
2.00
2.865496091236286
E
-
002
2.01
20
4.245213110963431
E
-
003
2.00
1.631472721918670
E
-
002
1.96
P
2
5
6.684874322990170
E
-
003
—
3.175103718965400
E
-
002
—
10
8.507593551652111
E
-
004
2.97
3.971586800980554
E
-
003
3.00
15
2.529184588627819
E
-
004
2.99
1.219811466072696
E
-
003
2.91
20
1.068313817637621
E
-
004
3.00
5.116463574011057
E
-
004
3.02Table 8
Spatial accuracy test for the fractional order equation with the forcing term (42). α
=
0.8, M
=
1
0
4, T
=
1.
N
L
2-Error
Order
L
∞-Error
Order
P
0
5
0.265437954148352
—
0.624343681308530
—
10
0.129345622072736
1.04
0.313690722834305
1.00
15
8.582518377117876
E
-
002
1.01
0.209302410813452
1.00
20
6.426312693843939
E
-
002
1.01
0.157022185190115
1.00
P
1
5
6.738947359969441
E
-
002
—
0.249987673730495
—
10
1.695600756044601
E
-
002
1.99
6.469244216241810
E
-
002
1.95
15
7.544198389450005
E
-
003
2.00
2.865465999545802
E
-
002
2.01
20
4.245194275362235
E
-
003
2.00
1.631464536623917
E
-
002
1.96
P
2
5
6.684430317357719
E
-
003
—
3.174891822494708
E
-
002
—
10
8.507458300947125
E
-
004
2.97
3.971511081252594
E
-
003
3.00
15
2.529167702581917
E
-
004
2.99
1.219806186412653
E
-
003
2.91
20
1.068311670673490
E
-
004
3.00
5.116439829187727
E
-
004
3.02Table 9
Spatial accuracy test for the fractional order equation with the forcing term (42). α
=
0.9, M
=
1
0
4, T
=
1.
N
L
2-Error
Order
L
∞-Error
Order
P
0
5
0.265347755163247
—
0.624121255088276
—
10
0.129335227490933
1.04
0.313664302744807
1.00
15
8.582209936957509
E
-
002
1.01
0.209294528357869
1.00
20
6.426180057018088
E
-
002
1.01
0.157018789299366
1.00
P
1
5
6.738379605827277
E
-
002
—
0.249956615433034
—
10
1.695569809610291
E
-
002
1.99
6.468959039144961
E
-
002
1.95
15
7.544136529032524
E
-
003
2.00
2.865315994055062
E
-
002
2.01
20
4.245173657220278
E
-
003
2.00
1.631331043368178
E
-
002
1.96
P
2
5
6.684004911093908
E
-
003
—
3.174613402126625
E
-
002
—
10
8.507316247783669
E
-
004
2.97
3.971434165464522
E
-
003
3.00
15
2.529116662199529
E
-
004
2.99
1.219542692523057
E
-
003
2.91
20
1.068228730392493
E
-
004
3.00
5.116410932367235
E
-
004
3.02Table 10
Temporal accuracy test for the problem (1) with the forcing term (42) when N
=
100, T
=
1, and k
=
2.
Δ
t
L
2-Error
Order
L
1-Error
Order
α
=
0.5
0.04
6.361850145394739
E
-
005
—
5.728850142441670
E
-
005
—
0.02
2.281473637510036
E
-
005
1.48
2.054898402168884
E
-
005
1.48
0.01
8.204587459292379
E
-
006
1.48
7.392564492689201
E
-
006
1.47
0.005
3.053359784205133
E
-
006
1.43
2.746671417581169
E
-
006
1.43
α
=
0.7
0.04
1.767309239658665
E
-
004
—
1.591263361456412
E
-
004
—
0.02
7.281830517994945
E
-
005
1.28
6.557148539250536
E
-
005
1.28
0.01
3.042807494504462
E
-
005
1.26
2.740365369523097
E
-
005
1.26
0.005
1.317960860254862
E
-
005
1.21
1.187237948375165
E
-
005
1.21


## 5. Conclusion

In this paper, an implicit fully discrete local discontinuous Galerkin (LDG) finite element method is presented for solving a class of time-fractional diffusion equation. Numerical examples show that the combination of the backward differentiation in time and local discontinuous Galerkin (LDG) finite element method in space leads to an approximation of order(
(
Δ
t
)
2
-
α
+
h
k
+
1
) for smooth enough solution. The scheme can be extended to solve the two or higher dimensional case easily, and the theoretical results are also valid. The results show that the LDG method is a powerful and efficient technique in solving this class of problems with fractional derivatives.



---

*Source: 898217-2014-08-19.xml*