Abstract
A standard finite difference method on a uniform mesh is used to solve a timefractional convectiondiffusion initialboundary value problem. Such problems typically exhibit a mild singularity at the initial time
t
=
0
. It is proved that the rate of convergence of the maximum nodal error on any subdomain that is bounded away from
t
=
0
is higher than the rate obtained when the maximum nodal error is measured over the entire spacetime domain. Numerical results are provided to illustrate the theoretical error bounds.
1 Introduction
In this paper, we examine the convergence rate of numerical approximations to a timefractional convectiondiffusion problem using a standard finite difference method on a uniform mesh. Initialboundary value problems of this type, where the time derivative is fractional, have solutions that are mildly singular at the initial time
t
=
0
; that is, their temporal derivatives are unbounded on the closed spacetime domain, but are bounded on any subdomain that is bounded away from
t
=
0
. It is shown in [9, 10] that the case of solutions with bounded temporal derivatives on the closed spacetime domain is very special and that the weakly singular solutions examined here are much more typical of how solutions to this class of problems behave. As one would expect, the rate of convergence of the computed numerical approximations is affected adversely by the presence of large temporal derivatives at
t
=
0
.
This paper is a companion paper to [10], where it was shown that the convergence rate of the same finite difference scheme on a uniform mesh was
O
(
N

α
)
, where
α
∈
(
0
,
1
)
is the order of the fractional derivative and the mesh spacing in time is
O
(
N

1
)
. Results related to the main result in [10] are available in [7, 8], using a finite element framework. In contrast to [10], we shall prove here for the same scheme on a uniform mesh that the convergence rate of the numerical solution is
O
(
N

1
)
on any subdomain that is bounded away from
t
=
0
.
Our analysis is carried out in the discrete
L
∞
norm; an analogous convergence result in the
L
2
norm was derived in [4]. Using an alternative formulation of the continuous problem, the phenomenon of higherorder convergence at some fixed distance away from the initial singularity is examined in [6] for a homogeneous version of (2.1) in the case of nonsmooth initial data.
Notation.
In this paper C denotes a generic constant that depends on the data of the boundary value problem (2.1) but is independent of T and of any mesh used to solve (2.1) numerically. Note that C can take different values in different places. For all
x
∈
ℝ
the ceiling function
⌈
x
⌉
is the smallest integer greater than or equal to x. For any continuous function
z
:
Q
→
ℝ
with
Q
⊂
ℝ
2
and any mesh function
z
m
n
with
n
=
0
,
1
,
…
,
N
and
m
=
0
,
1
,
…
,
M
, we set
∥
z
∥
:=
max
(
x
,
t
)
∈
Q
¯

z
(
x
,
t
)

and
∥
z
n
∥
:=
max
0
≤
m
≤
M

z
m
n

.
2 The Continuous Problem
Consider the initialboundary value problem
(2.1)
(2.1a)
L
u
:=
D
t
α
u

p
(
x
)
∂
2
u
∂
x
2
+
q
(
x
)
∂
u
∂
x
+
r
(
x
)
u
=
f
(
x
,
t
)
for
(
x
,
t
)
∈
Q
:=
(
0
,
l
)
×
(
0
,
T
]
, with initial and boundary conditions
(2.1b)
u
(
0
,
t
)
=
u
(
l
,
t
)
=
0
for
t
∈
(
0
,
T
]
,
(2.1c)
u
(
x
,
0
)
=
ϕ
(
x
)
for
x
∈
[
0
,
l
]
.
Here
0
<
α
<
1
,
p
(
x
)
≥
p
0
>
0
on
[
0
,
l
]
, the functions
p
(
x
)
,
q
(
x
)
and
r
(
x
)
are smooth on
[
0
,
l
]
and are assumed to satisfy the constraint
(2.1d)
r
(
x
)

q
′
(
x
)
+
p
′′
(
x
)
2
≥
0
for
(
x
,
t
)
∈
Q
.
The initial condition ϕ is also smooth on
[
0
,
l
]
and the function f is smooth on
Q
¯
. Furthermore, in (2.1a)
D
t
α
denotes the Caputo fractional derivative which is defined [1] by
D
t
α
g
(
x
,
t
)
:=
[
J
1

α
(
∂
g
∂
t
)
]
(
x
,
t
)
for
0
≤
x
≤
l
,
0
<
t
≤
T
,
where
(
J
1

α
g
)
(
x
,
t
)
:=
[
1
Γ
(
1

α
)
∫
s
=
0
t
(
t

s
)

α
g
(
x
,
s
)
𝑑
s
]
for
(
x
,
t
)
∈
Q
¯
is the Riemann–Liouville fractional integral operator of order
1

α
.
There is no loss of generality in assuming homogeneous boundary conditions in (2.1b), because inhomogeneous boundary conditions are easily made homogeneous by a simple change of variable.
Under the transformation
y
(
x
,
t
)
:=
u
(
x
,
t
)
p
(
0
)
p
(
x
)
e

1
2
∫
0
x
q
(
s
)
p
(
s
)
𝑑
s
,
problem (2.1) becomes
(2.2)
(2.2a)
D
t
α
y

p
(
x
)
∂
2
y
∂
x
2
+
r
1
(
x
)
y
=
f
1
(
x
,
t
)
for
(
x
,
t
)
∈
Q
:=
(
0
,
l
)
×
(
0
,
T
]
,
(2.2b)
y
(
0
,
t
)
=
y
(
l
,
t
)
=
0
for
t
∈
(
0
,
T
]
,
(2.2c)
y
(
x
,
0
)
=
ϕ
1
(
x
)
:=
ϕ
(
x
)
p
(
0
)
p
(
x
)
e

1
2
∫
0
x
q
(
s
)
p
(
s
)
𝑑
s
for
x
∈
[
0
,
l
]
,
(2.2d)
where
r
1
(
x
)
:=
r
(
x
)

(
q
′
(
x
)
+
p
′′
(
x
)
)
2
+
(
q
(
x
)
+
p
′
(
x
)
)
2
4
p
,
(2.2e)
f
1
(
x
,
t
)
:=
f
(
x
,
t
)
p
(
0
)
p
(
x
)
e

1
2
∫
0
x
q
(
s
)
p
(
s
)
𝑑
s
.
Note that no firstorder derivative in space appears in (2.2a), and (2.1d) implies that
r
1
≥
0
. Consequently, (2.2) belongs to the class of problems analysed in [10]. In [10] it was assumed that
p
(
x
)
was a positive constant, but the analysis of [10] can be extended to the case of a smooth variable positive coefficient
p
(
x
)
. Thus after placing suitable regularity and compatibility conditions on the data of (2.2), one can invoke [10, Theorem 2.1] to conclude that (2.2) has a unique solution y whose derivatives satisfy certain bounds. Transforming back to the original problem (2.1), under certain conditions on its data one obtains the following bounds on the derivatives of u:
(2.3)
(2.3a)

∂
k
u
∂
x
k
(
x
,
t
)

≤
C
for
k
=
0
,
1
,
2
,
3
,
4
,
(2.3b)

∂
ℓ
u
∂
t
ℓ
(
x
,
t
)

≤
C
(
1
+
t
α

ℓ
)
for
ℓ
=
1
,
2
,
for all
(
x
,
t
)
∈
[
0
,
l
]
×
(
0
,
T
]
.
In [10, Theorem 2.1] the estimates in (2.3) are proved assuming that
ϕ
1
∈
D
(
ℒ
5
/
2
)
,
(
f
1
)
(
⋅
,
t
)
∈
D
(
ℒ
5
/
2
)
,
(
f
1
)
t
(
⋅
,
t
)
and
(
f
1
)
t
t
(
⋅
,
t
)
are in
D
(
ℒ
1
/
2
)
for each
t
∈
(
0
,
T
]
and
∥
(
f
1
)
(
⋅
,
t
)
∥
ℒ
5
/
2
+
∥
(
f
1
)
t
(
⋅
,
t
)
∥
ℒ
1
/
2
+
t
ρ
∥
(
f
1
)
t
t
(
⋅
,
t
)
∥
ℒ
1
/
2
≤
C
1
for all
t
∈
(
0
,
T
]
and some constant
ρ
<
1
, where
C
1
is a constant independent of t and
∥
⋅
∥
ℒ
γ
is the norm associated with the vector space
D
(
ℒ
γ
)
. This space is defined by
D
(
ℒ
γ
)
:=
{
g
∈
L
2
(
0
,
l
)
:
∑
i
=
1
∞
λ
i
2
γ

(
g
,
ψ
i
)

2
<
∞
}
,
γ
≥
0
,
where
(
⋅
,
⋅
)
is the inner product in the Hilbert space
L
2
(
0
,
l
)
and
{
(
λ
i
,
ψ
i
)
:
i
=
1
,
2
,
…
}
are the eigenvalues and normalised eigenfunctions of the Sturm–Liouville twopoint boundary value problem
ℒ
ψ
i
:=

p
ψ
i
′′
+
c
ψ
i
=
λ
i
ψ
i
on
(
0
,
l
)
,
ψ
i
(
0
)
=
ψ
i
(
l
)
=
0
.
3 The Discrete Problem
The solution of problem (2.1) is approximated by the solution of a finite difference scheme on a mesh
{
(
x
m
,
t
n
)
:
m
=
0
,
1
,
…
,
M
,
n
=
0
,
1
,
…
,
N
}
, that is uniform in both space and time. Let M and N be positive integers. Set
h
=
l
M
and
x
m
:=
m
h
for
m
=
0
,
1
,
…
,
M
. Set
t
n
=
n
τ
=
n
T
N
for
n
=
0
,
1
,
…
,
N
. The nodal approximation to the solution u computed at the mesh point
(
x
m
,
t
n
)
is denoted by
u
m
n
.
The first and secondorder spatial derivatives are discretised using standard approximations:
∂
u
∂
x
(
x
m
,
t
n
)
≈
D
x
0
u
m
n
:=
u
m
+
1
n

u
m

1
n
2
h
,
∂
2
u
∂
x
2
(
x
m
,
t
n
)
≈
δ
x
2
u
n
m
:=
u
m
+
1
n

2
u
m
n
+
u
m

1
n
h
2
.
The Caputo fractional derivative
D
t
α
u
, which can be written as
D
t
α
u
(
x
m
,
t
n
)
=
1
Γ
(
1

α
)
∑
k
=
0
n

1
∫
s
=
t
k
t
k
+
1
(
t
n

s
)

α
∂
u
(
x
m
,
s
)
∂
s
𝑑
s
,
is approximated by the classical L1 approximation
(3.1)
D
N
α
u
m
n
:=
1
Γ
(
1

α
)
∑
k
=
0
n

1
u
m
k
+
1

u
m
k
τ
∫
s
=
t
k
t
k
+
1
(
t
n

s
)

α
𝑑
s
(3.1a)
=
1
Γ
(
2

α
)
∑
k
=
0
n

1
(
u
m
k
+
1

u
m
k
)
d
n

k
,
where
(3.1b)
d
k
:=
k
1

α

(
k

1
)
1

α
,
k
≥
1
.
Thus, we approximate (2.1) by the discrete problem
(3.2)
(3.2a)
L
M
,
N
u
m
n
:=
D
N
α
u
m
n

p
(
x
m
)
δ
x
2
u
m
n
+
q
(
x
m
)
D
x
0
u
m
n
+
r
(
x
m
)
u
m
n
=
f
(
x
m
,
t
n
)
for
1
≤
m
≤
M

1
,
1
≤
n
≤
N
,
(3.2b)
u
0
n
=
u
M
n
=
0
for
0
<
n
≤
N
,
(3.2c)
u
m
0
=
ϕ
(
x
m
)
for
0
≤
m
≤
M
.
This discretisation of (2.1) is standard.
To ensure the stability of the discrete operator
L
M
,
N
by imposing the correct sign pattern in the associated matrix, we make the nonrestrictive assumption that N satisfies
l
∥
q
∥
2
p
0
<
N
.
After some minor modifications in the proof of [10, Theorem 5.2] to handle the term
q
(
x
m
)
D
x
0
u
m
n
, it follows that the solution
u
m
n
of scheme (3.2) satisfies the error bound
(3.3)
max
(
x
m
,
t
n
)
∈
Q
¯

u
(
x
m
,
t
n
)

u
m
n

≤
C
(
h
2
+
N

α
)
for some constant C. In particular, the method has the low order of convergence
O
(
N

α
)
in time when α is small. In the present paper we shall consider the rate of convergence in a subdomain
[
0
,
l
]
×
[
κ
,
T
]
, where κ is a fixed positive value.
4 Error Analysis
The structure of our error analysis is the standard finite difference technique of estimating the truncation error at each mesh point, then invoking a stability argument to derive an error bound for the computed solution
u
m
n
. In this analysis the truncation error bound (4.2) indicates that the truncation error decreases as one moves further away from the initial time
t
=
0
. The stability bound (4.5) shows that the error at any discrete time level depends on a weighted sum of the truncation errors at all the previous time levels.
The estimate of the truncation error in space is standard: using (2.3a), one gets
(4.1)
∂
u
∂
x
(
x
m
,
t
n
)
=
D
x
0
u
(
x
m
,
t
n
)
+
O
(
h
2
)
,
∂
2
u
∂
x
2
(
x
m
,
t
n
)
=
δ
x
2
u
(
x
m
,
t
n
)
+
O
(
h
2
)
.
The truncation error in time is more tricky to estimate and this is done in the next lemma.
Lemma 1.
Assume that u satisfies (2.3). Then there exists a positive constant C such that for each mesh point
(
x
m
,
t
n
)
∈
Q
one has
(4.2)

(
D
N
α

D
t
α
)
u
(
x
m
,
t
n
)

≤
C
n

min
{
2

α
,
α
+
1
}
.
Proof.
We modify the argument of [10, Lemma 5.1]. By (3.1a) and the definition of
D
t
α
, for each mesh point
(
x
m
,
t
n
)
∈
Q
one has
(
D
N
α

D
t
α
)
u
(
x
m
,
t
n
)
=
∑
k
=
0
n

1
T
n
k
,
where for
n
=
1
,
2
,
…
,
N
and
k
=
0
,
1
,
…
,
n

1
we define the truncation error in the kth time cell
[
t
k
,
t
k
+
1
]
to be
(4.3)
T
n
k
:=
1
Γ
(
1

α
)
∫
s
=
t
k
t
k
+
1
(
t
n

s
)

α
[
u
(
x
m
,
t
k
+
1
)

u
(
x
m
,
t
k
)
τ

∂
u
∂
s
(
x
m
,
s
)
]
𝑑
s
.
The following four bounds are established in [10, equations (5.9), (5.10), (5.11) and (5.14)]:
(4.4)
(4.4a)
∑
k
=
1
⌈
n
2
⌉

1

T
n
k

≤
C
n

(
α
+
1
)
for
1
≤
k
<
n

1
,
(4.4b)
∑
k
=
⌈
n
2
⌉
n

2

T
n
k

≤
C
n

(
2

α
)
for
1
≤
k
<
n

1
,
(4.4c)

T
10

≤
C
,
(4.4d)

T
n
,
n

1

≤
C
n

(
2

α
)
.
It remains to bound

T
n
0

for
n
>
1
. We sharpen the bound [10, (5.12)] of this term. An integration by parts in (4.3) yields
T
n
0
=

α
Γ
(
1

α
)
∫
s
=
0
t
1
(
t
n

s
)

α

1
(
ϕ

ψ
)
(
x
m
,
s
)
𝑑
s
,
where
ϕ
(
x
m
,
s
)
:=
s
[
u
(
x
m
,
t
1
)

u
(
x
m
,
0
)
τ
]
and
ψ
(
x
m
,
s
)
:=
u
(
x
m
,
s
)

u
(
x
m
,
0
)
.
For
0
≤
s
≤
τ
, it is clear that

ϕ
(
x
m
,
s
)

≤

u
(
x
m
,
τ
)

u
(
x
m
,
0
)

and

ψ
(
x
m
,
s
)

≤
∫
0
τ

u
t
(
x
m
,
t
)

𝑑
t
. Thus, we see that

ϕ
(
x
m
,
s
)

+

ψ
(
x
m
,
s
)

≤
2
∫
0
τ

u
t
(
x
m
,
t
)

𝑑
t
≤
C
∫
0
τ
(
1
+
t
α

1
)
𝑑
t
≤
C
τ
α
.
Hence

T
n
0

≤
C
τ
α
∫
s
=
0
t
1
(
t
n

s
)

α

1
𝑑
s
≤
C
τ
α
[
(
t
n

t
1
)

α

t
n

α
]
=
C
[
(
n

1
)

α

n

α
]
≤
C
n

(
α
+
1
)
,
by the Mean Value Theorem. Combine this bound with (4.4) to complete the proof. ∎
Observe that
min
{
2

α
,
α
+
1
}
>
1
in (4.2) for all values of
α
∈
(
0
,
1
)
; thus this bound is sharper than the truncation error bound of
O
(
n

α
)
proved in [10, Lemma 5.1]. This improvement is critical in establishing our main result later.
Next, we derive some new information about the stability constants that appear in [10, Section 4]. It follows from [10, Lemma 4.2] that the computed solution
u
m
n
of (3.2) satisfies
(4.5)
∥
u
(
x
m
,
t
n
)

u
m
n
∥
≤
τ
α
Γ
(
2

α
)
∑
j
=
1
n
σ
n

j
∥
(
L
M
,
N
(
u
(
x
m
,
t
j
)
)
j

u
j
∥
for
n
=
1
,
2
,
…
,
N
,
where the positive weights
σ
i
are defined for
i
=
0
,
1
,
2
,
…
,
n

1
by the recurrence relation
(4.6)
σ
0
:=
1
,
σ
i
:=
∑
k
=
1
i
(
d
k

d
k
+
1
)
σ
i

k
for
i
=
1
,
2
,
…
.
Note that when the mesh is uniform, the weights
θ
n
,
j
in [10, Lemma 4.2] are the same as the weights
σ
n

j
defined in (4.6).
Lemma 2.
The coefficients
σ
i
satisfy
σ
i
<
(
i
+
1
)
α

1
for
i
=
1
,
2
,
…
.
Proof.
First,
σ
1
=
(
d
1

d
2
)
σ
0
=
2

2
1

α
<
2

1
+
α
as
0.5
w
+
2
w

1

2
>
0
for all
w
∈
(
1
,
2
)
, so the lemma is true when
i
=
1
. The proof is completed by induction. Assume that
σ
j
<
(
j
+
1
)
α

1
for
j
=
1
,
2
,
…
,
i

1
. We want to prove that
σ
i
<
(
i
+
1
)
α

1
. It is easy to check that
d
k

d
k
+
1
>
0
for all k. Using this inequality and the inductive hypothesis, we require the inequality
∑
k
=
1
i
(
d
k

d
k
+
1
)
(
i
+
1

k
)
α

1
<
(
i
+
1
)
α

1
,
which is established in [5, Lemma 3.2]. ∎
The next result, which is a variant of [10, Lemma 4.3], bounds a weighted sum of the
σ
n

j
that will be used in the proof of Theorem 4.
Lemma 3.
Let the parameter β satisfy
β
>
1
. Then, for
n
=
1
,
2
,
…
,
N
, one has
τ
α
∑
j
=
1
n
j

β
σ
n

j
≤
C
N

1
T
t
n
α

1
.
Proof.
By Lemma 2, we have
(4.7)
τ
α
∑
j
=
1
n
j

β
σ
n

j
≤
τ
α
∑
j
=
1
n
j

β
(
n
+
1

j
)
α

1
≤
τ
α
[
(
n
2
)
α

1
∑
j
=
1
⌈
n
2
⌉
j

β
+
(
n
2
)

(
β

α
)
∑
j
=
⌈
n
2
⌉
+
1
n
j

α
(
n
+
1

j
)
α

1
]
.
But for
s
≥
j

1
and
j
≤
n
, one has
(
n
+
1

s
)
α

1
≥
(
n
+
2

j
)
α

1
≥
2
α

1
(
n
+
1

j
)
α

1
.
Hence,
∑
j
=
⌈
n
2
⌉
+
1
n
j

α
(
n
+
1

j
)
α

1
≤
∑
j
=
⌈
n
2
⌉
+
1
n
(
n
+
1

j
)
α

1
∫
s
=
j

1
j
s

α
𝑑
s
≤
∑
j
=
⌈
n
2
⌉
+
1
n
2
1

α
∫
s
=
j

1
j
(
n
+
1

s
)
α

1
s

α
𝑑
s
≤
2
1

α
∫
s
=
0
n
+
1
(
n
+
1

s
)
α

1
s

α
𝑑
s
=
2
1

α
Γ
(
α
)
Γ
(
1

α
)
,
by [1, Theorem D.6]. Substituting this inequality into (4.7) and using
t
n
=
n
τ
and
β
>
1
, we get
τ
α
∑
j
=
1
n
j

β
σ
n

j
≤
C
τ
α
n
α

1
∑
j
=
1
⌈
n
2
⌉
j

β
+
C
τ
α
n

(
β

α
)
≤
C
τ
α
n
α

1
+
C
τ
α
n

(
β

α
)
=
C
t
n
α

1
(
τ
+
τ
n

(
β

1
)
)
≤
C
t
n
α

1
N

1
T
.
This completes the proof. ∎
We can now prove our main result.
Theorem 4.
Assume that u satisfies (2.3). Then, for
n
=
1
,
2
,
3
,
…
,
N
, the solution
u
m
n
of scheme (3.2) satisfies
(4.8)
max
0
≤
m
≤
M

u
(
x
m
,
t
n
)

u
m
n

≤
C
(
T
α
h
2
+
T
N

1
t
n
α

1
)
for some constant C.
Proof.
Fix
(
x
m
,
t
n
)
∈
Q
. By (4.1) and Lemma 1, the truncation error at
(
x
m
,
t
n
)
satisfies
∥
(
L
M
,
N
(
u
(
x
m
,
t
n
)
)
n

u
n
)
∥
≤
C
(
h
2
+
n

min
{
2

α
,
1
+
α
}
)
.
By (4.5) we then obtain
max
0
≤
m
≤
M

u
(
x
m
,
t
n
)

u
m
n

≤
C
τ
α
∑
j
=
1
n
h
2
σ
n

j
+
C
τ
α
∑
j
=
1
n
j

min
{
2

α
,
1
+
α
}
σ
n

j
.
Invoking Lemma 3 (with
β
=
min
{
2

α
,
1
+
α
}
) for the
j

min
{
2

α
,
1
+
α
}
term and [10, Lemma 4.3] (with
β
=
0
) for the term involving
h
2
, we obtain (4.8). ∎
The bound in (4.8) implies that for any fixed
κ
>
0
one has
(4.9)
max
(
x
m
,
t
n
)
∈
Q
¯
∩
{
t
n
≥
κ
>
0
}

(
x
m
,
t
n
)

u
m
n

≤
C
T
α
(
h
2
+
N

1
)
.
That is, on any subdomain that is bounded away from
t
=
0
, we observe an improved rate of convergence in time compared with the rate of convergence (in time) of
N

α
on
Q
¯
that is given by (3.3).
5 Numerical Results
In this section we give numerical results for the numerical method (3.2) applied to two particular examples from the problem class (2.1). In the first example the exact solution of the problem is known; in the second example it is unknown, so we estimate the order of convergence using the doublemesh principle [2]. In these numerical experiments we always take
N
=
M
. Hence the bounds in (3.3) and (4.8) imply that the spatial error term
C
h
2
will be dominated by the temporal error term
C
N

α
or
C
N

1
.
Example 5.1.
Consider the constant coefficient homogeneous problem
D
t
α
u

∂
2
u
∂
x
2
+
∂
u
∂
x
=
0
for
(
x
,
t
)
∈
Q
=
(
0
,
π
)
×
(
0
,
1
]
,
with initial condition
u
(
x
,
0
)
=
e
x
2
sin
x
,
0
<
x
<
π
, and boundary conditions
u
(
0
,
t
)
=
u
(
π
,
t
)
=
0
,
0
≤
t
≤
1
. The exact solution of this problem is
u
(
x
,
t
)
=
E
α
(

1.25
t
α
)
e
x
2
sin
x
,
where
E
α
is the MittagLeffler function which is defined [1] by
E
α
(
z
)
:=
∑
k
=
0
∞
z
k
Γ
(
α
k
+
1
)
.
In Figure 1 we display the computed solutions with scheme (3.2) for
α
=
0.4
,
0.8
and
N
=
M
=
32
and we observe that the solution has an initial layer at
t
=
0
, which becomes sharper as the parameter α decreases.
Figure 1Example 5.1: Computed solutions with scheme (3.2) for
N
=
M
=
32
.
For Example 5.1 we computed the maximum errors
e
M
,
N
:=
max
(
x
m
,
t
n
)
∈
Q
′

u
(
x
m
,
t
n
)

u
m
n

and the orders of convergence
p
M
,
N
:=
log
2
(
e
M
,
N
e
2
M
,
2
N
)
,
where
Q
′
can be the entire domain
Q
¯
or the subdomain
Q
¯
*
:=
[
0
,
π
]
×
[
0.1
,
1
]
. The numerical results in
Q
¯
(see Table 1) show that scheme (3.2) is
O
(
N

α
)
convergent there (which agrees with [10, Theorem 5.2]), while it is
O
(
N

1
)
convergent in the subdomain
Q
¯
*
(see Table 2), which indicates that the error bound (4.9) is sharp.
Considering the convergence in time, identified by the factor
t
n
α

1
in the error bound (4.8), the error
max
m

u
(
x
m
,
t
n
)

u
m
n

is compared with the error bound
N

1
t
n
α

1
in Figure 2, for
α
=
0.2
and
α
=
0.8
. These plots indicate that the exponent
α

1
in the error bound is sharp for small values of α. However, in the case of larger values of α close to one, the maximum error decreases at a faster rate than
α

1
, as
t
n
increases.
Table 1
Example 5.1: Maximum errors and orders of convergence for scheme (3.2) in the domain
Q
¯
.
α 
N
=
M
=
128

N
=
M
=
256

N
=
M
=
512

N
=
M
=
1024

N
=
M
=
2048

0.4 
8.438E2 
6.714E2 
5.282E2 
4.120E2 
3.191E2 

0.330 
0.346 
0.359 
0.368 

0.6 
3.759E2 
2.512E2 
1.672E2 
1.109E2 
7.342E3 

0.581 
0.588 
0.592 
0.595 

0.8 
1.121E2 
6.401E3 
3.666E3 
2.102E3 
1.206E3 

0.809 
0.804 
0.803 
0.802 

Table 2
Example 5.1: Maximum errors and orders of convergence for scheme (3.2) in the subdomain
Q
¯
*
.
α 
N
=
M
=
128

N
=
M
=
256

N
=
M
=
512

N
=
M
=
1024

N
=
M
=
2048

0.4 
1.024E2 
4.966E3 
2.436E3 
1.214E3 
6.050E4 

1.044 
1.027 
1.005 
1.005 

0.6 
1.300E2 
6.432E3 
3.190E3 
1.595E3 
7.965E4 

1.015 
1.012 
1.000 
1.002 

0.8 
9.844E3 
5.123E3 
2.644E3 
1.361E3 
6.963E4 

0.942 
0.954 
0.959 
0.966 

Figure 2Example 5.1: Loglog plot of the error bound
N

1
t
n
α

1
(
⋄
) and the maximum errors
max
0
≤
m
≤
M

u
(
x
m
,
t
n
)

u
m
n

(
∘
) at each time level
t
=
t
n
,
n
=
1
,
2
,
…
,
N
, generated by scheme (3.2) for
N
=
M
=
100
.
Example 5.2.
Consider the variable coefficient inhomogeneous problem
(5.1)
(5.1a)
D
t
α
u

∂
2
u
∂
x
2
+
(
1
+
x
2
)
∂
u
∂
x
+
(
1
+
x
)
u
=
4
π
2
x
(
π

x
)
for
(
x
,
t
)
∈
Q
,
with
Q
=
(
0
,
π
)
×
(
0
,
1
]
and
(5.1b)
u
(
x
,
0
)
=
0
for
0
<
x
<
π
,
u
(
0
,
t
)
=
u
(
π
,
t
)
=
0
for
0
≤
t
≤
1
.
Figure 3 displays the computed solution for
α
=
0.4
,
0.8
and
N
=
M
=
64
and we observe that the solution again exhibits an initial layer at
t
=
0
.
Figure 3Example 5.2: Computed solutions with scheme (3.2) for
N
=
M
=
64
.
The exact solution of Example 5.2 is unknown and we shall estimate the order of convergence using the twomesh principle [2]. Let
u
m
n
be the computed solution with scheme (3.2) on the mesh
{
(
x
m
,
t
n
)
}
for
m
=
0
,
1
,
…
,
M
,
n
=
0
,
1
,
…
,
N
. To estimate the order of convergence, we compute a new approximation
z
m
/
2
n
2
using the same scheme defined on the finer mesh
{
(
x
m
/
2
,
t
n
/
2
)
}
for
m
=
0
,
1
,
…
,
2
M
,
n
=
0
,
1
,
…
,
2
N
, where
x
m
+
1
/
2
:=
1
2
(
x
m
+
1
+
x
m
)
and
t
n
+
1
/
2
:=
1
2
(
t
n
+
1
+
t
n
)
/
2
. We then compute the twomesh differences
d
M
,
N
:=
max
(
x
m
,
t
n
)
∈
Q
′

u
m
n

z
m
n

and hence the estimated orders of convergence
q
M
,
N
:=
log
2
(
d
M
,
N
d
2
M
,
2
N
)
.
Tables 3 and 4 give the maximum twomesh differences and their corresponding orders of convergence for Example 5.2 in the domain
Q
¯
and the subdomain
Q
¯
*
. The numerical results in both cases are again in agreement with Theorem 4: the order of convergence improves from
O
(
N

α
)
on
Q
¯
to
O
(
N

1
)
on
Q
¯
*
.
Table 3
Example 5.2: Maximum twomesh differences and orders of convergence for scheme (3.2) in the domain
Q
¯
.
α 
N
=
M
=
128

N
=
M
=
256

N
=
M
=
512

N
=
M
=
1024

N
=
M
=
2048

0.4 
1.031E2 
8.673E3 
7.123E3 
5.740E3 
4.558E3 

0.250 
0.284 
0.311 
0.333 

0.6 
4.935E3 
3.338E3 
2.234E3 
1.486E3 
9.857E4 

0.564 
0.579 
0.588 
0.593 

0.8 
1.661E3 
9.441E4 
5.368E4 
3.060E4 
1.748E4 

0.815 
0.815 
0.811 
0.808 

Table 4
Example 5.2: Maximum twomesh differences and orders of convergence for scheme (3.2) in the subdomain
Q
¯
*
.
α 
N
=
M
=
128

N
=
M
=
256

N
=
M
=
512

N
=
M
=
1024

N
=
M
=
2048

0.4 
5.849E4 
2.783E4 
1.351E4 
6.711E5 
3.337E5 

1.072 
1.042 
1.010 
1.008 

0.6 
1.148E3 
5.457E4 
2.628E4 
1.291E4 
6.356E5 

1.073 
1.054 
1.025 
1.023 

0.8 
1.335E3 
6.752E4 
3.387E4 
1.703E4 
8.531E5 

0.984 
0.995 
0.992 
0.997 

In [10], numerical results were given for the particular case of a fractional reactiondiffusion equation (i.e., with
q
≡
0
in (2.1)) showing that the scheme also converges with order α when the whole domain is considered. Additional numerical results that illustrate the improved rate of convergence away from
t
=
0
are given in [3].
Funding source: National Natural Science Foundation of China
Award Identifier / Grant number: 91430216
Award Identifier / Grant number: NSAF U1530401
Funding source: Ministerio de Economía y Competitividad
Award Identifier / Grant number: MTM201675139R
References
[1] K. Diethelm, The Analysis of Fractional Differential Equations, Lecture Notes in Math. 2004, Springer, Berlin, 2010. Search in Google Scholar
[2] P. A. Farrell, A. F. Hegarty, J. J. H. Miller, E. O’Riordan and G. I. Shishkin, Robust Computational Techniques for Boundary Layers, Appl. Math. Math. Comput. 16, Chapman & Hall/CRC, Boca Raton, 2000. Search in Google Scholar
[3] J. L. Gracia, E. O’Riordan and M. Stynes, Convergence outside the initial layer for a numerical method for the timefractional heat equation, Numerical Analysis and Its Applications, Lecture Notes in Comput. Sci. 10187, Springer, Cham (2017), 82–94. Search in Google Scholar
[4] B. Jin, R. Lazarov and Z. Zhou, An analysis of the L1 scheme for the subdiffusion equation with nonsmooth data, IMA J. Numer. Anal. 36 (2016), no. 1, 197–221. Search in Google Scholar
[5] B. Jin and Z. Zhou, An analysis of Galerkin proper orthogonal decomposition for subdiffusion, ESAIM Math. Model. Numer. Anal. 51 (2017), no. 1, 89–113. Search in Google Scholar
[6] W. McLean and K. Mustapha, Timestepping error bounds for fractional diffusion problems withnonsmooth initial data, J. Comput. Phys. 293 (2015), 201–217. Search in Google Scholar
[7] K. Mustapha, An implicit finitedifference timestepping method for a subdiffusion equation, with spatial discretization by finite elements, IMA J. Numer. Anal. 31 (2011), no. 2, 719–739. Search in Google Scholar
[8] K. Mustapha and W. McLean, Piecewiselinear, discontinuous Galerkin method for a fractional diffusion equation, Numer. Algorithms 56 (2011), no. 2, 159–184. Search in Google Scholar
[9] M. Stynes, Too much regularity may force too much uniqueness, Fract. Calc. Appl. Anal. 19 (2016), no. 6, 1554–1562. Search in Google Scholar
[10] M. Stynes, E. O’Riordan and J. L. Gracia, Error analysis of a finite difference method on graded meshes for a timefractional diffusion equation, SIAM J. Numer. Anal. 55 (2017), no. 2, 1057–1079. Search in Google Scholar
Received: 2017127
Revised: 2017527
Accepted: 201766
Published Online: 201776
Published in Print: 201811
© 2017 Walter de Gruyter GmbH, Berlin/Boston