See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/327790093

Subdiffusion with a time-dependent coefficient: analysis and numerical

solution

Preprint · September 2018

CITATIONS READS

0 714

3 authors:

Bangti Jin Buyang Li

University College London The Hong Kong Polytechnic University
233 PUBLICATIONS 6,798 CITATIONS 130 PUBLICATIONS 2,873 CITATIONS

SEE PROFILE SEE PROFILE

Zhi Zhou
The Hong Kong Polytechnic University
113 PUBLICATIONS 3,409 CITATIONS

SEE PROFILE

All content following this page was uploaded by Buyang Li on 01 December 2018.

The user has requested enhancement of the downloaded file.

 SUBDIFFUSION WITH A TIME-DEPENDENT COEFFICIENT: ANALYSIS AND
NUMERICAL SOLUTION

BANGTI JIN, BUYANG LI , AND ZHI ZHOU

Abstract. In this work, a complete error analysis is presented for fully discrete solutions of the sub-
diffusion equation with a time-dependent diffusion coefficient, obtained by the Galerkin finite element
method with conforming piecewise linear finite elements in space and backward Euler convolution quad-
rature in time. The regularity of the solutions of the subdiffusion model is proved for both nonsmooth
initial data and incompatible source term. Optimal-order convergence of the numerical solutions is
established using the proven solution regularity and a novel perturbation argument via freezing the
diffusion coefficient at a fixed time. The analysis is supported by numerical experiments.

1. Introduction
d
Let Ω ⊂ R (d ≥ 1) be a convex polyhedral domain with a boundary ∂Ω. Consider the following
fractional-order parabolic problem for the function u(x, t):
 α
 ∂t u(x, t) − ∇ · (a(x, t)∇u(x, t)) = f (x, t) (x, t) ∈ Ω × (0, T ],

(1.1) u(x, t) = 0 (x, t) ∈ ∂Ω × (0, T ],

u(x, 0) = u0 (x) x ∈ Ω,


where T > 0 is a fixed final time, f ∈ L∞ (0, T ; L2 (Ω)) and u0 ∈ L2 (Ω) are given source term and initial
data, respectively, and a(x, t) ∈ Rd×d is a symmetric matrix-valued diffusion coefficient such that for
some constant λ ≥ 1
(1.2) λ−1 |ξ|2 ≤ a(x, t)ξ · ξ ≤ λ|ξ|2 , ∀ ξ ∈ Rd , ∀ (x, t) ∈ Ω × (0, T ],
(1.3) |∂t a(x, t)| + |∇x a(x, t)| + |∇x ∂t a(x, t)| ≤ c, ∀ (x, t) ∈ Ω × (0, T ].
The notation ∂tα u(t) denotes the Caputo derivative in time of order α ∈ (0, 1), defined by [16, p. 70]
Z t
1
(1.4) ∂tα u(x, t) = (t − s)−α ∂s u(x, s)ds.
Γ(1 − α) 0
The literature on the numerical analysis of the subdiffusion problem is vast; see [21, 11, 9, 15] for a
rather incomplete list and the overview [10] (and the references therein). The work [11] analyzed two
spatially semidiscrete schemes, i.e., Galerkin finite element method (FEM) and lumped mass method,
and derived nearly optimal order error estimates for the homogeneous problem. The inhomogeneous case
was analyzed in [9]. See [15] for a finite volume element discretization, and [14] for a unified approach.
There are a number of fully discrete schemes, e.g., convolution quadrature [36, 12], piecewise polynomial
interpolation [33, 21, 2, 7, 35, 32], discontinuous Galerkin method [26, 27]; and some of them have an
O(τ ) rate for nonsmooth data, with τ being time step size. However, all these works analyzed only the
case that the diffusion coefficient a is independent of the time t. These works mostly employ Laplace
transform and its discrete analogue for analysis, which are not directly applicable to the case of a time-
dependent coefficient. Recently, Mustapha [28] analyzed the spatially semidiscrete Galerkin FEM for
(1.1) using a novel energy argument, and proved optimal-order convergence rates for both smooth and
nonsmooth initial data (with a zero source term) based on certain assumptions on the regularity of the
PDE’s solution.

2010 Mathematics Subject Classification. Primary: 65M30, 65M15, 65M12.

Key words and phrases. subdiffusion, time-dependent coefficient, Galerkin finite element method, convolution quadra-
ture, perturbation argument, error estimate.
1
 2
In this article, using a novel perturbation argument, we present a new approach to analyze a fully
discrete scheme for problem (1.1) based on the Galerkin FEM in space and backward Euler (BE) convo-
lution quadrature in time, covering initial data and source term simultaneously. The main contributions
of this paper are as follows. First, we give a complete existence, uniqueness and regularity theory for
problem (1.1) in Theorems 2.1–2.3, which are crucial to the error analysis. Second, we derive sharp error
estimates for the spatially semidiscrete Galerkin FEM. This is achieved by combining error estimates for
a time-independent coefficient and a perturbation argument in time. Third, we derive nearly sharp error
estimates for the fully discrete method. All error estimates are given directly in terms of the regularity
of the initial data and source term, under mild regularity assumptions on the diffusion coefficient a(x, t)
that are weaker than the assumptions in [28]; see Remark 2.2 for the precise statement.
There are a few relevant works on standard parabolic problems with a time-dependent coefficient
[24, 30, 31, 19]. For example, Luskin and Rannacher [24] proved optimal order error estimates for both
spatially semidiscrete and fully discrete problems (by BE method) using a novel energy argument, and
Sammon [30] analyzed fully discrete schemes with linear multistep methods. Our error analysis relies
crucially on a perturbation argument, using basic estimates given in Lemmas 3.1 and 3.2, which are of
independent interest. Generally, the idea of freezing coefficients and perturbation in time has been proved
very useful in combination with energy estimates [31] and Lp estimates [1, 18, 20]. In this work, we have
successfully adapted the idea to the subdiffuion model.
The rest of the paper is organized as follows. In Section 2, we discuss temporal and spatial regularity
of the solution for nonsmooth problem data. Then in Section 3, we prove optimal-order convergence of
the spatially semidiscrete Galerkin FEM for both homogeneous and inhomogeneous problems. In Section
4, we present the error analysis for the fully discrete FEM and prove first-order convergence in time.
Last, in Section 5, we present numerical examples to support the theoretical analysis. Throughout, the
notation c, with or without a subscript, denotes a generic positive constant, which may differ at each
occurrence, but is always independent of the mesh size h and step size τ .

2. Regularity theory
In this section we investigate the regularity of the solutions of problem (1.1). For any function f (x, t)
defined on Ω × (0, T ), we denote by f (t) the function f (·, t). Let −∆ : H01 (Ω) ∩ H 2 (Ω) → L2 (Ω) be the
negative Laplacian operator with a zero Dirichlet boundary condition, and {(λj , ϕj )} be its eigenvalues
ordered nondecreasingly (with multiplicity counted) and the corresponding eigenfunctions normalized in
r
the L2 (Ω) norm. For any r ≥ 0, we denote the space Ḣ r (Ω) = {v ∈ L2 (Ω) : (−∆) 2 v ∈ L2 (Ω)}, with the
norm [34, Chapter 3]
∞
X
kvk2Ḣ r (Ω) = λrj (v, ϕj )2 .
j=1
0 2 1
Then we have Ḣ (Ω) = L (Ω), Ḣ (Ω) = H01 (Ω), and Ḣ 2 (Ω) = H 2 (Ω) ∩ H01 (Ω).

2.1. Subdiffusion with a time-independent coefficient. First we recall basic properties of the sub-
diffusion model with a time-independent diffusion coefficient, i.e., a(x, t) = a(x). Accordingly, we denote
by A : H01 (Ω) ∩ H 2 (Ω) → L2 (Ω) an elliptic operator, defined by
Aφ(x) := −∇ · (a(x)∇φ(x)),
and consider the problem
(2.1) ∂tα u(t) + Au(t) = f (t) t ∈ (0, T ], with u(0) = u0 .
This problem has been studied in [3, 4, 23, 25, 29, 13]. The following maximal Lp -regularity holds [3].
Lemma 2.1. If u0 = 0 and f ∈ Lp (0, T ; L2 (Ω)) with 1 < p < ∞, then problem (2.1) has a unique
solution u ∈ Lp (0, T ; Ḣ 2 (Ω)) such that ∂tα u ∈ Lp (0, T ; L2 (Ω)) and
kukLp (0,T ;Ḣ 2 (Ω)) + k∂tα ukLp (0,T ;L2 (Ω)) ≤ ckf kLp (0,T ;L2 (Ω)) ,
where the constant c does not depend on f and T .
 3
By means of Laplace transform, the solution u(t) can be represented by [13, Section 4]
Z t
(2.2) u(t) = F (t)u0 + E(t − s)f (s)ds,
0
where the solution operators F (t) and E(t) are defined by
Z
1
(2.3) F (t) := ezt z α−1 (z α + A)−1 dz,
2πi Γθ,δ
Z
1
(2.4) E(t) := ezt (z α + A)−1 dz,
2πi Γθ,δ
with integration over a contour Γθ,δ ⊂ C (oriented with an increasing imaginary part):
Γθ,δ = {z ∈ C : |z| = δ, | arg z| ≤ θ} ∪ {z ∈ C : z = ρe±iθ , ρ ≥ δ}.
Throughout, we fix θ ∈ ( π2 , π) so that z α ∈ Σαθ ⊂ Σθ := {0 6= z ∈ C : arg(z) ≤ θ}, for all z ∈ Σθ . The
next lemma gives smoothing properties of F (t) and E(t), which follow from the resolvent estimate
(2.5) k(z + A)−1 k ≤ cφ |z|−1 , ∀z ∈ Σφ , ∀ φ ∈ (0, π),
where k · k denotes the operator norm from L (Ω) to L2 (Ω).
2

Lemma 2.2. The operators F and E defined in (2.3) and (2.4) satisfy the following properties.
(i) t−α kA−1 (F (t) − I)k + kF (t) − Ik ≤ c, ∀ t ∈ (0, T ];
(ii) t1−α kE(t)k + t2−α kE 0 (t)k + tkAE(t)k ≤ c, ∀ t ∈ (0, T ];
(iii) tα kAF (t)k + t1−βα kA−β F 0 (t)k ≤ c, ∀ t ∈ (0, T ], β ∈ [0, 1].
Proof. Parts (i) and (ii) have been shown in [13, Lemma 3.4]. By letting δ = t−1 in Γθ,δ and z =
s cos ϕ + is sin ϕ, using (2.5), we have (with |dz| being the arc length of Γθ,δ )
Z Z
1
kAF (t)k = ezt z α−1 A(z α + A)−1 dz ≤ c e<(z)t |z|α−1 |dz|
2πi Γθ,δ Γθ,δ
Z ∞ Z θ
≤c est cos θ sα−1 ds + c ecos ϕ δ α dϕ ≤ ct−α .
δ −θ
0 1
e z (z α + A) dz, and thus by the estimate (2.5),
zt α
R
Next, direct computation gives F (t) = 2πi Γθ,δ
Z
kF 0 (t)k ≤ c e<(z)t |z|α |z|−α |dz| ≤ ct−1 ,
Γθ,δ

which shows the Rassertion for β = 0. Meanwhile, R by ztthe identity z α (z α + A)−1 = I − A(z α + A)−1 , we
0 1 zt α α 1 α
have F (t) = 2πi Γθ,δ e z (z + A) dz = − 2πi Γθ,δ e A(z + A) dz, and thus
Z
kA−1 F 0 (t)k ≤ c e<(z)t |z|−α |dz| ≤ ct−1+α .
Γθ,δ

This shows the assertion for β = 1. Then the desired bound on t1−βα kA−β F 0 (t)k in part (iii) follows by
interpolation. 
2.2. Regularity theory for subdiffusion with a time-dependent coefficient. Now we develop
the regularity theory for the case of a time-dependent diffusion coefficient. The work [37] gave some
interior Hölder estimates for bounded measurable coefficients. Recently, Kubica et al [17] showed the
unique existence by approximating the coefficients by smooth functions, and derived several regularity
estimates. We shall provide a different approach to derive regularity estimates in Sobolev spaces, which
are essential for the error analysis in Sections 3 and 4.
We define a time-dependent elliptic operator A(t) : Ḣ 2 (Ω) → L2 (Ω) by
A(t)φ = −∇ · (a(x, t)∇φ), ∀φ ∈ Ḣ 2 (Ω).
Under condition (1.3), the following estimate holds:
(2.6) k(A(t) − A(s))vkL2 (Ω) ≤ c|t − s|kvkH 2 (Ω) .
 4
First we give the existence, uniqueness and regularity of solutions to problem (1.1) with u0 = 0.
Theorem 2.1. Under conditions (1.2)-(1.3), with u0 = 0 and f ∈ Lp (0, T ; L2 (Ω)), 1/α < p < ∞, prob-
lem (1.1) has a unique solution u ∈ C([0, T ]; L2 (Ω)) ∩ Lp (0, T ; Ḣ 2 (Ω)) such that ∂tα u ∈ Lp (0, T ; L2 (Ω)).
Proof. For any θ ∈ [0, 1], consider the following fractional-order parabolic problem
(2.7) ∂tα u(t) + A(θt)u(t) = f (t), t ∈ (0, T ], with u(0) = 0,
and define a set
D = {θ ∈ [0, 1] : (2.7) has a solution u ∈ Lp (0, T ; Ḣ 2 (Ω)) such that ∂tα u ∈ Lp (0, T ; L2 (Ω))}.
Lemma 2.1 implies 0 ∈ D and so D 6= ∅.
For any θ ∈ D, by rewriting (2.7) as
(2.8) ∂tα u(t) + A(θt0 )u(t) = f (t) + (A(θt0 ) − A(θt))u(t), t ∈ (0, T ], with u(0) = 0,
and by applying Lemma 2.1 in the time interval (0, t0 ), we obtain
k∂tα ukLp (0,t0 ;L2 (Ω)) + kukLp (0,t0 ;H 2 (Ω))
≤ckf kLp (0,t0 ;L2 (Ω)) + ck(A(θt0 ) − A(θt))u(t)kLp (0,t0 ;L2 (Ω))
(2.9) ≤ckf kLp (0,t0 ;L2 (Ω)) + ck(t0 − t)u(t)kLp (0,t0 ;H 2 (Ω)) ,
where the last line follows from (2.6). Let g(t) = kukpLp (0,t;H 2 (Ω)) , which satisfies g 0 (t) = ku(t)kpH 2 (Ω) .
Then (2.9) and integration by parts imply
Z t0
g(t0 ) ≤ ckf kpLp (0,t0 ;L2 (Ω)) + c (t0 − t)p g 0 (t)dt
0
Z t0
p
= ckf kLp (0,t0 ;L2 (Ω)) + cp (t0 − t)p−1 g(t)dt
0
Z t0
p
≤ ckf kLp (0,t0 ;L2 (Ω)) + c g(t)dt,
0
which implies (via the standard Gronwall’s inequality)
g(t0 ) ≤ ckf kpLp (0,t0 ;L2 (Ω)) , i.e., kukLp (0,t0 ;H 2 (Ω)) ≤ ckf kLp (0,t0 ;L2 (Ω)) .
Substituting the last inequality into (2.9) yields
(2.10) k∂tα ukLp (0,t0 ;L2 (Ω)) + kukLp (0,t0 ;H 2 (Ω)) ≤ ckf kLp (0,t0 ;L2 (Ω)) .
Since the estimate (2.10) is independent of θ ∈ D, D is a closed subset of [0, 1].
Now we show that D is also open with respect to the subset topology of [0, 1]. In fact, if θ0 ∈ D, then
problem (2.7) can be rewritten as
(2.11) ∂tα u(t) + A(θ0 t)u(t) + (A(θt) − A(θ0 t))u(t) = f (t), t ∈ (0, T ], with u(0) = 0,
which is equivalent to
h i
(2.12) 1 + (∂tα + A(θ0 t))−1 (A(θt) − A(θ0 t)) u(t) = (∂tα + A(θ0 t))−1 f (t).

It follows from (2.10) that the operator (∂tα + A(θ0 t))−1 (A(θt) − A(θ0 t)) is small in the sense that
k(∂tα + A(θ0 t))−1 (A(θt) − A(θ0 t))kLp (0,T ;H 2 (Ω))→Lp (0,T ;H 2 (Ω)) ≤ c|θ − θ0 |.
Thus for θ sufficiently close to θ0 , the operator 1 + (∂tα + A(θ0 t))−1 (A(θt) − A(θ0 t)) is invertible on
Lp (0, T ; Ḣ 2 (Ω)), which implies θ ∈ D. Thus D is open with respect to the subset topology of [0, 1]. Since
D is both closed and open respect to the subset topology of [0, 1], D = [0, 1]. Further, note that for
1/α < p < ∞, the inequality (2.10) and the condition u(0) = 0 directly imply u ∈ C([0, T ]; L2 (Ω)) [13,
Theorem 2.1], which completes the proof of the theorem. 
The following generalized Gronwall’s inequality is useful ([5, Lemma 6.3] and [8, Exercise 3, p. 190]).
 5
Lemma 2.3. Let the function ϕ(t) ≥ 0 be continuous for 0 < t ≤ T . If
Z t
−1+α
ϕ(t) ≤ at +b (t − s)−1+β ϕ(s)ds, 0 < t ≤ T,
0
for some constants a, b ≥ 0, α, β > 0, then there is a constant c = c(b, T, α, β) such that
ϕ(t) ≤ cat−1+α , 0 < t ≤ T.
Next we give the spatial regularity of the solution u for the case f = 0.
Theorem 2.2. Under conditions (1.2)-(1.3), with u0 ∈ Ḣ β (Ω), 0 ≤ β ≤ 2, and f = 0, problem (1.1) has
a unique solution u ∈ C([0, T ]; L2 (Ω)) ∩ C((0, T ]; Ḣ 2 (Ω)) such that ∂tα u ∈ C((0, T ]; L2 (Ω)), and
ku(t)kH 2 (Ω) ≤ ct−(1−β/2)α ku0 kḢ β (Ω) .
Proof. The existence and uniqueness of a solution can be proved in the same way as Theorem 2.1 based
on the a priori estimate below. We rewrite problem (1.1) as
∂tα u(t) + A(t0 )u(t) = (A(t0 ) − A(t))u(t) + f (t), t ∈ (0, T ], with u(0) = u0 ,
Then the solution u(t) can be represented by
Z t Z t
(2.13) u(t) = F (t; t0 )u0 + E(t − s; t0 )(A(t0 ) − A(s))u(s)ds + E(t − s; t0 )f (s)ds,
0 0
where the operators F (t; t0 ) and E(t; t0 ) are defined respectively by
Z Z
1 −1 1
F (t; t0 ) := zt α−1 α
e z (z + A(t0 )) dz and E(t; t0 ) := ezt (z α + A(t0 ))−1 dz.
2πi Γθ,δ 2πi Γθ,δ
In the case f = 0, applying A(t0 ) to both sides of (2.13) yields
Z t
A(t0 )u(t) =A(t0 )F (t; t0 )u0 + A(t0 )E(t − s; t0 )(A(t0 ) − A(s))u(s)ds.
0
Then conditions (1.2)-(1.3) and Lemma 2.2(ii) imply
Z t0
ku(t0 )kH 2 (Ω) ≤ ckA(t0 )F (t0 ; t0 )u0 kL2 (Ω) + c kA(t0 )E(t0 − s; t0 )kk(A(t0 ) − A(s))u(s)kL2 (Ω) ds
0
Z t0
−(1−β/2)α
≤ ct0 ku0 kḢ β (Ω) + c (t0 − s)−1 (t0 − s)ku(s)kH 2 (Ω) ds
0
Z t0
−(1−β/2)α
= ct0 ku0 kḢ β (Ω) + c ku(s)kH 2 (Ω) ds, ∀ t0 ∈ (0, T ].
0
The desired estimate follows from the generalized Gronwall’s inequality in Lemma 2.3. It remains to
show u ∈ C([0, T ]; L2 (Ω)) ∩ C((0, T ]; H 2 (Ω)). Indeed, note that (by fixing t0 = 0)
Z t
ku(t) − u0 kL2 (Ω) ≤ kF (t; 0)u0 − u0 kL2 (Ω) + (t − s)−α sku(s)kH 2 (Ω) ds,
0
which together with the bound on ku(s)kH 2 (Ω) implies
lim ku(t) − u0 kL2 (Ω)
t→0+
Z t
≤ lim+ kF (t; 0)u0 − u0 kL2 (Ω) + lim+ c (t − s)α−1 s1−(1−β/2)α dsku0 kḢ β (Ω) = 0.
t→0 t→0 0
i.e., limt→0+ u(t) = u0 in L2 (Ω). The rest of the assertion follows similarly. This completes the proof. 
To analyze the temporal regularity, we first give three technical lemmas.
Lemma 2.4. Let conditions (1.2) and (1.3) be fulfilled, and u be the solution to problem (1.1) with
u0 ∈ L2 (Ω) and f = 0. Then there holds
d t
Z
k E(t − s; t0 )(A(t0 ) − A(s))u(s)ds t=t kL2 (Ω) ≤ cku0 kL2 (Ω) .
dt 0 0
 d
Rt 6
Proof. Let I = dt 0
E(t − s; t0 )(A(t0 ) − A(s))u(s)ds|t=t0 . Then
 Z t0 +ε
1
I = lim E(t0 + ε − s; t0 )(A(t0 ) − A(s))u(s)ds
ε→0 ε 0
Z t0 
(2.14) − E(t0 − s; t0 )(A(t0 ) − A(s))u(s)ds =: lim Λ(ε).
0 ε→0

If ε > 0, then
Z t0 +ε
1
Λ(ε) = E(t0 + ε − s; t0 )(A(t0 ) − A(s))u(s)ds
ε t0
Z t0
E(t0 + ε − s; t0 ) − E(t0 − s; t0 )
(2.15) + (A(t0 ) − A(s))u(s)ds =: I+ + II+ .
0 ε
By applying Lemma 2.2(ii), (2.6) and Theorem 2.2, we deduce
Z t0 +ε
kI+ kL2 (Ω) ≤ cε−1 kE(t0 + ε − s; t0 )kk(A(t0 ) − A(s))u(s)kL2 (Ω) ds
t0
Z t0 +ε
≤ cε−1 |t0 + ε − s|α−1 |t0 − s|ku(s)kH 2 (Ω) ds
t0
Z t0 +ε
≤ cε−1 |t0 + ε − s|α−1 |t0 − s|s−α ku0 kL2 (Ω) ds
t0
Z t0 +ε
≤ cε−1 |t0 + ε − s|α−1 (ε)t−α α −α
0 ku0 kL2 (Ω) ds ≤ cε t0 ku0 kL2 (Ω) ,
t0

and similarly,
Z t0 Z 1
kII+ kL2 (Ω) = E 0 (t0 + θε − s; t0 )(A(t0 ) − A(s))u(s) dθds
0 0 L2 (Ω)
Z 1 Z t0
≤c (t0 + θε − s)α−2 (t0 − s)ku(s)kH 2 (Ω) dsdθ
0 0
Z t0
≤c (t0 − s)α−1 ku(s)kH 2 (Ω) ds
0
Z t0
≤c (t0 − s)α−1 s−α ku0 kL2 (Ω) ds ≤ cku0 kL2 (Ω) .
0

If −t0 < ε < 0, then

Z t0
Λ(ε) = ε−1 E(t0 − s; t0 )(A(t0 ) − A(s))u(s)ds
t0 −|ε|
t0 −|ε|
E(t0 − s; t0 ) − E(t0 − |ε| − s; t0 )
Z
+ (A(t0 ) − A(s))u(s)ds =: I− + II− ,
0 ε
and similarly, we obtain
kI− kL2 (Ω) ≤ cεα (t0 + ε)−α ku0 kL2 (Ω) and kII− kL2 (Ω) ≤ cku0 kL2 (Ω) .
Combining the preceding estimates yields the assertion. 

Lemma 2.5. Let conditions (1.2) and (1.3) be fulfilled, and u be the solution to problem (1.1) with
Rt
f ∈ C([0, T ]; L2 (Ω)), 0 (t − s)α−1 kf 0 (s)kL2 (Ω) ds < ∞ and u0 = 0. Then there holds
Z t Z t
(t − s)α−1 ku(s)kH 2 (Ω) ds ≤ ckf (0)kL2 (Ω) + c (t − s)α−1 kf 0 (s)kL2 (Ω) ds.
0 0
 7
Proof. By the solution representation (2.13) with u0 = 0, we have
Z t0 Z t0
A(t0 )u(t0 ) = A(t0 )E(t0 − s; t0 )f (s)ds + A(t0 )E(t0 − s; t0 )(A(s) − A(t0 ))u(s)ds := I + II.
0 0
It follows directly from the definition of the operators E(s; t0 ) and F (s; t0 ) that the identity A(t0 )E(s; t0 ) =
d d
− ds F (s; t0 ) = ds (I − F (s; t0 )) holds. So upon changing variables and integration by parts, we obtain
Z t0 Z t0
d
I= A(t0 )E(s; t0 )f (t0 − s)ds = (I − F (s; t0 ))f (t0 − s)ds
0 0 ds
Z t0
d
= (F (t0 ; t0 ) − I)f (0) − (I − F (s; t0 )) f (t0 − s)ds,
0 ds
where we have used the identity F (0; t0 ) = I. Thus, by Lemma 2.2(i), we obtain
Z t0
kIkL2 (Ω) ≤ ckf (0)kL2 (Ω) + c kf 0 (s)kL2 (Ω) ds.
0
Similarly, by Lemma 2.2(ii) and (2.6), for the term II, we have
Z t0 Z t0
kIIkL2 (Ω) ≤ c (t0 − s)−1 |t0 − s|kA(t0 )u(s)kL2 (Ω) ds = c kA(t0 )u(s)kL2 (Ω) ds.
0 0
Rt
Let g(t) = 0
(t − s)α−1 ku(s)kH 2 (Ω) ds. Then the last two estimates together give
Z t
g(t) ≤ c (t − s)α−1 (kIkL2 (Ω) + kIIkL2 (Ω) ) ds
0
Z t Z s Z ξ
(t − s)α−1 kf (0)kL2 (Ω) + kf 0 (ξ)kL2 (Ω) dξ +


≤c ku(ξ)kH 2 (Ω) dξ ds
0 0 0
Z t Z t
α α 0
≤ ct kf (0)kL2 (Ω) + c (t − s) kf (s)kL2 (Ω) ds + c g(s) ds,
0 0
where the last line follows directly from the semigroup property of Riemann-Liouville integral and change
of integration orders. Now Gronwall’s inequality gives
Z t
g(t) ≤ ctα kf (0)kL2 (Ω) + c (t − s)α kf 0 (s)kL2 (Ω) ds,
0
from which the desired assertion follows directly. 
Lemma 2.6. Let conditions (1.2) and (1.3) be fulfilled, and u be the solution to problem (1.1) with
Rt
f ∈ C([0, T ]; L2 (Ω)), 0 (t − s)α−1 kf 0 (s)kL2 (Ω) ds < ∞ and u0 = 0. Then there holds
d t
Z Z t0
k E(t − s; t0 )(A(t0 ) − A(s))u(s)ds t=t0 kL2 (Ω) ≤ ckf (0)kL2 (Ω) + c (t0 − s)α−1 kf 0 (s)kL2 (Ω) ds.
dt 0 0

Proof. For any small ε > 0, we employ the splitting (2.14). By Lemma 2.2(ii) and (2.6), we bound the
term I+ by
Z t0 +ε
kI+ kL2 (Ω) ≤ ε−1 kE(t0 + ε − s; t0 )kk(A(t0 ) − A(s))u(s)kL2 (Ω) ds
t0
Z t0 +ε
≤ cε−1 |t0 + ε − s|α−1 |t0 − s|ku(s)kH 2 (Ω) ds
t0
Z t0 +ε
≤c |t0 + ε − s|α−1 ku(s)kH 2 (Ω) ds.
t0
This estimate, Hölder’s inequality and Theorem 2.1 directly imply limε→0+ kI+ kL2 (Ω) = 0. For the term
II+ , Lemma 2.2(ii) gives
Z t0 Z 1
kII+ kL2 (Ω) = E 0 (t0 + θε − s; t0 )(A(t0 ) − A(s))u(s) dθds
0 0 L2 (Ω)
 Z 1 Z t0 8
≤c (t0 + θε − s)α−2 (t0 − s)ku(s)kH 2 (Ω) dsdθ
0 0
Z t0
≤c (t0 − s)α−1 ku(s)kH 2 (Ω) ds,
0

which together with Lemma 2.5 yields the assertion for ε > 0. Similar estimates hold for the case ε < 0,
and this completes the proof of the lemma. 

Remark 2.1. Note that the bound on II+ in Lemma 2.4 blows up for α → 1− :
Z t0
(t0 − s)α−1 s−α ds = B(α, 1 − α),
0

and in view of the asymptotics B(α, 1 − α) = O((1 − α)−1 ) as α → 1− , it blows up at a rate 1/(1 − α).
Actually, this can be avoided by the following alternative argument:
Z t0 Z 1
kII+ kL2 (Ω) = E 0 (t0 + θε − s; t0 )A(t0 )1/2 A(t0 )1/2 (1 − A(t0 )−1 A(s))u(s) dθds
0 0 L2 (Ω)
Z 1 Z t0
≤c (t0 + θε − s)α/2−2 (t0 − s)kA(t0 )1/2 u(s)kL2 (Ω) dsdθ
0 0
Z t0
≤c (t0 − s)α/2−1 ku(s)kH 1 (Ω) ds
0
Z t0
≤c (t0 − s)α/2−1 s−α/2 ku0 kL2 (Ω) ds ≤ cku0 kL2 (Ω) ,
0

where the first inequality is due to (2.6), Corollary 3.1 below and interpolation. The same argument can
be applied to the term II+ in Lemma 2.6. Thus, the involved constants are bounded for α → 1− .
Now we can give the temporal regularity of the solution u.
Theorem 2.3. Let conditions (1.2)-(1.3) be fulfilled, and u be the solution to problem (1.1).
(i) For u0 ∈ Ḣ β (Ω), 0 ≤ β ≤ 2, and f = 0, then
ku0 (t)kL2 (Ω) ≤ ct−(1−αβ/2) ku0 kḢ β (Ω) .
Rt
(ii) For u0 = 0, f ∈ C([0, T ]; L2 (Ω)) and 0 (t − s)α−1 kf 0 (s)kL2 (Ω) ds < ∞, then
Z t
0 −(1−α)
ku (t)kL2 (Ω) ≤ ct kf (0)kL2 (Ω) + c (t − s)α−1 kf 0 (s)kL2 (Ω) ds.
0
p 2
(iii) For u0 = 0 and f ∈ L (0, T ; L (Ω)) with 2/α < p < ∞, then
ku(t)kH 1 (Ω) ≤ ckf kLp (0,t;L2 (Ω)) .
Proof. The proof employs the solution representation (2.13). By Lemma 2.2(iii), we have
d −(1−αβ/2)
F (t; t0 )u0 t=t0 L2 (Ω) ≤ ct0 ku0 kḢ β (Ω) .
dt
This and Lemma 2.4 yield the assertion in part (i).
To show part (ii), differentiating (2.13) with respect to t yields
d t d t
Z Z
(2.16) u0 (t0 ) = E(t − s; t0 )(A(t0 ) − A(s))u(s)ds + E(s; t0 )f (t − s)ds .
dt 0 t=t0 dt 0 t=t0

In view of the identity

Z t Z t
d
E(s; t0 )f (t − s)ds = E(t; t0 )f (0) + E(s; t0 )f 0 (t − s)ds,
dt 0 0
 9
by Lemma 2.2(ii), we have
d t t
Z Z
E(s; t0 )f (t − s)ds ≤ kE(t; t0 )f (0)kL2 (Ω) + kE(s; t0 )f 0 (t − s)kL2 (Ω) ds,
dt 0 L2 (Ω) 0
Z t
−(1−α)
(2.17) ≤ ct kf (0)kL2 (Ω) + c sα−1 kf 0 (t − s)kL2 (Ω) ds.
0
This and Lemma 2.6 complete the proof of part (ii).
Last, for the choice 2/α < p < ∞, Lemma 2.1 implies
u ∈ Lp (0, T ; H 2 (Ω)) ∩ W α,p (0, T ; L2 (Ω)) ,→ W α/2,p (0, T ; (L2 (Ω), H 2 (Ω))1/2 )
= W α/2,p (0, T ; H 1 (Ω)) ,→ C([0, T ]; H 1 (Ω)),
where (L2 (Ω), H 2 (Ω))1/2 denotes the complex interpolation space between L2 (Ω) and H 2 (Ω), and the
last embedding is a consequence of [13, equation (2.3)]. Then the proof of Theorem 2.3 is complete. 
Remark 2.2. In the error analysis, the work [28] requires the following conditions on the coefficient
a(x, t): a(x, t), ∂t a(x, t) ∈ L∞ (0, T ; W 1,∞ (Ω)) and ∂tt
2
a(x, t) ∈ L∞ (0, T ; L∞ (Ω)), which are more strin-
gent than (1.3). Further, the work [28] has assumed the following regularity on the solution u to the
homogeneous problem: for 0 ≤ p ≤ q ≤ 2,
ku(t)kḢ q (Ω) + tku0 (t)kḢ q (Ω) ≤ ct−(q−p)α/2 ku0 kḢ p (Ω) .
In contrast, for the homogeneous problem, we proved the following estimates under assumption (1.3):
t(1−β/2)α ku(t)kH 2 (Ω) + t1−αβ/2 ku0 (t)kL2 (Ω) ≤ cku0 kḢ β (Ω)
and similar estimates for the inhomogeneous problem. It is worth noting that unlike the argument in [28],
the error analysis below does not need the regularity ku0 (t)kḢ 2 (Ω) , which allows us to relax the regularity
assumption on the coefficient a(x, t).
Remark 2.3. Our discussions focus on the low regularity in space, i.e., u(t) ∈ H 2 (Ω), which is sufficient
for the error analysis of the piecewise linear FEM in Section 3. These results cannot be further improved
for u0 ∈ L2 (Ω) or f ∈ Lp (0, T ; L2 (Ω)), due to the limited smoothing properties of the solution operators
(at most of order two in space). For smoother problem data, one may expect higher spatial regularity of
the solution. For example, for the homogeneous problem with a time-independent elliptic operator, there
holds for any β ≥ 0 [29]
ku(t)kḢ 2+β (Ω) ≤ ct−α ku0 kḢ β (Ω) , t > 0.
Naturally, one may expect similar estimates for the case of a time-dependent elliptic operator, provided
both the domain Ω and the coefficient a(x, t) are sufficiently smooth. Further, note that the regularity
analysis extends straightforwardly to the slightly more general elliptic operators with the potential and
convective terms, provided that the coefficients in the lower-order terms have suitable regularity.

3. Semi-discrete Galerkin finite element method

In this part we investigate the semidiscrete Galerkin FEM. Let Th be a shape regular quasi-uniform
triangulation of the domain Ω into simplicial elements, and h be the maximal diameter of the elements.
Let Sh ⊂ H01 (D) be the space of continuous piecewise linear functions over the triangulation Th . Then
we define the L2 (Ω) orthogonal projection Ph : L2 (Ω) → Sh by
(Ph ϕ, χ) = (ϕ, χ) ∀ ϕ ∈ L2 (Ω), ∀ χ ∈ Sh .
The operator Ph satisfies the following error estimate
kPh ϕ − ϕkL2 (Ω) + hk∇(Ph ϕ − ϕ)kL2 (Ω) ≤ chq kϕkH q (Ω) , ϕ ∈ Ḣ q (Ω), q = 1, 2.
The spatially semidiscrete FEM for problem (1.1) reads: find uh (t) ∈ Sh such that
(3.1) (∂tα uh (t), χ) + (a(·, t)∇uh (t), ∇χ) = (f (·, t), χ), ∀χ ∈ Sh , t ∈ (0, T ], with uh (0) = Ph u0 .
 10
Then we define a time-dependent operator Ah (t) : Sh → Sh by
(Ah (t)vh , χ) = (a(·, t)∇vh , ∇χ), ∀ vh , χ ∈ Sh .
Under condition (1.2), Ah (t) : Sh → Sh is bounded and invertible on Sh , and problem (3.1) can be
rewritten as
(3.2) ∂tα uh (t) + Ah (t)uh (t) = Ph f (t), ∀ t ∈ (0, T ], with uh (0) = Ph u0 .
3.1. Perturbation lemmas. In this part we give two crucial perturbation results. We need a time-
dependent Ritz projection operator Rh (t) : H01 (Ω) → Sh defined by
(3.3) (a(·, t)∇Rh (t)ϕ, ∇χ) = (a(·, t)∇ϕ, ∇χ), ∀ϕ ∈ H01 (Ω), χ ∈ Sh .
The operator Rh (t) satisfies the following approximation property [24, p. 99]:
(3.4) kRh (t)ϕ − ϕkL2 (Ω) + hk∇(Rh (t)ϕ − ϕ)kL2 (Ω) ≤ chq kϕkH q (Ω) , ϕ ∈ Ḣ q (Ω), q = 1, 2.
Lemma 3.1. Under conditions (1.2)-(1.3), the following estimate holds:
k(I − Ah (t)−1 Ah (s))vh kL2 (Ω) ≤ c|t − s|kvh kL2 (Ω) , ∀ vh ∈ Sh .
Proof. For any given vh ∈ Sh , let ϕh = Ah (s)vh and wh = Ah (t)−1 ϕh . Then
(Ah (t)wh , χ) = (ϕh , χ) = (Ah (s)vh , χ), ∀ χ ∈ Sh ,
which implies
(a(·, t)∇wh , ∇χ) = (a(·, s)∇vh , ∇χ), ∀ χ ∈ Sh .
Consequently,
(a(·, t)∇(wh − vh ), ∇χ) = ((a(·, s) − a(·, t))∇vh , ∇χ), ∀ χ ∈ Sh .
Let φ ∈ H01 (Ω) be the weak solution of the elliptic problem
(3.5) (a(·, t)∇φ, ∇ξ) = ((a(·, s) − a(·, t))∇vh , ∇ξ), ∀ ξ ∈ H01 (Ω).
By Lax-Milgram theorem, φ satisfies the following a priori estimate:
kφkH 1 (Ω) ≤ ck(a(·, s) − a(·, t))∇vh kL2 (Ω) ≤ c|t − s|kvh kH 1 (Ω) .
Thus, with the Ritz projection Rh (t), cf. (3.3), we have wh − vh = Rh (t)φ.
By the error estimate (3.4) and the inverse inequality,
kwh − vh − φkL2 (Ω) ≤ chkφkH 1 (Ω) ≤ ch|t − s|kvh kH 1 (Ω)
≤ c|t − s|kvh kL2 (Ω) .
Thus, the triangle inequality implies
(3.6) kwh − vh kL2 (Ω) ≤ c|t − s|kvh kL2 (Ω) + kφkL2 (Ω) .
For any ϕ ∈ L2 (Ω), let ξ ∈ Ḣ 2 (Ω) be the solution of the elliptic problem
−∇ · (a(·, t)∇ξ) = ϕ.
Then kξkH 2 (Ω) ≤ ckϕkL2 (Ω) . By substituting ξ into (3.5), we obtain
|(φ, ϕ)| = |(a(·, t)∇φ, ∇ξ)|
= |((a(·, s) − a(·, t))∇vh , ∇ξ)|
= |(vh , ∇ · (a(·, s) − a(·, t))∇ξ)|
≤ c|t − s|kvh kL2 (Ω) kξkH 2 (Ω)
≤ c|t − s|kvh kL2 (Ω) kϕkL2 (Ω) .
This implies (via duality)
|(φ, ϕ)|
kφkL2 (Ω) = sup ≤ c|t − s|kvh kL2 (Ω) .
ϕ∈L2 (Ω) kϕkL2 (Ω)
 11
Substituting the last inequality back into (3.6), we deduce
kwh − vh kL2 (Ω) ≤ c|t − s|kvh kL2 (Ω) .
This completes the proof of Lemma 3.1. 
Remark 3.1. Note that the semidiscrete operator Ah (t) is self-adjoint. Then Lemma 3.1 together with
a duality argument yields
k(I − Ah (s)Ah (t)−1 )vh kL2 (Ω) ≤ c|t − s|kvh kL2 (Ω) , ∀ vh ∈ Sh .
Consequently,
k(Ah (t) − Ah (s))vh kL2 (Ω) ≤ k(I − Ah (s)Ah (t)−1 )Ah (t)vh kL2 (Ω)
≤ c|t − s|kAh (t)vh kL2 (Ω) .
Further, the interpolation between β = 0, 1 yields
kAβh (t)(I − Ah (t)−1 Ah (s))vh kL2 (Ω) ≤ c|t − s|kAβh (t)vh kL2 (Ω) .
The following result is the continuous analogue of Lemma 3.1, and it is independent interest.
Corollary 3.1. Under conditions (1.2)-(1.3), the following estimate holds:
k(I − A(t)−1 A(s))vkL2 (Ω) ≤ c|t − s|kvkL2 (Ω) , ∀ v ∈ H01 (Ω).
Proof. For any v ∈ Ḣ 2 (Ω), there holds Ah (s)Rh (s)v = Ph A(s)v [34, equation (1.34), p. 11]. By the
standard error estimates for Galerkin FEM,
kAh (t)−1 Ah (s)Rh (s)v − A(t)−1 A(s)vkL2 (Ω)
=k(Ah (t)−1 Ph − A(t)−1 )A(s)vkL2 (Ω)
≤ch2 kA(s)vkL2 (Ω) ≤ ch2 kvkH 2 (Ω) .
Then by Lemma 3.1 and the triangle inequality, we deduce
k(I − A(t)−1 A(s))vkL2 (Ω) ≤ k(I − Ah (t)−1 Ah (s))Rh (s)vkL2 (Ω) + ch2 kvkH 2 (Ω)
≤ c|t − s|kRh (s)vkL2 (Ω) + ch2 kvkH 2 (Ω)
≤ c|t − s|kvkL2 (Ω) + c|t − s|kRh (s)v − vkL2 (Ω) + ch2 kvkH 2 (Ω)
≤ c|t − s|kvkL2 (Ω) + (c|t − s| + c)h2 kvkH 2 (Ω) , ∀ v ∈ Ḣ 2 (Ω).
Then the assertion follows by letting h → 0 and noting that the space Ḣ 2 (Ω) is dense in H01 (Ω). 
Lemma 3.2. Under conditions (1.2)-(1.3), the following estimate holds:
(3.7) k(Rh (t) − Rh (s))vkL2 (Ω) ≤ ch2 |t − s|kvkH 2 (Ω) , ∀ v ∈ Ḣ 2 (Ω).
Proof. By the definition of Ritz projection, cf. (3.3), the difference ηh = Rh (t)v − Rh (s)v ∈ Sh satisfies
(a(·, s)∇ηh , ∇χ) = ((a(·, t) − a(·, s))∇(v − Rh (t)v), ∇χ), ∀ χ ∈ Sh .
Let η ∈ H01 (Ω) be the weak solution of the elliptic problem
(3.8) (a(·, s)∇η, ∇ξ) = ((a(·, t) − a(·, s))∇(v − Rh (t)v), ∇ξ), ∀ ξ ∈ H01 (Ω).
By the definition of Rh (s), cf. (3.3), ηh = Rh (s)η and by the error estimate (3.4), there holds
kηh − ηkL2 (Ω) ≤ chkηkH 1 (Ω) ≤ chk(a(·, t) − a(·, s))∇(v − Rh (t)v)kL2 (Ω)
≤ ch2 |t − s|kvkH 2 (Ω) .
The triangle inequality implies
(3.9) kηh kL2 (Ω) ≤ ch2 |t − s|kvkH 2 (Ω) + kηkL2 (Ω) .
Next we use a duality argument to bound kηkL2 (Ω) . For any ϕ ∈ L2 (Ω), let ξ ∈ Ḣ 2 (Ω) be the solution of
the elliptic problem
−∇ · (a(·, s)∇ξ) = ϕ.
 12
Upon substituting ξ into (3.8), we obtain
|(η, ϕ)| = |(a(·, s)∇η, ∇ξ)| = |((a(·, t) − a(·, s))∇(v − Rh (t)v), ∇ξ)|
= |(v − Rh (t)v, ∇ · ((a(·, t) − a(·, s))∇ξ)|
≤ ckv − Rh (t)vkL2 (Ω) |t − s|kξkH 2 (Ω)
≤ ch2 kvkH 2 (Ω) |t − s|kϕkL2 (Ω) ,
which implies (via duality)
kηkL2 (Ω) ≤ ch2 kvkH 2 (Ω) |t − s|.
Substituting the above inequality into (3.9) yields Lemma 3.2. 
3.2. Semidiscrete scheme and error estimates. By the discrete maximal Lp -regularity, one can show
the existence and uniqueness of a FEM solution uh (t). We also have the following stability estimates.
The proof is identical with that for Theorems 2.1–2.3, using the estimates in Section 3.1, and hence it is
omitted.
Theorem 3.1. Let conditions (1.2)-(1.3) be fulfilled, and uh be the solution to problem (3.1).
(i) For u0 ∈ Ḣ β (Ω), 0 ≤ β ≤ 2, and f = 0, then
kAh uh (t)kL2 (Ω) ≤ ct−(1−β/2)α ku0 kḢ β (Ω) and ku0h (t)kL2 (Ω) ≤ ct−(1−αβ/2) ku0 kḢ β (Ω) .
Rt
(ii) For u0 = 0, f ∈ C([0, T ]; L2 (Ω)) and 0 (t − s)α−1 kf 0 (s)kL2 (Ω) ds < ∞, then
kAh uh kLp (0,T ;L2 (Ω)) + k∂tα uh kLp (0,T ;L2 (Ω)) ≤ ckf kLp (0,T ;L2 (Ω)) ,
Z t
ku0h (t)kL2 (Ω) ≤ ct−(1−α) kf (0)kL2 (Ω) + c (t − s)α−1 kf 0 (s)kL2 (Ω) ds,
0
Z t Z t
(t − s)α−1 kAh uh (s)kL2 (Ω) ds ≤ ckf (0)kL2 (Ω) + c (t − s)α−1 kf 0 (s)kL2 (Ω) ds.
0 0
p 2
(iii) For u0 = 0 and f ∈ L (0, T ; L (Ω)) with p ∈ (2/α, ∞), then
kuh (t)kH 1 (Ω) ≤ ckf kLp (0,t;L2 (Ω)) .
Now we derive error estimates for the semidiscrete solution uh . Problem (3.1) can be rewritten as
∂tα uh (t) + Ah (t0 )uh (t) = Ph f (t) + (Ah (t0 ) − Ah (t))uh (t), t ∈ (0, T ], with uh (0) = Ph u0 ,
whose solution is given by
Z t  
(3.10) uh (t) = Fh (t; t0 )Ph u0 + Eh (t − s; t0 ) Ph f (s) + (Ah (t0 ) − Ah (s))uh (s) ds,
0
where the semidiscrete solution operators Fh (t; t0 ) and Eh (t; t0 ) are defined respectively by
Z Z
1 1
Fh (t; t0 ) := ezt z α−1 (z α + Ah (t0 ))−1 dz and Eh (t; t0 ) := ezt (z α + Ah (t0 ))−1 dz.
2πi Γθ,δ 2πi Γθ,δ
Let eh = Ph u − uh . Then by (2.13) and (3.10), eh can be represented by
Z t
eh (t) =(Ph F (t; t0 )u0 − Fh (t; t0 )Ph u0 ) + (Ph E(t − s; t0 ) − Eh (t − s; t0 )Ph )f (s)ds
0
Z t
+ (Ph E(t − s; t0 ) − Eh (t − s; t0 )Ph )(A(t0 ) − A(s))u(s)ds
0
Z t  
+ Eh (t − s; t0 ) Ph (A(t0 ) − A(s))u(s) − (Ah (t0 ) − Ah (s))uh (s) ds
0
4
X
(3.11) =: Ii (t).
i=1
 13
The terms I1 (t) and I2 (t) represent the errors for the homogeneous and inhomogeneous problems with
a time-independent operator A(t0 ), respectively, which have been analyzed: [11, Theorem 3.7] implies
−(1−β/2)α 2
(3.12) kI1 (t0 )kL2 (Ω) ≤ ct0 h ku0 kḢ β (Ω) , β ∈ [0, 2],
and by the argument in [9], there holds (with `h = log(1 + 1/h))
(3.13) kI2 (t0 )kL2 (Ω) ≤ ch2 `2h kf kL∞ (0,T ;L2 (Ω)) .
It remains to bound the two terms I3 (t) and I4 (t), which are given below. We shall discuss the
homogeneous and inhomogeneous problems separately.
Lemma 3.3. Under conditions (1.2) and (1.3), for u0 ∈ L2 (Ω) and f = 0, for the term I3 (t), there holds
kI3 (t0 )kL2 (Ω) ≤ ch2 ku0 kL2 (Ω) .
Proof. By the definitions of the operators E(t; t0 ) and Eh (t; t0 ), we have
Z t0 Z
kI3 (t0 )kL2 (Ω) ≤c |ez(t0 −s) | (z α + A(t0 ))−1 − (z α + Ah (t0 ))−1 Ph L2 (Ω)
k(A(t0 ) − A(s))u(s)k |dz|ds.
0 Γθ,δ

By condition (1.2), for any z ∈ Γθ,δ , we have [6, p. 820]

k(z α + A(t0 ))−1 − (z α + Ah (t0 ))−1 Ph ≤ ch2 ,
where the constant c is independent of z. Meanwhile, condition (1.3) implies
k(A(t0 ) − A(s))u(s)kL2 (Ω) ≤ c|t0 − s|ku(s)kH 2 (Ω) .
Thus by Theorem 2.2,
Z t0 Z t0
2 −1 2
kI3 (t0 )kL2 (Ω) ≤ ch (t0 − s) (t0 − s)ku(s)kH 2 (Ω) ds ≤ ch ku(s)kH 2 (Ω) ds
0 0
Z t0
(3.14) = ch2 s−α ku0 kL2 (Ω) ds ≤ ch2 ku0 kL2 (Ω) .
0
This completes the proof of the lemma. 
Lemma 3.4. Under conditions (1.2) and (1.3), for u0 ∈ L2 (Ω) and f = 0, for the term I4 (t), there holds
Z t0
2
kI4 (t0 )kL2 (Ω) ≤ ch ku0 kL2 (Ω) + c keh (s)kL2 (Ω) ds.
0

Proof. Let eh = Ph u − uh . Using the identity Ph A(s) = Ah (s)Rh (s) [34, (1.34), p. 11] and the triangle
inequality, we derive
Z t0  
kI4 (t0 )kL2 (Ω) = Eh (t0 − s; t0 ) (Ah (t0 )Rh (t0 ) − Ah (s)Rh (s))u(s) − (Ah (t0 ) − Ah (s))uh (s) ds
0 L2 (Ω)
Z t0
≤ Eh (t0 − s; t0 )(Ah (t0 ) − Ah (s))eh (s)ds
0 L2 (Ω)
Z t0  
+ Eh (t0 − s; t0 ) Ah (t0 )(Rh (t0 ) − Ph )u(s) − Ah (s)(Rh (s) − Ph )u(s) ds
0 L2 (Ω)
=: I4,1 (t0 ) + I4,2 (t0 ).
For the term I4,1 (t0 ), by Lemmas 2.2(ii) and 3.1, we have
Z t0
I4,1 (t0 ) = Ah (t0 )Eh (t0 − s; t0 )(I − Ah (t0 )−1 Ah (s))eh (s)ds
0 L2 (Ω)
Z t0
≤ kAh (t0 )Eh (t0 − s; t0 )kk(I − Ah (t0 )−1 Ah (s))eh (s)kL2 (Ω) ds
0
Z t0 Z t0
≤c (t0 − s)−1 (t0 − s)keh (s)kL2 (Ω) ds = c keh (s)kL2 (Ω) ds.
0 0
 14
For the term I4,2 (t0 ), by the triangle inequality, we further split it into
Z t0
I4,2 (t0 ) ≤ Eh (t0 − s; t0 )Ah (t0 )(Rh (t0 ) − Rh (s))u(s)ds
0 L2 (Ω)
Z t0 
+ Eh (t0 − s; t0 )(Ah (t0 ) − Ah (s))(Rh (s) − Ph )u(s) ds =: I04,2 (t0 ) + I004,2 (t0 ).
0 L2 (Ω)

Now by Lemmas 2.2(ii) and 3.2 and Theorem 2.2, we bound I04,2 (t0 ) by
Z t0
I04,2 (t0 ) ≤ kEh (t0 − s; t0 )Ah (t0 )kk(Rh (t0 ) − Rh (s))u(s)kL2 (Ω) ds
0
Z t0
≤c (t0 − s)−1 (t0 − s)h2 ku(s)kH 2 (Ω) ds
0
Z t0
≤ ch2 s−α ku0 kL2 (Ω) ds ≤ ch2 ku0 kL2 (Ω) .
0
00
Likewise, by Lemma 3.1 and Theorem 2.2, we bound I4,2 (t0 ) by
Z t0 
I004,2 (t0 ) = Ah (t0 )Eh (t0 − s; t0 )(I − Ah (t0 )−1 Ah (s))(Rh (s) − Ph )u(s) ds
0 L2 (Ω)
Z t0
≤ kAh (t0 )Eh (t0 − s; t0 )kk(I − Ah (t0 )−1 Ah (s))(Rh (s) − Ph )u(s)kL2 (Ω) ds
0
Z t0
≤c (t0 − s)−1 (t0 − s)k(Rh (s) − Ph )u(s)kL2 (Ω) ds
0
Z t0 Z t0
≤ ch2 ku(s)kH 2 (Ω) ds ≤ ch2 s−α ku0 kL2 (Ω) ds ≤ ch2 ku0 kL2 (Ω) .
0 0
The desired assertion follows by combining the preceding estimates. 
Now we can state the main result of this part, i.e., error estimate on the semidiscrete solution uh .
Theorem 3.2. Under conditions (1.2) and (1.3), for u0 ∈ L2 (Ω) and f = 0, there holds
ku(t) − uh (t)kL2 (Ω) ≤ ch2 t−α ku0 kL2 (Ω) .
Proof. Substituting (3.12) and Lemmas 3.3 and 3.4 into (3.11) yields
Z t0
kPh u(t0 ) − uh (t0 )kL2 (Ω) ≤ ct−α
0 h2
ku k 2
0 L (Ω) + c kPh u(s) − uh (s)kL2 (Ω) ds, ∀ t0 ∈ (0, T ].
0
By Gronwall’s inequality from Lemma 2.3, we obtain
kPh u(t) − uh (t)kL2 (Ω) ≤ ct−α h2 ku0 kL2 (Ω) , ∀ t ∈ (0, T ].
By the approximation property of Ph and Theorem 2.2, we have
ku(t0 ) − Ph u(t0 )kL2 (Ω) ≤ ch2 ku(t0 )kH 2 (Ω) ≤ ct−α 2
0 h ku0 kL2 (Ω) .

The last two estimates together imply the desired result. 

A similar error estimate holds for the inhomogeneous problem.
Theorem 3.3. Under conditions (1.2) and (1.3), for u0 = 0 and f ∈ L∞ (0, T ; L2 (Ω)), there holds
ku(t) − uh (t)kL2 (Ω) ≤ ch2 `2h kf kL∞ (0,t;L2 (Ω)) , with `h = log(1 + 1/h).
Proof. The proof is similar to Theorem 3.2, in view of (3.13), and the following estimates:
kI3 (t0 )kL2 (Ω) ≤ ch2 kf kL∞ (0,t0 ;L2 (Ω)) ,
Z t0
kI4 (t0 )kL2 (Ω) ≤ ch2 kf kL∞ (0,t0 ;L2 (Ω)) + c keh (s)kL2 (Ω) ds,
0
 15
which follow similarly as Lemmas 3.3 and 3.4. Actually, the first follows from (3.14) and Theorem 2.1 by
Z t0
2
kI3 (t0 )kL2 (Ω) ≤ ch ku(s)kH 2 (Ω) ds ≤ ch2 kf kL∞ (0,t0 ;L2 (Ω)) .
0
Similarly, the second follows from the expressions of I04,1 and I004,2 in Lemma 3.4, and Theorem 2.1. 
Remark 3.2. We have only discussed discretization by piecewise linear finite elements. It is of much
interest to extend the analysis to high-order finite elements. This seems missing even for the case of
a time-independent diffusion coefficient when problem data are nonsmooth, partly due to the limited
smoothing property of the solution operators [10].

4. Time discretization
Now we study the time discretization of problem (1.1). We divide the time interval [0, T ] into a uniform
grid, with tn = nτ , n = 0, . . . , N , and τ = T /N being the time step size. Then we approximate the
Riemann-Liouville fractional derivative
d t
Z
1
R α
∂t ϕ(t) = (t − s)−α ϕ(s)ds
Γ(1 − α) dt 0
by the backward Euler (BE) convolution quadrature (with ϕj = ϕ(tj )) [22, 12]:
n
X ∞
X
R α
∂t ϕ(tn ) ≈ τ −α bj ϕn−j := ∂¯τα ϕn , with bj ξ j = (1 − ξ)α .
j=0 j=0

The fully discrete scheme for problem (1.1) reads: find ∈ Sh such that unh
(4.1) ¯α n 0 n
∂ (u − u ) + Ah (tn )u = Ph f (tn ), n = 1, 2, . . . , N,
τ h h h

with the initial condition u0h

= Ph u0 ∈ Sh . Similar to the semidiscrete case, for a given m ∈ N with
1 ≤ m ≤ N , we rewrite (4.1) as
(4.2) ∂¯α (un − u0 ) + Ah (tm )un = Ph f (tn ) + (Ah (tm ) − Ah (tn ))un .
τ h h h h
By means of discrete Laplace transform, the fully discrete solution um
h ∈ Sh is given by
m
X
(4.3) um m 0
h = Fτ,m uh + τ
m−k
Eτ,m [Ph f (tk ) + (Ah (tm ) − Ah (tk ))ukh ],
k=1
nn
where the fully discrete operators Fτ,m
and Eτ,m are respectively defined by (with δτ (ξ) = (1 − ξ)/τ )
Z
1
(4.4) n
Fτ,m = eznτ δτ (e−zτ )α−1 (δτ (e−zτ )α + Ah (tm ))−1 dz,
2πi Γτθ,δ
Z
1
(4.5) n
Eτ,m = eznτ (δτ (e−zτ )α + Ah (tm ))−1 dz,
2πi Γτθ,δ
with the contour Γτθ,δ := {z ∈ Γθ,δ : |=(z)| ≤ π/τ } (oriented with an increasing imaginary part).
The next lemma gives elementary properties of the kernel δτ (e−zτ ).
Lemma 4.1. For any θ ∈ (π/2, π), there exists θ0 ∈ (π/2, π) and positive constants c, c1 , c2 (independent
of τ ) such that for all z ∈ Γτθ,δ
c1 |z| ≤ |δτ (e−zτ )| ≤ c2 |z|, δτ (e−zτ ) ∈ Σθ0 ,
|δτ (e−zτ ) − z| ≤ cτ |z|2 , |δτ (e−zτ )α − z α | ≤ cτ |z|1+α .
By the solution representations (3.10) and (4.3), the temporal error em m
h = uh − uh (tm ) satisfies
m
X m Z tk
X
em m 0 m−k



h = (F h (tm ; t m )P h u0 − F u
τ,m h + τ Eτ,m P h f (t k ) − Eh (tm − s)Ph f (s)ds
k=1 k=1 tk−1
m
X m Z tk
X
m−k
(Ah (tm ) − Ah (tk ))ukh −


+ τ Eτ,m Eh (tm − s)(Ah (tm ) − Ah (s))uh (s)ds
k=1 k=1 tk−1
 3 16
X
(4.6) = Im
i .
i=1

For the first two terms, there hold [12, Theorem 3.5]
−(1−αβ/2)
kIm
1 kL2 (Ω) ≤ cτ tm ku0 kḢ β (Ω) , β ∈ [0, 2],
Z tm
kIm
2 kL2 (Ω) ≤ cτ t−(1−α)
m kf (0)k 2
L (Ω) + cτ (tm − s)α−1 kf 0 (s)kL2 (Ω) ds.
0
To estimate Im
3 ,
n
we need two preliminary bounds on the operator Eτ,m .
m−k
Lemma 4.2. For the operator Eτ,m defined in (4.5), there holds for any β ∈ [0, 1]
Z tk
[τ Aβh (tm )Eτ,m
m−k
− Aβh (tm )Eh (tm − s; tm ) ds] ≤ cτ 2 (tm − tk + τ )−(2−(1−β)α) .
tk−1

Proof. First we consider the case β = 0. By the definition of the operator Eh (t; tm ), we have
Z tk Z Z tk
1
Eh (tm − s; tm ) ds = (z α + Ah (tm ))−1 ez(tm −s) ds dz
tk−1 2πi Γθ,δ tk−1
Z
1
= ez(tm −tk ) z −1 (ezτ − 1)(z α + Ah (tm ))−1 dz.
2πi Γθ,δ
This and the defining relation (4.5) yield
Z tk
m−k
τ Eτ,m − Eh (tm − s; tm ) ds
tk−1
Z
1
ez(tm −tk ) τ (δτ (e−zτ )α + Ah (tm ))−1 − z −1 (ezτ − 1)(z α + Ah (tm ))−1 dz
 
=
2πi Γτθ,δ
Z
1
− ez(tm −tk ) z −1 (ezτ − 1)(z α + Ah (tm ))−1 dz := I + II.
2πi Γθ,δ \Γτθ,δ

For k < m, let δ = (tm − tk + τ )−1 and z = s cos ϕ + is sin ϕ. By Lemma 4.1 and (2.5), we obtain
τ (δτ (e−zτ )α + Ah (tm ))−1 − z −1 (ezτ − 1)(z α + Ah (tm ))−1 ≤ cτ 2 |z|−α+1 , ∀z ∈ Γτθ,δ .
Then the bound on the term I follows by
π
Z τ sin θ
Z θ
kIk ≤ cτ 2 es(tm −tk ) cos θ s−α+1 ds + cτ 2 ecos ϕ δ −α+2 dϕ ≤ cτ 2 (tm − tk + τ )α−2 .
δ −θ

Similarly, Taylor expansion of ezτ , (2.5) and Lemma 4.1 bound the term II by
Z Z ∞
kIIk ≤ cτ |ez(tm −tk ) ||z|−α |dz| ≤ cτ es(tm −tk ) cos θ s−α ds
π
Γθ,δ \Γτθ,δ τ sin θ
Z ∞
s(tm −tk ) cos θ 1−α
≤ cτ 2
e s ds ≤ cτ 2 (tm − tk )−(2−α) .
π
τ sin θ

For k = m, there hold

π
Z τ sin θ
Z θ
kIk ≤ cτ 2 s−α+1 ds + cτ 2 δ −α+2 dϕ ≤ cτ α ,
δ −θ
Z Z ∞
kIIk ≤ c |z|−α−1 |dz| ≤ c s−α−1 ds ≤ cτ α .
π
Γθ,δ \Γτθ,δ τ sin θ

The proof for the case β = 1 is analogous, and the intermediate case β ∈ (0, 1) follows by interpolation. 
n
The next result gives the smoothing property of the operator Eτ,m .
 n 17
Lemma 4.3. For the operator Eτ,m defined in (4.5), there holds
n
kAh (tm )Eτ,m k ≤ c(tn + τ )−1 , n = 0, 1, . . . , N.
Proof. Upon letting δ = (tn + τ )−1 in Γτθ,δ and z = s cos ϕ + is sin ϕ, by (2.5) and Lemma 4.1, we have
Z
1
n
kAh (tm )Eτ,m k= eztn Ah (tm )(δτ (e−zτ )α + Ah (tm ))−1 dz
2πi Γτθ,δ
π
Z τ sin Z θ
θ
≤c estn cos θ ds + c ecos ϕ (tn + τ )−1 dϕ ≤ c(tn + τ )−1 .
(tn +τ )−1 −θ

This completes the proof of the lemma. 

Below we analyze the scheme (4.1) for the homogeneous and inhomogeneous problems separately.
4.1. Error estimate for the homogeneous problem. First we analyze the homogeneous problem. It
suffices to bound the term Im
3 in the splitting (4.6).

Lemma 4.4. Under conditions (1.2)-(1.3), for u0 ∈ L2 (Ω) and f = 0, there holds
m
X
−1
kIm
3 kL2 (Ω) ≤ cτ log(1 + tm /τ )tm ku0 kL2 (Ω) + cτ kekh kL2 (Ω) .
k=1

Proof. Let ekh = ukh− uh (tk ), and Q(t) = (Ah (tm ) − Ah (t))uh (t). Then we split the summand of Im
3 into
Z tk
m−k
τ Eτ,m (Ah (tm ) − Ah (tk ))ukh − Eh (tm − s; tm )Q(s) ds
tk−1
Z tk
m−k m−k k


= τ Eτ,m (Ah (tm ) − Ah (tk ))eh + (τ Eτ,m − Eh (tm − s; tm ) ds)Q(tk )
tk−1
Z tk
+ Eh (tm − s; tm )(Q(tk ) − Q(s)) ds =: Ik + IIk + IIIk .
tk−1

It remains to bound the terms Ik , IIk and IIIk . First, Lemmas 4.3 and 3.1 bound the term kIk kL2 (Ω) by:
m−k
kIk kL2 (Ω) = τ kAh (tm )Eτ,m (I − A−1 k
h (tm )Ah (tk ))eh kL2 (Ω)

≤ cτ (tm − tk + τ )−1 k(I − A−1 k k

h (tm )Ah (tk ))eh kL2 (Ω) ≤ cτ keh kL2 (Ω) .

Second, by Lemma 4.2 (with β = 0) and Remark 3.1, we bound the term IIk by
Z tk
m−k
kIIk kL (Ω) ≤ kτ Eτ,m −
2 Eh (tm − s; tm ) dskkQ(tk )kL2 (Ω)
tk−1

≤ cτ 2 (tm − tk + τ )α−1 kAh (tm )uh (tk )kL2 (Ω) ,

and consequently, by Theorem 3.1(i), we deduce
m
X m
X
kIIk kL2 (Ω) ≤ cτ 2 (tm − tk + τ )α−1 kAh (tm )uh (tk )kL2 (Ω)
k=1 k=1
m
X
≤ cτ 2 ku0 kL2 (Ω) (tm − tk + τ )α−1 t−α
k ≤ cτ ku0 kL2 (Ω) ,
k=1
where the last line follows from the inequality
m
X
τ (tm − tk + τ )α−1 t−α
k ≤ c.
k=1
Last, for the third term IIIk , with k = 1, by Lemma 3.1, we have
Z τ Z τ
kIII1 kL2 (Ω) ≤ kEh (tm − s; tm )Q(τ )kL2 (Ω) ds + kEh (tm − s; tm )Q(s)kL2 (Ω) ds
0 0
 Z τ 18
= kAh (tm )Eh (tm − s; tm )Ah (tm )−1 Q(τ )kL2 (Ω) ds
0
Z τ
+ kAh (tm )Eh (tm − s; tm )Ah (tm )−1 Q(s)kL2 (Ω) ds
0
Z τ
≤c (tm − s)−1 ((tm − τ ) + (tm − s)) dsku0 kL2 (Ω) ≤ cτ ku0 kL2 (Ω) .
0

Meanwhile, for k > 1, we further split the term IIIk into

Z tk Z s
IIIk = Eh (tm − s; tm ) Q0 (ξ) dξ ds
tk−1 tk
Z tk Z s
= Eh (tm − s; tm ) (Ah (tm ) − Ah (ξ))u0h (ξ) dξ ds
tk−1 tk
Z tk Z s
− Eh (tm − s; tm ) A0h (ξ)uh (ξ) dξ ds =: IIIk,1 + IIIk,2 .
tk−1 tk

By Lemmas 3.1 and 2.2(ii), the term IIIk,1 for any k > 1 can be bounded by
Z tk Z tk
kIIIk,1 kL2 (Ω) ≤ kAh (tm )Eh (tm − s; tm )k k(I − Ah (tm )−1 Ah (ξ))u0h (ξ)kL2 (Ω) dξ ds
tk−1 s
Z tk Z tk
≤c (tm − s)−1 (tm − ξ)ku0h (ξ)kL2 (Ω) dξ ds
tk−1 s
Z tk Z tk
≤c (tm − s)−1 (tm − ξ)ξ −1 ku0 kL2 (Ω) dξ ds,
tk−1 s

where the last step is due to Theorem 3.1(i). Now we note the elementary inequality
Z tk Z tk Z tk Z ξ
(tm − s)−1 (tm − ξ)ξ −1 dξds = (tm − ξ)ξ −1 (tm − s)−1 dsdξ
tk−1 s tk−1 tk−1
Z tk Z tk
−1 −1
≤ (tm − ξ)(tm − ξ) τξ dξ = τ ξ −1 dξ.
tk−1 tk−1

Consequently,
Z tk
kIIIk,1 kL2 (Ω) ≤ cτ ξ −1 dξku0 kL2 (Ω) ,
tk−1
and
m
X Z tm
kIIIk,1 kL2 (Ω) ≤ cτ ξ −1 dξku0 kL2 (Ω) = cτ log(tm /τ ).
k=2 τ

Similarly, by Theorem 3.1(i), the term IIIk,2 for any k > 1 is bounded by
Z tk Z tk
kIIIk,2 kL2 (Ω) = kEh (tm − s; tm )k kA0h (ξ)uh (ξ)kL2 (Ω) dξ ds
tk−1 s
Z tk Z tk
≤c (tm − s)α−1 ξ −α dξ dsku0 kL2 (Ω)
tk−1 s
Z tk
≤ cτ (tm − s)α−1 s−α dsku0 kL2 (Ω) ,
tk−1
Rt Rt
where the last line follows from the inequality s k ξ −α dξ ≤ s−α s k dξ ≤ s−α τ. Thus,
Xm Z tm
kIIIk,2 kL2 (Ω) ≤ cτ (tm − s)α−1 s−α dsku0 kL2 (Ω) ≤ cku0 kL2 (Ω) .
k=1 0
 19
Hence, there holds
m
X
kIIIk kL2 (Ω) ≤ cτ (1 + log(tm /τ ))ku0 kL2 (Ω) .
k=1

Combining the preceding estimates completes the proof of the lemma. 

Now we can state an error estimate for the homogeneous problem.

Theorem 4.1. Under conditions (1.2)-(1.3), u0 ∈ L2 (Ω) and f = 0, there holds
−1
kum
h − uh (tm )kL2 (Ω) ≤ cτ tm log(1 + tm /τ ))ku0 kL2 (Ω) .

Proof. It follows from Lemma 4.4 that

m
X
−1
kem
h kL2 (Ω) ≤ cτ (tm + log(1 + tm /τ ))ku0 kL2 (Ω) + cτ kekh kL2 (Ω) .
k=1

The desired estimate follows from a variant of the discrete Gronwall’s inequality [34, p. 258]. 

Remark 4.1. The logarithmic factor log(1 + tm /τ ) is also present for the BE method for standard
parabolic problems with a time-dependent diffusion coefficient [24, Theorem 2, p. 95]. For u0 ∈ Ḣ β (Ω),
β ∈ (0, 2], it may be improved:
−(1−βα/2)
kum
h − uh (tm )kL2 (Ω) ≤ cτ tm ku0 kḢ β (Ω) .
In fact, the argument of Lemma 4.4 (together with Theorem 3.1(i)) implies
Z tk Z tk
kIIIk,1 kL2 (Ω) ≤ c (tm − s)−1 (tm − ξ)ku0h (ξ)kL2 (Ω) dξ ds
tk−1 s
Z tk Z tk
≤c (tm − s)−1 (tm − ξ)ξ −(1−αβ/2) ku0 kḢ β (Ω) dξ ds
tk−1 s
Z tk
≤ cτ ξ −(1−βα/2) dξku0 kḢ β (Ω) ,
tk−1

and
Z tk Z tk
kIIIk,2 kL2 (Ω) = kEh (tm − s; tm )k kA0h (ξ)uh (ξ)kL2 (Ω) dξ ds
tk−1 s
Z tk Z tk
≤c (tm − s)α−1 ξ −(1−β/2)α dξ dsku0 kḢ β (Ω)
tk−1 s
Z tk
≤ cτ (tm − s)α−1 s−(1−β/2)α dsku0 kḢ β (Ω) .
tk−1

Thus, for any β ∈ (0, 2], there holds

Xm Z tk
E(tm − s; tm )(Q(tk ) − Q(s)) ds ≤ cτ ku0 kḢ β (Ω) .
k=1 tk−1

Then the desired estimate follows by repeating the argument for Theorem 4.1.
4.2. Error estimate for the inhomogeneous problem. Now we give the temporal discretization
error for the inhomogeneous problem.
Rt
Theorem 4.2. Under conditions (1.2)-(1.3), u0 = 0, f ∈ C([0, T ]; L2 (Ω)) and 0 (t−s)α−1 kf 0 (s)kL2 (Ω) ds <
∞ for any 0 < t ≤ T , there holds
 Z tm 
−(1−α)
kumh − u (t )k 2
h m L (Ω) ≤ cτ tm kf (0)k 2
L (Ω) + (tm − s)α−1 kf 0 (s)kL2 (Ω) ds .
0
 20
Proof. The argument is similar to Theorem 4.1, and thus we only sketch the main steps. It suffices to
bound the terms IIk and IIIk in Lemma 4.4. By Theorem 3.1(iii), Remark 3.1 and Lemma 4.2 (with
β = 1/2), the following estimate holds:
m m
1/2
X X
kIIk kL2 (Ω) ≤ cτ 2 (tm − tk + τ )α/2−1 kAh (tm )uh (tk )kL2 (Ω)
k=1 k=1
Xm
≤ cτ 2 (tm − tk + τ )α/2−1 kf kL∞ (0,tm ;L2 (Ω)) ≤ cτ kf kL∞ (0,tm ;L2 (Ω)) .
k=1
Meanwhile, upon noting tm ≤ T , we have
Z tm
kf kL∞ (0,tm ;L2 (Ω)) ≤ kf (0)kL2 (Ω) + kf 0 (s)kL2 (Ω) ds
0
Z tm
≤ kf (0)kL2 (Ω) + cT (tm − s)α−1 kf 0 (s)kL2 (Ω) ds.
0
Next, the two terms IIIk,1 and IIIk,2 can be bounded respectively by
Z tk Z tk
−1
kIIIk,1 kL2 (Ω) ≤ c (tm − s) (tm − ξ)ku0h (ξ)kL2 (Ω) dξ ds
tk−1 s
Z tk Z ξ
=c (tm − ξ)ku0h (ξ)kL2 (Ω) (tm − s)−1 ds dξ
tk−1 tk−1
Z tk
≤ cτ ku0h (ξ)kL2 (Ω) dξ,
tk−1

and
Z tk Z tk
kIIIk,2 kL2 (Ω) ≤ kEh (tm − s; tm )k kA0h (ξ)uh (ξ)kL2 (Ω) dξ ds
tk−1 s
Z tk Z tk
≤c (tm − s) α−1
kA0h (ξ)uh (ξ)kL2 (Ω) dξ ds
tk−1 s
Z tk Z ξ
=c kA0h (ξ)uh (ξ)kL2 (Ω) (tm − s)α−1 dsdξ
tk−1 tk−1
Z tk
≤ cτ (tm − ξ)α−1 kA0h (ξ)uh (ξ)kL2 (Ω) dξ.
tk−1

Then by Theorem 3.1(ii), we have

m
X Z tm
kIIIk,1 kL2 (Ω) ≤ cτ ku0h (s)kL2 (Ω) ds
k=1 0
Z tm Z tm Z s
≤ cτ kf (0)kL2 (Ω) ds +(s − ξ)α−1 kf 0 (ξ)kL2 (Ω) dξds
0 0 0
 Z tm Z tm Z tm 
= cτ kf (0)kL2 (Ω) ds + kf 0 (ξ)kL2 (Ω) (s − ξ)α−1 dsdξ
0 0 ξ
 Z tm 
≤ cτ kf (0)kL2 (Ω) + (tm − ξ)α−1 kf 0 (ξ)kL2 (Ω) dξ ,
0
and similarly, by Theorem 3.1(ii), we deduce
Xm Z tm
kIIIk,2 kL2 (Ω) ≤ cτ (tm − s)α−1 kA0h (s)uh (s)kL2 (Ω) ds
k=1 0
 Z tm 
≤ cτ kf (0)kL2 (Ω) + (tm − s)α−1 kf 0 (s)kL2 (Ω) ds .
0
 21

These estimates together with discrete Gronwall’s inequality complete the proof. 
Remark 4.2. The proof techniques in this section apply to other first-order methods, e.g., L1 scheme
[21], and similar error estimates can also be derived for these methods.
Remark 4.3. We briefly comment on the dependence of the constant c in error estimates on the fractional
order α. At a few occasions, it can blow up as α → 1− ; see e.g., I3 (t0 ) in Lemma 3.3, I4,2 (t0 ) in Lemma
3.4 and IIIk,2 in Lemma 4.4. This phenomenon does not fully agree with the results for the continuous
model. Such a blowup phenomenon appears also in some existing error analysis; see, e.g., [28, eq. (2.2)]
and [32, Lemma 4.3], and it is of interest to further refine the estimates to fill in the gap.

5. Numerical results
Now we present numerical examples to verify the theoretical results in Sections 3 and 4. We consider
problem (1.1) with a time-dependent elliptic operator A(t) = −(2 + cos(t))∆ on the domain Ω = (0, 1)
and the following two sets of problem data:
(a) u0 (x) = x−1/4 ∈ H 1/4− (Ω) with  ∈ (0, 1/4) and f ≡ 0.
(b) u0 (x) = 0 and f = et (1 + χ(0,1/2) (x)).
Unless otherwise specified, the final time T is fixed at T = 1.
We divide the domain Ω into M subintervals of equal length h = 1/M . The numerical solutions are
computed using the Galerkin FEM in space, and the backward Euler (BE) CQ or L1 scheme in time. To
evaluate the convergence, we compute the spatial error es and temporal error et , respectively, defined by
es (tN ) = kuh (tN ) − u(tN )kL2 (Ω) and et (tN ) = kuN
h − uh (tN )kL2 (Ω) .
Since the exact solution is unavailable, we compute reference solutions on a finer mesh: for the error es ,
we take the time step τ = 1/10000 and mesh size h = 1/1280, and for the error et , take h = 1/100 and
τ = 1/10000, unless otherwise specified.
First we examine the spatial convergence of the semidiscrete Galerkin scheme (3.2). The spatial errors
for case (a) are shown in Table 1, which indicates a steady O(h2 ) rate for the semidiscrete scheme (3.2),
just as predicted by Theorem 3.2. The O(h2 ) rate holds for all three fractional orders and different
terminal times. Since the initial data is nonsmooth, the spatial error es (tN ) decreases with the time tN ,
which is in good agreement with the regularity result in Theorem 2.2. To further illuminate the precise
dependence of the spatial error es (tN ) on tN , in Table 2, we present the error es as the time tN → 0
−(2−β)α/2 2
for case (a). By repeating the argument for Theorem 3.2, there holds es (tN ) ≤ ctN h kvkḢ β (Ω) ,
0 ≤ β ≤ 2. For case (a), this estimate predicts an exponent 7α/8 for the dependence on the time tN ,
which gives the numbers shown in the bracket in Table 2. Table 2 indicates that the empirical rate agrees
excellently with the predicted one, fully confirming the analysis. Similar observations hold also for the
numerical results for the inhomogeneous problem in case (b), cf. Table 3. These results fully support the
error analysis of the semidiscrete scheme in Section 3.

Table 1. Spatial errors es for example (a) with τ = 1/10000 and h = 1/M .
N
T 10 20 40 80 160 rate
α
0.25 1.44e-5 3.62e-6 9.06e-7 2.26e-7 5.62e-8 2.00 (2.00)
1 0.50 1.02e-5 2.56e-6 6.40e-7 1.60e-7 3.97e-8 2.00 (2.00)
0.75 5.18e-6 1.30e-6 3.25e-7 8.12e-8 2.02e-8 2.00 (2.00)
0.25 6.26e-5 1.57e-5 3.93e-6 9.80e-7 2.44e-7 2.01 (2.00)
10−3 0.50 2.12e-4 5.31e-5 1.33e-5 3.32e-6 8.24e-7 2.01 (2.00)
0.75 5.99e-4 1.50e-4 3.75e-5 9.36e-6 2.33e-6 2.01 (2.00)

Next we turn to the temporal convergence, and present numerical results for both BE and L1 schemes,
cf. Remark 4.2. The temporal errors et for case (a) at two time instances are given in Table 4, which
 Table 2. Spatial errors es for example (a) with h = 1/200 and N = 10000, at T = 10−k . 22

k
2 3 4 5 6 7 rate
α
0.25 2.40e-6 4.04e-6 6.58e-6 1.05e-5 1.65e-5 2.66e-5 0.21 (0.22)
0.5 5.28e-6 1.31e-5 3.36e-5 9.02e-5 2.40e-4 6.40e-4 0.43 (0.44)
0.75 9.95e-6 3.90e-5 1.70e-4 7.45e-4 3.26e-3 1.39e-2 0.64 (0.66)

Table 3. Spatial errors es for example (b) at T = 1 with τ = 1/10000 and h = 1/M .
M
10 20 40 80 160 rate
α
0.25 2.03e-4 5.06e-5 1.27e-5 3.16e-6 7.85e-7 2.01 (2.00)
0.50 2.08e-4 5.19e-5 1.30e-5 3.24e-6 8.04e-7 2.01 (2.00)
0.75 2.13e-4 5.32e-5 1.33e-5 3.32e-6 8.25e-7 2.01 (2.00)

indicate an O(τ ) convergence rate for both time stepping schemes. Further, the accuracy of both schemes
is largely comparable. The convergence is very steady for both schemes, and the convergence rate is
independent of the fractional order α and the final time tN (so long as it is fixed). Further, it is observed
that the error et decreases with the time tN . To show the dependence of the temporal error et (tN ) with
the time tN , in Table 5, we present et (tN ) as the time tN tends to zero. In view of Remark 4.1, there
−(1−βα/2)
holds et (tN ) ≤ cτ tN ku0 kḢ β (Ω) , 0 < β ≤ 2. This estimate predicts a decay O(N −α/8 ) for case
(a), which agrees excellently with the empirical rate (in the bracket) in Table 5, thereby confirming the
sharpness of the error estimate. These observations hold also for the inhomogeneous problem in case
(b), cf. Table 6. These numerical results fully support the error analysis of the fully discrete scheme in
Section 4.
Table 4. Temporal errors et for example (a) with h = 1/100 and τ = T /N .
N
T method 100 200 400 800 1600 rate
α
0.25 5.43e-5 2.71e-5 1.35e-5 6.76e-6 3.38e-6 1.00 (1.00)
BE 0.50 9.49e-5 4.73e-5 2.36e-5 1.18e-5 5.90e-6 1.00 (1.00)
1 0.75 9.01e-5 4.49e-5 2.24e-5 1.12e-5 5.59e-6 1.00 (1.00)
0.25 4.35e-5 2.17e-5 1.08e-5 5.41e-6 2.70e-6 1.00 (1.00)
L1 0.50 6.33e-5 3.15e-5 1.57e-5 7.84e-6 3.92e-6 1.00 (1.00)
0.75 5.12e-5 2.54e-5 1.26e-5 6.29e-6 3.14e-6 1.01 (1.00)
0.25 2.00e-4 9.99e-5 4.99e-5 2.49e-5 1.25e-5 1.00 (1.00)
BE 0.50 8.16e-4 4.08e-4 2.04e-4 1.02e-4 5.10e-5 1.00 (1.00)
10−3 0.75 7.58e-4 3.79e-4 1.89e-4 9.46e-5 4.73e-5 1.00 (1.00)
0.25 1.69e-4 8.43e-5 4.21e-5 2.10e-5 1.05e-5 1.00 (1.00)
L1 0.50 8.08e-4 3.99e-4 1.98e-4 9.84e-5 4.90e-5 1.01 (1.00)
0.75 8.28e-4 4.11e-4 2.04e-4 1.02e-4 5.07e-5 1.01 (1.00)

Table 5. Temporal errors et for example (a) with α = 0.5, h = 10−3 and N = 5, at T = 10−k .
k
3 4 5 6 7 8 rate
α
0.5 BE 1.65e-2 1.06e-2 8.79e-3 7.45e-3 6.33e-3 5.39e-3 0.07 (0.06)
L1 1.91e-2 1.41e-2 1.06e-2 8.95e-3 7.59e-3 6.46e-3 0.07 (0.06)
0.8 BE 1.61e-2 1.23e-2 9.52e-3 7.44e-3 5.84e-3 4.56e-3 0.11 (0.10)
L1 1.83e-2 1.40e-2 1.08e-2 8.46e-3 6.64e-3 5.19e-3 0.11 (0.10)
 Table 6. Temporal errors et for example (b) at T = 1 with h = 1/100 and τ = T /N . 23

N
method 100 200 400 800 1600 rate
α
0.25 3.26e-6 1.63e-6 8.15e-7 4.07e-7 2.04e-7 1.00 (1.00)
BE 0.50 4.76e-6 2.37e-6 1.18e-6 5.92e-7 2.96e-7 1.00 (1.00)
0.75 2.76e-6 1.37e-6 6.84e-7 3.41e-7 1.71e-7 1.00 (1.00)
0.25 2.33e-6 1.17e-6 5.85e-7 2.93e-7 1.47e-7 1.00 (1.00)
L1 0.50 3.25e-6 1.64e-6 8.29e-7 4.18e-7 2.10e-7 0.99 (1.00)
0.75 1.85e-6 9.89e-7 5.22e-7 2.72e-7 1.41e-7 0.94 (1.00)

Acknowledgements
The research of B. Li is partially supported by a Hong Kong RGC grant (Project No. 15300817), and
that of Z. Zhou by a start-up grant from the Hong Kong Polytechnic University and Hong Kong RGC
grant No. 25300818.

References
[1] G. Akrivis, B. Li, and C. Lubich. Combining maximal regularity and energy estimates for time discretizations of
quasilinear parabolic equations. Math. Comp., 86:1527–1552, 2017.
[2] A. A. Alikhanov. A new difference scheme for the time fractional diffusion equation. J. Comput. Phys., 280:424–438,
2015.
[3] E. G. Bajlekova. Fractional Evolution Equations in Banach Spaces. PhD thesis, Eindhoven University of Technology,
Eindhoven, Netherlands, 2001.
[4] S. D. Eidelman and A. N. Kochubei. Cauchy problem for fractional diffusion equations. J. Diff. Eq., 199(2):211–255,
2004.
[5] C. M. Elliott and S. Larsson. Error estimates with smooth and nonsmooth data for a finite element method for the
Cahn-Hilliard equation. Math. Comp., 58(198):603–630, S33–S36, 1992.
[6] H. Fujita and T. Suzuki. Evolution problems. In Handbook of Numerical Analysis, Vol. II, pages 789–928. North-
Holland, Amsterdam, 1991.
[7] G.-H. Gao, Z.-Z. Sun, and H.-W. Zhang. A new fractional numerical differentiation formula to approximate the Caputo
fractional derivative and its applications. J. Comput. Phys., 259:33–50, 2014.
[8] D. Henry. Geometric Theory of Semilinear Parabolic Equations. Springer-Verlag, Berlin-New York, 1981.
[9] B. Jin, R. Lazarov, J. Pasciak, and Z. Zhou. Error analysis of semidiscrete finite element methods for inhomogeneous
time-fractional diffusion. IMA J. Numer. Anal., 35(2):561–582, 2015.
[10] B. Jin, R. Lazarov, and Z. Zhou. Numerical methods for time-fractional evolution equations with nonsmooth data: a
concise overview. Preprint, arXiv:1805.11309.
[11] B. Jin, R. Lazarov, and Z. Zhou. Error estimates for a semidiscrete finite element method for fractional order parabolic
equations. SIAM J. Numer. Anal., 51(1):445–466, 2013.
[12] B. Jin, R. Lazarov, and Z. Zhou. Two fully discrete schemes for fractional diffusion and diffusion-wave equations with
nonsmooth data. SIAM J. Sci. Comput., 38(1):A146–A170, 2016.
[13] B. Jin, B. Li, and Z. Zhou. Numerical analysis of nonlinear subdiffusion equations. SIAM J. Numer. Anal., 56(1):1–23,
2018.
[14] S. Karaa. Semidiscrete finite element analysis of time fractional parabolic problems: a unified approach. SIAM J.
Numer. Anal., 56(3):1673–1692, 2018.
[15] S. Karaa, K. Mustapha, and A. K. Pani. Finite volume element method for two-dimensional fractional subdiffusion
problems. IMA J. Numer. Anal., 37(2):945–964, 2017.
[16] A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo. Theory and Applications of Fractional Differential Equations.
Elsevier Science B.V., Amsterdam, 2006.
[17] A. Kubica and M. Yamamoto. Initial-boundary value problems for fractional diffusion equations with time-dependent
coeffcients. Fract. Calc. Appl. Anal., 21(2):276–311, 2018.
[18] P. C. Kunstmann, B. Li, and C. Lubich. Runge-Kutta time discretization of nonlinear parabolic equations studied via
discrete maximal parabolic regularity. Found. Comput. Math., 2017, DOI: 10.1007/s10208-017-9364-x.
[19] H. Lee, J. Lee, and D. Sheen. Laplace transform method for parabolic problems with time-dependent coefficients. SIAM
J. Numer. Anal., 51(1):112–125, 2013.
[20] B. Li and W. Sun. Regularity of the diffusion-dispersion tensor and error analysis of FEMs for a porous media flow.
SIAM J. Numer. Anal., 53:1418–1437, 2015.
[21] Y. Lin and C. Xu. Finite difference/spectral approximations for the time-fractional diffusion equation. J. Comput.
Phys., 225(2):1533–1552, 2007.
[22] C. Lubich. Discretized fractional calculus. SIAM J. Math. Anal., 17(3):704–719, 1986.
 24
[23] Y. Luchko. Maximum principle for the generalized time-fractional diffusion equation. J. Math. Anal. Appl., 351(1):218–
223, 2009.
[24] M. Luskin and R. Rannacher. On the smoothing property of the Galerkin method for parabolic equations. SIAM J.
Numer. Anal., 19(1):93–113, 1982.
[25] W. McLean. Regularity of solutions to a time-fractional diffusion equation. ANZIAM J., 52(2):123–138, 2010.
[26] W. McLean and K. Mustapha. Convergence analysis of a discontinuous Galerkin method for a sub-diffusion equation.
Numer. Algorithms, 52(1):69–88, 2009.
[27] W. McLean and K. Mustapha. Time-stepping error bounds for fractional diffusion problems with non-smooth initial
data. J. Comput. Phys., 293:201–217, 2015.
[28] K. Mustapha. FEM for time-fractional diffusion equations, novel optimal error analyses. Math. Comp., 87(313):2259–
2272, 2018.
[29] K. Sakamoto and M. Yamamoto. Initial value/boundary value problems for fractional diffusion-wave equations and
applications to some inverse problems. J. Math. Anal. Appl., 382(1):426–447, 2011.
[30] P. Sammon. Fully discrete approximation methods for parabolic problems with nonsmooth initial data. SIAM J.
Numer. Anal., 20(3):437–470, 1983.
[31] G. Savaré. A(Θ)-stable approximations of abstract Cauchy problems. Numer. Math., 65(3):319–335, 1993.
[32] M. Stynes, E. O’Riordan, and J. L. Gracia. Error analysis of a finite difference method on graded meshes for a time-
fractional diffusion equation. SIAM J. Numer. Anal., 55(2):1057–1079, 2017.
[33] Z.-Z. Sun and X. Wu. A fully discrete difference scheme for a diffusion-wave system. Appl. Numer. Math., 56(2):193–209,
2006.
[34] V. Thomée. Galerkin Finite Element Methods for Parabolic Problems. Springer-Verlag, Berlin, second edition, 2006.
[35] Y. Yan, M. Khan, and N. J. Ford. An analysis of the modified L1 scheme for time-fractional partial differential equations
with nonsmooth data. SIAM J. Numer. Anal., 56(1):210–227, 2018.
[36] S. B. Yuste and L. Acedo. An explicit finite difference method and a new von Neumann-type stability analysis for
fractional diffusion equations. SIAM J. Numer. Anal., 42(5):1862–1874, 2005.
[37] R. Zacher. A De Giorgi–Nash type theorem for time fractional diffusion equations. Math. Ann., 356(1):99–146, 2013.

Department of Computer Science, University College London, Gower Street, London, WC1E 2BT, UK.
E-mail address: b.jin@ucl.ac.uk,bangti.jin@gmail.com

Department of Applied Mathematics, The Hong Kong Polytechnic University, Kowloon, Hong Kong
E-mail address: bygli@polyu.edu.hk

Department of Applied Mathematics, The Hong Kong Polytechnic University, Kowloon, Hong Kong
E-mail address: zhizhou@polyu.edu.hk

View publication stats