# A mixed FEM for a time-fractional Fokker-Planck model††thanks: The work of
Naveed Ahmed has been funded through the GUST Internal Seed Grant, case number
8. The work of Samir Karaa was supported by Sultan Qaboos University under
Grant IG/SCI/MATH/21/01.

Samir Karaa, Kassem Mustapha, Naveed Ahmed Department of Mathematics, Sultan
Qaboos University, Al-Khod 123, Muscat, Oman, Email<EMAIL_ADDRESS>of Mathematics and Statistics, University of New South Wales, Sydney,
Australia, Email<EMAIL_ADDRESS>for Applied
Mathematics and Bioinformatics (CAMB), Department of Mathematics and Natural
Sciences, Gulf University for Science and Technology, Hawally, Kuwait,
<EMAIL_ADDRESS>
###### Abstract

We propose and analyze a mixed finite element method for the spatial
approximation of a time-fractional Fokker–Planck equation in a convex
polyhedral domain, where the given driving force is a function of space.
Taking into account the limited smoothing properties of the model, and
considering an appropriate splitting of the errors, we employed a sequence of
clever energy arguments to show optimal convergence rates with respect to both
approximation properties and regularity results. In particular, error bounds
for both primary and secondary variables are derived in $L^{2}$-norm for cases
with smooth and nonsmooth initial data. We further investigate a fully
implicit time-stepping scheme based on a convolution quadrature in time
generated by the backward Euler method. Our main result provides pointwise-in-
time optimal $L^{2}$-error estimates for the primary variable. Numerical
examples are then presented to illustrate the theoretical contributions.

Keywords. fractional Fokker–Planck equation, mixed finite element method,
convolution quadrature, error analysis, smooth and nonsmooth initial data.
AMS Subject Classification. 65M60, 65M12, 65M15, 65M70 and 35S10

##  1 Introduction

Consider a bounded convex polygonal domain $\Omega$ in $\mathbb{R}^{2}$ with a
boundary $\partial\Omega$, and let $T>0$ be a fixed value. The primary
objective of this paper is to develop and analyze a mixed finite element
method (FEM) for solving the time-fractional Fokker-Planck model problem [27,
10]: with $0<\alpha<1,$

$\partial_{t}u-\nabla\cdot\bigl{(}^{R}\partial_{t}^{1-\alpha}\kappa\nabla
u-{\bf F}\,^{R}\partial_{t}^{1-\alpha}u\bigr{)}=0,\quad\mbox{ in
}\Omega\times(0,T],$ (1)

subject to homogeneous Dirichlet boundary conditions, with $u(x,0)=u_{0}(x)$
at the initial time $t=0.$ Here $\partial_{t}=\partial/\partial t$, the time-
independent diffusivity $\kappa\in L^{\infty}(\Omega)$ and bounded below by
$c_{0}>0$. In (1), the Riemann-Liouville fractional derivative
${}^{R}\partial_{t}^{1-\alpha}=\frac{d}{dt}\mathcal{I}^{\alpha}$, where the
time-fractional integral

$\mathcal{I}^{\alpha}{\varphi}(t):=\int_{0}^{t}\omega_{\alpha}(t-s)\varphi(s)\,ds\quad\text{with}\quad\omega_{\alpha}(t):=\frac{t^{\alpha-1}}{\Gamma(\alpha)}.$

We assume that the time-independent driving forces ${\bf F}$ and
$\nabla\cdot{\bf F}$ are bounded in $\Omega$. Since ${\bf F}$ is independent
of $t$, applying $\mathcal{I}^{1-\alpha}$ to (1) and using
$\mathcal{I}^{1-\alpha}(^{R}\partial_{t}^{1-\alpha}u)=u$ yield the alternative
formulation

$\partial_{t}^{\alpha}u+{\cal L}u=0\quad\mbox{ on }\Omega\times(0,T],$ (2)

where the second-order non-selfadjoint operator ${\cal
L}u=-\nabla\cdot(\kappa\nabla u)+\nabla\cdot({\bf F}u)$. In (2),
$\partial_{t}^{\alpha}u(t):=\mathcal{I}^{1-\alpha}{\partial_{t}u}(t)$ is the
time-fractional Caputo derivative of order $\alpha$ ($0<\alpha<1$). As
$\alpha$ approaches $1^{-}$, $\partial_{t}^{\alpha}u$ converges to
$\partial_{t}u$, and consequently, problem (1) (equivalently (2)) reduces to
the standard Fokker-Planck equation. For a comprehensive analysis of the well-
posedness of the model problem (1) with a general driving force ${\bf F}$, we
refer the reader to [25].

Various authors have studied the numerical solution of (2), mostly for a 1D
spatial domain and a smooth continuous solution (which is not practically the
case). For instance, Deng [6] transferred (2) to a system of fractional
ordinary differential equations (ODEs) by discretizing the spatial
derivatives, and then applied a predictor-corrector approach combined with the
method of lines. Cao et al. [3] adopted a similar approach and solved the
resulting system using a second-order backward Euler scheme. Chen et al. [4]
investigated three different implicit finite difference techniques, in each of
which the diffusion term was approximated by the standard second-order
difference approximation at the advanced time level. In contrast, Jiang [12]
established monotonicity properties of the numerical solutions obtained by
using these schemes, and so showed that the time-stepping preserves the non-
negativity of the solution. Jiang and Xu [13] proposed a finite volume method,
and Yang et al. [34] a spectral collocation method. Fairweather et al. [8]
developed an orthogonal spline collocation method in space combined with the
backward Euler method in time. Saadmandi [31] studied a collocation method
based on shifted Legendre polynomials in time and Sinc functions in space.
Vong and Wang [33] have analyzed a high-order, compact difference scheme. Cui
[5] applied a high-order approximation for the time-fractional derivative
combined with a compact exponential finite difference scheme.

For problem (1) with a time-space dependent ${\bf F}$, first-order backward
Euler and second-order corrected Crank-Nicolson time-stepping schemes were
studied in [19, 11] and [29], respectively. The stability and error analyses
were carried out under certain regularity assumptions on the given data
excluding the case of nonsmooth initial data, this is in addition to some
restrictions on the range of the fractional exponent $\alpha.$ For the spatial
discretization, a continuous, piecewise-linear Galerkin FEM was considered in
a series of papers, see the recent papers [20, 26] and related references
therein. In [26], based on novel energy arguments, $\alpha$-robust stability,
and optimal convergence results were investigated for both cases of smooth and
nonsmooth initial data.

In the special case ${\bf F}\equiv 0$ in (1), mixed FEM were studied only in
[35, 15, 16, 17]. In [35], the authors established mixed finite element error
analysis under higher (unpractical) regularity assumption on the continuous
solution $u$. Subsequently, Karaa [15, 16] has exploited the inverse of two
mixed operators and derived optimal error estimates allowing nonsmooth initial
data. The studies included both semidiscrete and fully discrete schemes. Later
on, an alternative approach was adopted in [17] to handle the case of the
time-dependent diffusivity coefficient. The error analyses in these papers are
not applicable for the case of non-zero ${\bf F}$, and extending it to cover
this case is not feasible.

The main aim of this work is: (1) develop a spatial mixed FEM for problem (2)
and derive $\alpha$-robust (that is, remains valid as $\alpha$ approaches $1$)
error estimates which are not only optimal in approximation property but also
in regularity of the solution; (2) combine our mixed FEM with a convolution
quadrature generated by the backward Euler method in time and show optimal
error estimates for the primary variable. Our convergence analysis approach
involves splitting the errors into two distinct components. Then, a suitable
exploitation of the results [15] and [16], and the use of a sequence of
challenging energy arguments in combination with general weakly singular
Gronwall inequalities (see [7, Theorem 3.1] and [7, Theorem 6.1] for the
continuous version and discrete version, respectively) allow us to prove the
required optimal error bounds for both smooth and nonsmooth initial data.

The paper is organized as follows. In Section 2, we introduce the necessary
notations and gather some technical estimates involving fractional integrals.
In Section 3, we define the semidiscrete mixed finite element scheme. Error
estimates follow in Theorem 5 showing optimal convergence rates for both
primary and flux variables. We combine our mixed FEM with a time-stepping
scheme that is based on a convolution quadrature in time generated by the
backward Euler method in Section 4. A pointwise-in-time optimal $L^{2}$-error
estimate is proved in Theorem 10. Finally, in Section 5, numerical experiments
are presented to confirm our theoretical findings.

##  2 Preliminaries

We shall first introduce notation and recall some preliminary results. Let
$(\cdot,\cdot)$ be the inner product in $L^{2}(\Omega)$ and $\|\cdot\|$ the
induced norm. For a nonnegative integer $m$, the space $H^{m}(\Omega)$ denotes
the standard Sobolev space with the usual norm $\|\cdot\|_{H^{m}(\Omega)}$.
The norm $\|\cdot\|_{\dot{H}^{r}(\Omega)}$ in the (fractional-order) Sobolev
space $\dot{H}^{r}(\Omega)$ is defined in the usual way via the Dirichlet
eigenfunctions of the self-adjoint elliptic operator
$-\nabla\cdot(\kappa\nabla)$ on $\Omega$, see [32]. Therein, for a constant
$\kappa,$ it is proved that when $r$ is a nonnegative integer,
$\dot{H}^{r}(\Omega)$ consists of all functions $v$ in $H^{r}(\Omega)$ which
satisfy the boundary conditions $\Delta^{j}v=0$ on $\partial\Omega$ for
$j<r/2$, and that the norms $\|\cdot\|_{\dot{H}^{r}(\Omega)}$ and
$\|\cdot\|_{r}$ are equivalent.

Regularity properties of the solution $u$ of the time-fractional problem (1)
play a key role in the error analysis of numerical methods, particularly,
since $u$ has singularity near $t=0$, even for smooth given data. From [24,
Theorems 11, 12 and 13] or [26, Theorem 3], we recall the following regularity
results: for $r\in[0,2]$, $u$ and its time derivative satisfy

$\|u(t)\|_{2}+t\|u^{\prime}(t)\|_{2}\leq
Ct^{-\alpha(2-r)/2}\|u_{0}\|_{\dot{H}^{r}(\Omega)},\quad t>0,$ (3)

where ′ denotes the time partial derivative. Properties in (3) are used in the
proofs of the estimates (18), (20), as well as in the proof of first and third
estimates of Theorem 7. For the error bound of the secondary variable, to show
a second order of accuracy, we require that

$\|u(t)\|_{3}\leq Ct^{-\alpha(3-r)/2}\|u_{0}\|_{\dot{H}^{r}(\Omega)},\quad
r\in[1,2].$ (4)

This additional assumption is used in (19) and in the proof of the second
estimate of Theorem 7. When ${\bf F}\equiv 0,$ (4) is proved in [23, Theorem
4.1].

Nevertheless, the above regularity estimate holds in the limiting case
$\alpha=1$, when the model problem (1) reduces to the classical Fokker–Planck
PDE; see [32, Lemmas 3.2 and 4.4]. In (3) and throughout the paper, $C$
denotes a generic positive constant that may depend on $T$, $\Omega$,
$\kappa$, and ${\bf F}$.

We show next some properties of $\mathcal{I}^{\alpha}$ for later use in the
subsequent sections.

######  Lemma 1.

If the function $\varphi\in L^{2}((0,t),L^{2}(\Omega))$, then

$\int_{0}^{t}\|\mathcal{I}^{\alpha}\varphi\|^{2}\,ds\leq\omega_{\alpha+1}(t)\int_{0}^{t}\omega_{\alpha}(t-s)\
\mathcal{I}(\|\varphi\|^{2})(s),\quad{\rm for~{}any}~{}~{}0\leq\alpha\leq 1,$
(5) $\int_{0}^{t}\|\mathcal{I}^{\alpha}\varphi\|^{2}\,ds\leq
Ct^{2(\alpha-\mu)}\int_{0}^{t}\|\mathcal{I}^{\mu}\varphi\|^{2}\,ds,\quad{\rm
for~{}any}~{}~{}0\leq\mu\leq\alpha\leq 1,$ (6)
$\mathcal{I}^{\alpha}(\|\mathcal{I}\varphi\|^{2})(t)\leq\omega_{2-\nu}(t)\mathcal{I}^{1+\alpha-\nu}(\|\mathcal{I}^{\nu}\varphi\|^{2})(t),\quad{\rm
for~{}any}~{}~{}0\leq\alpha,\nu\leq 1\,.$ (7)

###### Proof.

The proof of (5) is in [19, Lemma 2.3]. For the proof of (6), we refer to [20,
Lemma 3.1]. To show (7), we use first that
$\mathcal{I}\varphi(t)=\mathcal{I}^{1-\nu}(\mathcal{I}^{\nu}\varphi)(t)$ for
$0\leq\nu\leq 1,$ and notice that

$\|\mathcal{I}\varphi(t)\|=\|\mathcal{I}^{1-\nu}(\mathcal{I}^{\nu}\varphi)(t)\|\leq\mathcal{I}^{1-\nu}(\|\mathcal{I}^{\nu}\varphi\|)(t)\,,$

and hence, applying the Cauchy-Schwarz inequality and integrating, we obtain

$\|\mathcal{I}\varphi(t)\|^{2}\leq\int_{0}^{t}\omega_{1-\nu}(t-s)\,ds\int_{0}^{t}\omega_{1-\nu}(t-s)\|\mathcal{I}^{\nu}\varphi(s)\|^{2}\,ds=\omega_{2-\nu}(t)\mathcal{I}^{1-\nu}(\|\mathcal{I}^{\nu}\varphi\|^{2})(t)\,.$

Applying $\mathcal{I}^{\alpha}$ to both sides yields

$\mathcal{I}^{\alpha}(\|\mathcal{I}\varphi\|^{2})(t)\leq\omega_{2-\nu}(t)\mathcal{I}^{\alpha}\Big{(}\mathcal{I}^{1-\nu}(\|\mathcal{I}^{\nu}\varphi\|^{2})\Big{)}(t)=\omega_{2-\nu}(t)\mathcal{I}^{1+\alpha-\nu}(\|\mathcal{I}^{\nu}\varphi\|^{2})(t),$

and the proof of (7) is completed. ∎

Before ending this section, we recall the following identity which follows
from the generalized Leibniz formula:

$\partial_{t}^{\alpha}(t\varphi)=t\partial_{t}^{\alpha}\varphi+\alpha\mathcal{I}^{1-\alpha}\varphi+t\omega_{1-\alpha}(t)\varphi(0)\,.$
(8)

##  3 Mixed Finite Element Method

In this section, we consider the mixed form of problem (2). Introduce
${\boldsymbol{\sigma}}=-\kappa\nabla u+{\bf F}u$ so that
$\kappa^{-1}{\boldsymbol{\sigma}}+\nabla u=\beta u$ where
$\beta=\kappa^{-1}{\bf F}.$ Our model problem can now be formulated as

$\partial_{t}^{\alpha}u+\nabla\cdot{\boldsymbol{\sigma}}=0,\quad u=0\;\mbox{
on }\partial\Omega,\quad{\rm with}~{}~{}u(0)=u_{0}.$

Let $H(div;\Omega)=\\{{\bf w}\in(L^{2}(\Omega))^{2}:\nabla\cdot{\bf w}\in
L^{2}(\Omega)\\}$ be a Hilbert space equipped with the norm $\|{\bf
w}\|_{H(div;\Omega)}=(\|{\bf w}\|^{2}+\|\nabla\cdot{\bf
w}\|^{2})^{\frac{1}{2}}$. Then, with $V=L^{2}(\Omega)$ and ${\bf
W}=H(div;\Omega)$, the mixed formulation of (1) is defined as follows: find
$(u,{\boldsymbol{\sigma}}):(0,T]\to V\times{\bf W}$ such that

$\displaystyle(\partial_{t}^{\alpha}u,v)+(\nabla\cdot{\boldsymbol{\sigma}},v)$
$\displaystyle=$ $\displaystyle 0\quad\forall v\in V,$ (9)
$\displaystyle(\kappa^{-1}{\boldsymbol{\sigma}},{\bf w})-(u,\nabla\cdot{\bf
w})$ $\displaystyle=$ $\displaystyle(\beta u,{\bf w})\quad\forall{\bf
w}\in{\bf W},$ (10)

with $u(0)=u_{0}$. Note that the homogeneous Dirichlet boundary conditions of
$u$ on $\partial\Omega$ are implicitly contained in (10).

For the semidiscrete mixed formulation corresponding to (9)-(10), consider a
shape-regular partition ${\mathcal{T}}_{h}$ of the polygonal convex domain
$\bar{\Omega}$ into triangles $K$, each with diameter $h_{K}$. Additionally,
choose appropriate finite element subspaces $V_{h}$ and ${\bf W}_{h}$ of $V$
and ${\bf W}$, respectively, that satisfy the Ladyzenskaya-Babuska-Brezzi
(LBB) condition. As an example, let’s consider using the Raviart-Thomas (RT)
spaces [30] with an index $\ell\geq 0$, which are defined by:

$V_{h}=\\{v\in L^{2}(\Omega):\;v|_{K}\in P_{\ell}(K)\;\forall
K\in{\mathcal{T}}_{h}\\}$

and

${\bf W}_{h}=\\{{\bf v}\in{\bf W}:\;{\bf v}|_{K}\in RT_{\ell}(K)\;\forall
K\in{\mathcal{T}}_{h}\\},$

where
$RT_{\ell}(K)=(P_{\ell}(K))^{2}+{\mathcal{\mathbb{missing}}x}P_{\ell}(K),$
${\ell}\geq 0$, and $P_{\ell}(K)$ is the space of polynomials of degree
$\leq\ell$ on $K$. For each of these mixed spaces, there exists a projection
operator $\Pi_{h}:{\bf W}\to{\bf W}_{h}$ known as the Fortin projection, which
satisfies the property $\nabla\cdot\Pi_{h}=P_{h}(\nabla\cdot)$, where
$P_{h}:V\to V_{h}$ is the $L^{2}$-projection defined by $(P_{h}v-v,v_{h})=0$
for all $v_{h}\in V_{h}.$ Further, the following approximation properties
hold:

$||{\bf w}-\Pi_{h}{\bf w}||\leq Ch^{r}||\nabla\cdot{\bf w}||_{r-1}~{}~{}{\rm
and}~{}~{}\|v-P_{h}v\|\leq Ch^{r}\|v\|_{r}\;\forall~{}1\leq r\leq\ell+1.$ (11)

Additionally, it’s worth noting that $\Pi_{h}$ and $P_{h}$ satisfy the
following conditions:

$(\nabla\cdot(\Pi_{h}{\bf w}-{\bf w}),v_{h})=0\;\forall v_{h}\in
V_{h}\quad\text{and}\quad(P_{h}v-v,\nabla\cdot{\bf w}_{h})=0\;\forall{\bf
w}_{h}\in{\bf W}_{h}.$

For additional examples of these spaces, including Brezzi-Douglas-Marini
spaces and Brezzi-Douglas-Fortin-Marini spaces, refer to [2]. Due to the
limited smoothing properties of the continuous problem, we shall focus on the
case where $\ell=1$, as high-order elements are not relevant, as discussed in
[15].

The corresponding semidiscrete mixed finite element approximation of problem
(1) is to seek a pair $(u_{h},{\boldsymbol{\sigma}}_{h}):(0,T]\to
V_{h}\times{\bf W}_{h}$ such that

$\displaystyle(\partial_{t}^{\alpha}u_{h},v_{h})+(\nabla\cdot{\boldsymbol{\sigma}}_{h},v_{h})$
$\displaystyle=$ $\displaystyle 0\quad\forall v_{h}\in V_{h},$ (12)
$\displaystyle(\kappa^{-1}{\boldsymbol{\sigma}}_{h},{\bf
w}_{h})-(u_{h},\nabla\cdot{\bf w}_{h})$ $\displaystyle=$ $\displaystyle(\beta
u_{h},{\bf w}_{h})\quad\forall{\bf w}_{h}\in{\bf W}_{h},$ (13)

with $u_{h}(0)=u_{0h}$, where $u_{0h}$ an appropriate approximation of $u_{0}$
in $V_{h}$.

As both $V_{h}$ and ${\bf W}_{h}$ are finite-dimensional spaces, the system
(12)-(13) results in a system of time-fractional linear ODEs. Consequently,
the theory developed for the linear system in [18] guarantees the existence of
a unique solution to the semidiscrete system.

To derive optimal errors from the above mixed finite element discretization
for the case of smooth and nonsmooth initial data, we rely on the delicate
energy argument approach. Since the problem has a limited smoothing property,
integration in time with a $t$ type weight is an essential tool to provide
optimal error estimates. This idea has been used in [14] and [28] to derive
optimal error bounds of the standard continuous Galerkin method applied to
time-fractional diffusion equations with Riemann–Liouville and Caputo
derivatives, respectively. The approach in [14, 28] extends the argument in
[21] for parabolic problems to the time-fractional case. A similar approach
applied to mixed finite element methods for parabolic problems has also been
exploited in [9]. For convenience, define $e_{u}=u_{h}-u$ and
$e_{\boldsymbol{\sigma}}={\boldsymbol{\sigma}}_{h}-{\boldsymbol{\sigma}}$, and
introduce the intermediate solutions
$(\bar{u}_{h},\bar{\boldsymbol{\sigma}}_{h}):(0,T_{0}]\to V_{h}\times{\bf
W}_{h}$ defined by

$\displaystyle(\partial_{t}^{\alpha}\bar{u}_{h},v_{h})+(\nabla\cdot\bar{\boldsymbol{\sigma}}_{h},v_{h})$
$\displaystyle=$ $\displaystyle 0\quad\forall v_{h}\in V_{h},$ (14)
$\displaystyle(\kappa^{-1}\bar{\boldsymbol{\sigma}}_{h},{\bf
w}_{h})-(\bar{u}_{h},\nabla\cdot{\bf w}_{h})$ $\displaystyle=$
$\displaystyle(\beta u,{\bf w}_{h})\quad\forall{\bf w}_{h}\in{\bf W}_{h},$
(15)

with $\bar{u}_{h}(0)=u_{0h}$. Setting $\tilde{\phi}_{u}=\bar{u}_{h}-u$ and
$\bar{\phi}_{\boldsymbol{\sigma}}=\bar{\boldsymbol{\sigma}}_{h}-{\boldsymbol{\sigma}}$,
we split the errors $e_{u}=\tilde{\phi}_{u}+\bar{e}_{u}$ and
$e_{\boldsymbol{\sigma}}=\bar{\phi}_{\boldsymbol{\sigma}}+\bar{e}_{\boldsymbol{\sigma}}$,
where $\bar{e}_{u}:=u_{h}-\bar{u}_{h}$ and
$\bar{e}_{\boldsymbol{\sigma}}:={\boldsymbol{\sigma}}_{h}-\bar{\boldsymbol{\sigma}}_{h}$.
Thus, from (9)-(10) and (14)-(15), we obtain

$\displaystyle(\partial_{t}^{\alpha}\tilde{\phi}_{u},v_{h})+(\nabla\cdot\bar{\phi}_{\boldsymbol{\sigma}},v_{h})$
$\displaystyle=$ $\displaystyle 0\quad\forall v_{h}\in V_{h},$ (16)
$\displaystyle(\kappa^{-1}\bar{\phi}_{\boldsymbol{\sigma}},{\bf
w}_{h})-(\tilde{\phi}_{u},\nabla\cdot{\bf w}_{h})$ $\displaystyle=$
$\displaystyle 0\quad\forall{\bf w}_{h}\in{\bf W}_{h}.$ (17)

For the proof of the next theorem, we refer to [15, Theorems 5.4 and 5.6].

######  Theorem 2.

Let $u_{0}\in\dot{H}^{\delta}(\Omega)$ with $\delta\in[0,2]$. Then, for $t>0$,

$\|\tilde{\phi}_{u}(t)\|+t\|\tilde{\phi}_{u}^{\prime}(t)\|\leq
Ch^{r}t^{-\alpha(r-\delta)/2}\|u_{0}\|_{\dot{H}^{\delta}(\Omega)},\quad\mbox{for
}\delta\in[0,2],$ (18)

with $r=1,2,$ and

$\|\bar{\phi}_{\boldsymbol{\sigma}}(t)\|\leq
Ch^{2}t^{-\alpha(3-\delta)/2}\|u_{0}\|_{\dot{H}^{\delta}(\Omega)},\quad\mbox{for
}\delta\in[1,2],$ (19) $\|\bar{\phi}_{\boldsymbol{\sigma}}(t)\|\leq
Cht^{-\alpha(2-\delta)/2}\|u_{0}\|_{\dot{H}^{\delta}(\Omega)},\quad\mbox{for
}\delta\in[0,1).$ (20)

The main task now is to estimate $\bar{e}_{u}$ and
$\bar{e}_{\boldsymbol{\sigma}}$. We start by deriving a preliminary bound for
$\bar{e}_{\boldsymbol{\sigma}}$ in the next lemma.

######  Lemma 3.

For $t>0,$ we have

$\int_{0}^{t}\|\mathcal{I}^{\alpha}\bar{e}_{\boldsymbol{\sigma}}\|^{2}\leq
C\int_{0}^{t}\|\mathcal{I}^{\alpha}\tilde{\phi}_{u}\|^{2}\,ds\,.$

###### Proof.

From (12)-(13) and (14)-(15), we have

$\displaystyle(\partial_{t}^{\alpha}\bar{e}_{u},v_{h})+(\nabla\cdot\bar{e}_{\boldsymbol{\sigma}},v_{h})$
$\displaystyle=$ $\displaystyle 0\quad\forall v_{h}\in V_{h},$ (21)
$\displaystyle(\kappa^{-1}\bar{e}_{\boldsymbol{\sigma}},{\bf
w}_{h})-(\bar{e}_{u},\nabla\cdot{\bf w}_{h})$ $\displaystyle=$
$\displaystyle(\beta e_{u},{\bf w}_{h})\quad\forall{\bf w}_{h}\in{\bf W}_{h}.$
(22)

Now, applying $\mathcal{I}^{\alpha}$ to both sides of (21), then integrating
and using $\bar{e}_{u}(0)=0$,

$(\bar{e}_{u},v_{h})+(\nabla\cdot\mathcal{I}^{\alpha}\bar{e}_{\boldsymbol{\sigma}},v_{h})=0\,.$
(23)

Choosing $v_{h}=\mathcal{I}^{\alpha}\bar{e}_{u}$ in (23) and ${\bf
w}_{h}=\mathcal{I}^{\alpha}\bar{e}_{\boldsymbol{\sigma}}$ in (22) after
applying $\mathcal{I}^{\alpha}$ to both sides, then adding them up and
applying the Cauchy-Schwarz inequality,

$(\bar{e}_{u},\mathcal{I}^{\alpha}\bar{e}_{u})+\|\sqrt{\kappa^{-1}}\mathcal{I}^{\alpha}\bar{e}_{\boldsymbol{\sigma}}\|^{2}=(\beta\mathcal{I}^{\alpha}e_{u},\mathcal{I}^{\alpha}\bar{e}_{\boldsymbol{\sigma}})\leq
C\|\mathcal{I}^{\alpha}e_{u}\|^{2}+\frac{1}{2}\|\sqrt{\kappa^{-1}}\mathcal{I}^{\alpha}\bar{e}_{\boldsymbol{\sigma}}\|^{2}.$

Canceling the common terms, then integrating both sides, using
$\int_{0}^{t}(\bar{e}_{u},\mathcal{I}^{\alpha}\bar{e}_{u})\,ds\geq 0,$ and the
diffusivity coefficient assumption $\kappa>c_{0}>0,$ we reach

$\int_{0}^{t}\|\mathcal{I}^{\alpha}\bar{e}_{\boldsymbol{\sigma}}\|^{2}\,ds\leq
C\int_{0}^{t}\|\mathcal{I}^{\alpha}e_{u}\|^{2}\,ds\,.$ (24)

Applying $\mathcal{I}^{\alpha}$ to both sides of (22) and (23), then choosing
$v_{h}=\mathcal{I}^{\alpha}\bar{e}_{u}$ and ${\bf
w}_{h}=\mathcal{I}^{2\alpha}\bar{e}_{\boldsymbol{\sigma}}$. Adding them up to
obtain

$\|\mathcal{I}^{\alpha}\bar{e}_{u}\|^{2}+(\kappa^{-1}\mathcal{I}^{\alpha}(\mathcal{I}^{\alpha}\bar{e}_{\boldsymbol{\sigma}}),\mathcal{I}^{\alpha}\bar{e}_{\boldsymbol{\sigma}})=(\beta\mathcal{I}^{\alpha}e_{u},\mathcal{I}^{2\alpha}\bar{e}_{\boldsymbol{\sigma}})\,.$

Since $e_{u}=\bar{e}_{u}+\tilde{\phi}_{u}$,
$\|\mathcal{I}^{\alpha}e_{u}\|^{2}\leq
2\|\mathcal{I}^{\alpha}\bar{e}_{u}\|^{2}+2\|\mathcal{I}^{\alpha}\tilde{\phi}_{u}\|^{2}$.
Hence, by using the above equation, we get

$\|\mathcal{I}^{\alpha}e_{u}\|^{2}+2(\kappa^{-1}\mathcal{I}^{\alpha}(\mathcal{I}^{\alpha}\bar{e}_{\boldsymbol{\sigma}}),\mathcal{I}^{\alpha}\bar{e}_{\boldsymbol{\sigma}})\leq
2\|\mathcal{I}^{\alpha}\tilde{\phi}_{u}\|^{2}+2(\beta\mathcal{I}^{\alpha}e_{u},\mathcal{I}^{2\alpha}\bar{e}_{\boldsymbol{\sigma}})\,.$

Using
$2(\beta\mathcal{I}^{\alpha}e_{u},\mathcal{I}^{2\alpha}\bar{e}_{\boldsymbol{\sigma}})\leq\frac{1}{2}\|\mathcal{I}^{\alpha}e_{u}\|^{2}+C\|\mathcal{I}^{2\alpha}\bar{e}_{\boldsymbol{\sigma}}\|^{2}$,
then canceling the common terms and integrating both sides, and since
$\int_{0}^{t}(\kappa^{-1}\mathcal{I}^{\alpha}(\mathcal{I}^{\alpha}\bar{e}_{\boldsymbol{\sigma}}),\mathcal{I}^{\alpha}\bar{e}_{\boldsymbol{\sigma}})\,ds\geq
0$, we reach

$\int_{0}^{t}\|\mathcal{I}^{\alpha}e_{u}\|^{2}\,ds\leq
C\int_{0}^{t}\|\mathcal{I}^{\alpha}\tilde{\phi}_{u}\|^{2}\,ds+C\int_{0}^{t}\|\mathcal{I}^{2\alpha}\bar{e}_{\boldsymbol{\sigma}}\|^{2}\,ds\,.$

Using (5) with $\varphi=\mathcal{I}^{\alpha}\bar{e}_{\boldsymbol{\sigma}},$ we
have

$\int_{0}^{t}\|\mathcal{I}^{2\alpha}\bar{e}_{\boldsymbol{\sigma}}\|^{2}\,ds\leq\omega_{\alpha+1}(t)\int_{0}^{t}\omega_{\alpha}(t-s)\int_{0}^{s}\|\mathcal{I}^{\alpha}\bar{e}_{\boldsymbol{\sigma}}(q)\|^{2}dq\,ds,$

and then, with the help of (24), we deduce that

$\int_{0}^{t}\|\mathcal{I}^{\alpha}\bar{e}_{\boldsymbol{\sigma}}\|^{2}\,ds\leq
C\int_{0}^{t}\|\mathcal{I}^{\alpha}\tilde{\phi}_{u}\|^{2}\,ds+\omega_{\alpha+1}(t)\int_{0}^{t}\omega_{\alpha}(t-s)\int_{0}^{s}\|\mathcal{I}^{\alpha}\bar{e}_{\boldsymbol{\sigma}}(q)\|^{2}dq\,ds.$

An application of a weakly singular Gronwall inequality completes the proof. ∎

In the next lemma, we derive an upper bound for $\bar{e}_{u}$ and
$\bar{e}_{\boldsymbol{\sigma}}$. This bound leads to optimal convergence rates
in the $L^{2}(\Omega)$-norm of $e_{u}$ and $e_{\boldsymbol{\sigma}}$.

######  Lemma 4.

For $0<t\leq T$, we have

$\|t\bar{e}_{u}\|^{2}\leq C\max_{0\leq s\leq
t}\mathcal{I}^{\alpha}(\|s\tilde{\phi}_{u}\|^{2})+C\,t^{1-\alpha}\int_{0}^{t}\|\mathcal{I}^{\alpha}\tilde{\phi}_{u}\|^{2}\,ds,$

and

$t^{2}\|\bar{e}_{\boldsymbol{\sigma}}\|^{2}\leq C\max_{0\leq s\leq
t}\Big{(}\mathcal{I}^{\alpha}(\|\mathcal{I}^{1-\alpha}(s\tilde{\phi}_{u})^{\prime}\|^{2})+\mathcal{I}^{\alpha}(\|\mathcal{I}^{1-\alpha}\bar{e}_{u}\|^{2})(s)\Big{)}\,.$

###### Proof.

Multiply both sides of (21) and (22) by $t$ and use (8) to find that

$\displaystyle(\partial_{t}^{\alpha}(t\bar{e}_{u}),v_{h})+(\nabla\cdot(t\bar{e}_{\boldsymbol{\sigma}}),v_{h})$
$\displaystyle=$
$\displaystyle\alpha(\mathcal{I}^{1-\alpha}\bar{e}_{u},v_{h}),$ (25)
$\displaystyle(\kappa^{-1}(t\bar{e}_{\boldsymbol{\sigma}}),{\bf
w}_{h})-(t\bar{e}_{u},\nabla\cdot{\bf w}_{h})$ $\displaystyle=$
$\displaystyle(\beta(te_{u}),{\bf w}_{h}).$ (26)

Next choose $v_{h}=t\bar{e}_{u}$ and ${\bf
w}_{h}=t\bar{e}_{\boldsymbol{\sigma}}$ to obtain

$(\partial_{t}^{\alpha}(t\bar{e}_{u}),t\bar{e}_{u})+\|\sqrt{\kappa^{-1}}\,(t\bar{e}_{\boldsymbol{\sigma}})\|^{2}=(\beta(te_{u}),t\bar{e}_{\boldsymbol{\sigma}})+\alpha(\mathcal{I}^{1-\alpha}\bar{e}_{u},t\bar{e}_{u}).$

By integrating (21) in time followed by choosing $v_{h}=t\bar{e}_{u}$ then
choosing ${\bf w}_{h}=t\mathcal{I}\bar{e}_{\boldsymbol{\sigma}}$ in (22) and
adding them up. This leads to

$(\mathcal{I}^{1-\alpha}\bar{e}_{u},t\bar{e}_{u})=-(\kappa^{-1}t\bar{e}_{\boldsymbol{\sigma}},\mathcal{I}\bar{e}_{\boldsymbol{\sigma}})+(\beta
te_{u},\mathcal{I}\bar{e}_{\boldsymbol{\sigma}}),$ (27)

and hence,

$(\partial_{t}^{\alpha}(t\bar{e}_{u}),t\bar{e}_{u})+\|\sqrt{\kappa^{-1}}\,(t\bar{e}_{\boldsymbol{\sigma}})\|^{2}=(\alpha+1)(\beta(te_{u}),t\bar{e}_{\boldsymbol{\sigma}})-\alpha(\kappa^{-1}t\bar{e}_{\boldsymbol{\sigma}},\mathcal{I}\bar{e}_{\boldsymbol{\sigma}}).$

Using $\partial_{t}^{\alpha}(\|t\bar{e}_{u}\|^{2})\leq
2(\partial_{t}^{\alpha}(t\bar{e}_{u}),t\bar{e}_{u})$, which is due to [1,
Corollary 1], we get

$\partial_{t}^{\alpha}(\|t\bar{e}_{u}\|^{2})+2\|\sqrt{\kappa^{-1}}\,(t\bar{e}_{\boldsymbol{\sigma}})\|^{2}\leq
C\|te_{u}\|^{2}+\|\sqrt{\kappa^{-1}}\,(t\bar{e}_{\boldsymbol{\sigma}})\|^{2}+C\|\mathcal{I}\bar{e}_{\boldsymbol{\sigma}}\|^{2}.$
(28)

Canceling the similar terms, then using $e_{u}=\tilde{\phi}_{u}+\bar{e}_{u}$
and applying $\mathcal{I}^{\alpha}$,

$\displaystyle\|t\bar{e}_{u}\|^{2}$ $\displaystyle\leq
C\mathcal{I}^{\alpha}(\|t\tilde{\phi}_{u}\|^{2})+C\mathcal{I}^{\alpha}(\|t\bar{e}_{u}\|^{2})+C\mathcal{I}^{\alpha}(\|\mathcal{I}\bar{e}_{\boldsymbol{\sigma}}\|^{2})\,.$

By applying (7) with $\nu=\alpha$, we have
$\mathcal{I}^{\alpha}(\|\mathcal{I}\bar{e}_{\boldsymbol{\sigma}}\|^{2})(t)\leq\omega_{2-\alpha}(t)\mathcal{I}(\|\mathcal{I}^{\alpha}\bar{e}_{\boldsymbol{\sigma}}(s)\|^{2})(t),$
and so, using Lemma 3, we see that

$\displaystyle\|t\bar{e}_{u}\|^{2}$ $\displaystyle\leq C\max_{0\leq s\leq
t}\mathcal{I}^{\alpha}(\|s\tilde{\phi}_{u}\|^{2})+C\,t^{1-\alpha}\int_{0}^{t}\|\mathcal{I}^{\alpha}\tilde{\phi}_{u}\|^{2}\,ds+C\mathcal{I}^{\alpha}\Big{(}\|t\bar{e}_{u}\|^{2}\Big{)}.$

Consequently, an application of a weakly singular Gronwall inequality leads to
the first desired estimate for $\bar{e}_{u}$. We now show the estimate for
$\bar{e}_{\boldsymbol{\sigma}}$. To do so, we operate $\partial_{t}^{\alpha}$
on (26) so that

$(\kappa^{-1}\partial_{t}^{\alpha}(t\bar{e}_{\boldsymbol{\sigma}}),{\bf
w}_{h})-(\partial_{t}^{\alpha}(t\bar{e}_{u}),\nabla\cdot{\bf
w}_{h})=(\beta\partial_{t}^{\alpha}(te_{u}),{\bf w}_{h}).$ (29)

Next choose $v_{h}=\partial_{t}^{\alpha}(t\bar{e}_{u})$ in (25) and ${\bf
w}_{h}=t\bar{e}_{\boldsymbol{\sigma}}$ in (29), and arrive at

$\|\partial_{t}^{\alpha}(t\bar{e}_{u})\|^{2}+\frac{1}{2}\partial_{t}^{\alpha}\|\sqrt{\kappa^{-1}}t\bar{e}_{\boldsymbol{\sigma}}\|^{2}\leq
C\left(\|\partial_{t}^{\alpha}(t\tilde{\phi}_{u})\|^{2}+\|\mathcal{I}^{1-\alpha}\bar{e}_{u}\|^{2}+\|\sqrt{\kappa^{-1}}t\bar{e}_{\boldsymbol{\sigma}}\|^{2}\right)+\frac{1}{2}\|\partial_{t}^{\alpha}(t\bar{e}_{u})\|^{2},$

and thus,

$\partial_{t}^{\alpha}\|\sqrt{\kappa^{-1}}t\bar{e}_{\boldsymbol{\sigma}}\|^{2}\leq
C\left(\|\partial_{t}^{\alpha}(t\tilde{\phi}_{u})\|^{2}+\|\mathcal{I}^{1-\alpha}\bar{e}_{u}\|^{2}+\|\sqrt{\kappa^{-1}}t\bar{e}_{\boldsymbol{\sigma}}\|^{2}\right).$

Applying $\mathcal{I}^{\alpha}$ to both sides and using
$t\bar{e}_{\boldsymbol{\sigma}}(0)=0$, we get

$t^{2}\|\sqrt{\kappa^{-1}}\bar{e}_{\boldsymbol{\sigma}}\|^{2}\leq C\max_{0\leq
s\leq
t}\Big{(}\mathcal{I}^{\alpha}(\|\mathcal{I}^{1-\alpha}(s\tilde{\phi}_{u})^{\prime}\|^{2})+C\mathcal{I}^{\alpha}(\|\mathcal{I}^{1-\alpha}\bar{e}_{u}\|^{2})(s)\Big{)}\\\
+C\mathcal{I}^{\alpha}\left(\|\sqrt{\kappa^{-1}}t\bar{e}_{\boldsymbol{\sigma}}\|^{2}\right)\,.$

For completing the proof of the second desired estimate, apply a weakly
singular Gronwall inequality then use the assumption $\kappa>c_{0}>0$. ∎

We are ready now to derive optimal error estimates for the semidiscrete mixed
finite element problem with smooth and nonsmooth initial data.

######  Theorem 5.

Let the pairs $(u,{\boldsymbol{\sigma}})$ and
$(u_{h},{\boldsymbol{\sigma}}_{h})$ be the solutions of (9)-(10) and
(12)-(13), respectively. Assume $u_{0}\in\dot{H}^{\delta}(\Omega)$ with
$\delta\in[0,2]$, and $u_{0h}=P_{h}u_{0}$. Then the following error estimates
hold for $t>0$:

$\|(u_{h}-u)(t)\|\leq
Ch^{2}t^{-\alpha(2-\delta)/2}\|u_{0}\|_{\dot{H}^{\delta}(\Omega)},\quad\delta\in[0,2],$
(30)

and

$\|({\boldsymbol{\sigma}}_{h}-{\boldsymbol{\sigma}})(t)\|\leq
Ch^{2}t^{-\alpha(3-\delta)/2}\|u_{0}\|_{\dot{H}^{\delta}(\Omega)},\quad\delta\in[1,2],$
(31) $\|({\boldsymbol{\sigma}}_{h}-{\boldsymbol{\sigma}})(t)\|\leq
Cht^{-\alpha(2-\delta)/2}\|u_{0}\|_{\dot{H}^{\delta}(\Omega)},\quad\delta\in[0,1).$
(32)

###### Proof.

Using the estimate (18), we find that for $\delta\in[0,2]$,

$\mathcal{I}^{\alpha}(\|t\tilde{\phi}_{u}\|^{2})\leq
Ch^{4}\int_{0}^{t}(t-s)^{\alpha-1}s^{-\alpha(2-\delta)+2}\,ds\|u_{0}\|_{\dot{H}^{\delta}(\Omega)}\leq
Ch^{4}t^{\alpha(\delta-1)+2}\|u_{0}\|_{\dot{H}^{\delta}(\Omega)}^{2},$

and

$\int_{0}^{t}\|\mathcal{I}^{\alpha}\tilde{\phi}_{u}\|^{2}\,ds\leq
Ch^{4}\int_{0}^{t}s^{\alpha\delta}\,ds\|u_{0}\|_{\dot{H}^{\delta}(\Omega)}\leq
Ch^{4}t^{\alpha\delta+1}\|u_{0}\|_{\dot{H}^{\delta}(\Omega)}^{2}\,.$

Then, by Lemma 4, we deduce that for $\delta\in[0,2]$,

$\|\bar{e}_{u}\|^{2}\leq
Ch^{4}\Big{(}t^{\alpha(\delta-1)}+t^{\alpha\delta}\Big{)}\|u_{0}\|_{\dot{H}^{\delta}(\Omega)}^{2}\leq
Ch^{4}t^{-\alpha(1-\delta)}\|u_{0}\|_{\dot{H}^{\delta}(\Omega)}^{2}\,.$

Similarly, for $\delta\in[1,2]$, we have

$\mathcal{I}^{\alpha}\Big{(}\|\mathcal{I}^{1-\alpha}(t\tilde{\phi}_{u})^{\prime}\|^{2}\Big{)}\leq
Ch^{4}\int_{0}^{t}(t-s)^{\alpha-1}s^{-\alpha(4-\delta)+2}\|u_{0}\|_{\dot{H}^{\delta}(\Omega)}^{2}\leq
Ch^{4}t^{-\alpha(3-\delta)+2}\|u_{0}\|_{\dot{H}^{\delta}(\Omega)}^{2},$

and

$\mathcal{I}^{\alpha}\left(\|\mathcal{I}^{1-\alpha}\bar{e}_{u}\|^{2}\right)(t)\leq\mathcal{I}^{\alpha}((\mathcal{I}^{1-\alpha}\|\bar{e}_{u}\|)^{2})(t)\\\
\leq
Ch^{4}\int_{0}^{t}(t-s)^{\alpha-1}s^{\alpha(\delta-3)+2}\,ds\|u_{0}\|_{\dot{H}^{\delta}(\Omega)}\leq
Ch^{4}t^{\alpha(\delta-2)+2}\|u_{0}\|_{\dot{H}^{\delta}(\Omega)}^{2}\,.$

Then, by Lemma 4, we deduce $\|\bar{e}_{\boldsymbol{\sigma}}\|^{2}\leq
Ch^{4}t^{-\alpha(3-\delta)}\|u_{0}\|_{\dot{H}^{\delta}(\Omega)}^{2},$ for
$\delta\in[1,2].$ The desired estimates (30) and (31) follow now by using the
triangle inequality and the bounds in Theorem 2. The last inequality (32) can
be obtained analogously by using the estimate (18) with $r=1$ and
$\delta\in[0,1)$. ∎

######  Remark 6.

Considering the approximation properties given in (11), it can be deduced that
the estimates (30) and (31) are both of order $O(h)$ in the case of the lowest
RT elements ($\ell=0$).

##  4 Fully discrete schemes

This section is devoted to the analysis of a fully discrete scheme for problem
(9)-(10) based on a convolution quadrature (CQ) generated by the backward
Euler method. Divide the time interval $[0,T]$ into $N$ equal subintervals
with a time step size $\tau=T/N$, and let $t_{j}=j\tau$.

The time fractional Riemann–Liouville derivative
${{}^{R}}\partial_{t}^{\alpha}(\varphi_{n})$ can be conveniently discretized
using the BECQ (with $\varphi_{j}=\varphi(t_{j})$):

$\partial_{\tau}^{\alpha}\varphi_{n}:=\tau^{-\alpha}\sum_{j=0}^{n}a_{n-j}^{(\alpha)}\varphi_{j},\quad\mbox{where}\;\sum_{j=0}^{\infty}a_{j}^{(\alpha)}\xi^{j}=(1-\xi)^{\alpha},\quad
a_{j}^{(\alpha)}=(-1)^{j}\left(\begin{array}[]{c}\alpha\\\
j\end{array}\right).$

We approximate the Riemann–Liouville integral $\mathcal{I}^{\alpha}$ using the
discrete operator $\partial_{\tau}^{-\alpha}$, which is similarly defined.
Note that $\partial_{\tau}^{-1}\varphi_{n}=\tau\sum_{j=0}^{n}\varphi_{j}$.
Using the relation
${{}^{R}}\partial_{t}^{\alpha}\varphi(t)=\partial_{t}^{\alpha}(\varphi(t)-\varphi(0))$,
the proposed fully implicit scheme for (9)-(10) is to find a pair
$(U_{h}^{n},\Sigma_{h}^{n})\in V_{h}\times{\bf W}_{h}$ such that for $n\geq
1$,

$\displaystyle(\partial_{\tau}^{\alpha}U_{h}^{n},\chi_{h})+(\nabla\cdot\Sigma_{h}^{n},\chi_{h})$
$\displaystyle=$
$\displaystyle(\partial_{\tau}^{\alpha}U_{h}^{0},\chi_{h})\;\;\;\forall\chi_{h}\in
V_{h},$ (33) $\displaystyle(\kappa^{-1}\Sigma_{h}^{n},{\bf
w}_{h})-(U_{h}^{n},\nabla\cdot{\bf w}_{h})$ $\displaystyle=$
$\displaystyle(\beta U_{h}^{n},{\bf w}_{h})\;\;\;\forall{\bf w}_{h}\in{\bf
W}_{h},$ (34)

with $U_{h}^{0}=P_{h}u_{0}$. Recall that $u_{0}\in\dot{H}^{\delta}(\Omega)$
and $\delta\in[0,2]$. For the error analysis, we introduce the intermediate
discrete solution $(\bar{U}_{h}^{n},\bar{\Sigma}_{h}^{n})\in V_{h}\times{\bf
W}_{h}$ satisfying

$\displaystyle(\partial_{\tau}^{\alpha}\bar{U}_{h}^{n},\chi_{h})+(\nabla\cdot\bar{\Sigma}_{h}^{n},\chi_{h})$
$\displaystyle=$
$\displaystyle(\partial_{\tau}^{\alpha}U_{h}^{0},\chi_{h})\;\;\;\forall\chi_{h}\in
V_{h},$ (35) $\displaystyle(\kappa^{-1}\bar{\Sigma}_{h}^{n},{\bf
w}_{h})-(\bar{U}_{h}^{n},\nabla\cdot{\bf w}_{h})$ $\displaystyle=$
$\displaystyle(\beta u_{h}^{n},{\bf w}_{h})\;\;\;\forall{\bf w}_{h}\in{\bf
W}_{h},$ (36)

for $n\geq 1$, with $\bar{U}_{h}^{0}=U_{h}^{0}$. The first estimate is the
next theorem is proved in [16, Theorem 6.2]. The estimate for
$\bar{\Sigma}_{h}^{n}$ is proved in [16, Theorem 6.3] for initial data
$u_{0}\in L^{2}(\Omega)$, i.e., for the case $\delta=0$. For $\delta>0$, one
can directly obtain the corresponding estimates from the proof of the same
theorem.

######  Theorem 7.

Let $(\bar{U}_{h}^{n},\bar{\Sigma}_{h}^{n})$ and
$(u_{h},{\boldsymbol{\sigma}}_{h})$ be the solutions of (35)-(36) and
(12)-(13), respectively. Then the following error estimates hold for
$t_{n}>0$:

$\|\bar{U}_{h}^{n}-u_{h}(t_{n})\|\leq C(\tau
t_{n}^{\alpha\delta/2-1}+h^{2}t_{n}^{-\alpha(2-\delta)/2})\|u_{0}\|_{\dot{H}^{\delta}(\Omega)},\quad\delta\in[0,2],$

and

$\|\bar{\Sigma}_{h}^{n}-{\boldsymbol{\sigma}}_{h}(t_{n})\|\leq C(\tau
t_{n}^{\alpha(\delta-1)/2-1}+h^{2}t_{n}^{-\alpha(3-\delta)/2})\|u_{0}\|_{\dot{H}^{\delta}(\Omega)},\quad\delta\in[1,2],$
$\|\bar{\Sigma}_{h}^{n}-{\boldsymbol{\sigma}}_{h}(t_{n})\|\leq C(\tau
t_{n}^{\alpha(\delta-1)/2-1}+ht_{n}^{-\alpha(2-\delta)/2})\|u_{0}\|_{\dot{H}^{\delta}(\Omega)},\quad\delta\in[0,1).$

The main task now is to estimate the errors $\bar{\cal
E}_{u}^{n}:=U_{h}^{n}-\bar{U}_{h}^{n}$ and $\bar{\cal
E}_{\boldsymbol{\sigma}}^{n}:=\Sigma_{h}^{n}-\bar{\Sigma}_{h}^{n}$. From
(33)-(34) and (35)-(36), we note that $\bar{\cal E}_{u}^{n}$ and $\bar{\cal
E}_{\boldsymbol{\sigma}}^{n}$ satisfy

$\displaystyle(\partial_{\tau}^{\alpha}\bar{\cal
E}_{u}^{n},\chi_{h})+(\nabla\cdot\bar{\cal
E}_{\boldsymbol{\sigma}}^{n},v_{h})$ $\displaystyle=$ $\displaystyle
0\;\;\;\forall v_{h}\in V_{h},$ (37) $\displaystyle(\kappa^{-1}\bar{\cal
E}_{\boldsymbol{\sigma}}^{n},{\bf w}_{h})-(\bar{\cal
E}_{u}^{n},\nabla\cdot{\bf w}_{h})$ $\displaystyle=$ $\displaystyle(\beta{\cal
E}_{u}^{n},{\bf w}_{h})\;\;\;\forall{\bf w}_{h}\in{\bf W}_{h}.$ (38)

For convenience, we set $\bar{\Phi}_{u}^{n}=\bar{U}_{h}^{n}-u_{h}(t_{n})$ and
$\bar{\Phi}_{\boldsymbol{\sigma}}^{n}=\bar{\Sigma}_{h}^{n}-\sigma_{h}(t_{n})$.
In the next lemma, we derive a preliminary bound for $\bar{\cal
E}_{\boldsymbol{\sigma}}$. Before, we recall the following properties (see,
[22] and [36], respectively):

$\frac{\tau^{\alpha}}{2}\sum_{j=0}^{n}\|\bar{\cal
E}_{u}^{n}\|^{2}\leq\sum_{j=0}^{n}(\partial_{\tau}^{-\alpha}\bar{\cal
E}_{u}^{n},\bar{\cal E}_{u}^{n}),$ (39)

and

$\frac{1}{2}\partial_{\tau}^{\alpha}\|\bar{\cal
E}_{u}^{n}\|^{2}\leq(\partial_{\tau}^{\alpha}\bar{\cal E}_{u}^{n},\bar{\cal
E}_{u}^{n}).$ (40)

######  Lemma 8.

For $n\geq 1,$ we have
$\partial_{\tau}^{-\alpha}(\|\partial_{\tau}^{-1}\bar{\cal
E}_{\boldsymbol{\sigma}}\|^{2})^{n}\leq
C\,t_{n}^{1-\alpha}\partial_{\tau}^{-1}(\|\partial_{\tau}^{-\alpha}\bar{\Phi}_{u}\|^{2})^{n}\,.$

###### Proof.

Using the Cauchy-Schwarz inequality, an appropriate simplification yields

$\|\partial_{\tau}^{-1}\bar{\cal
E}_{\boldsymbol{\sigma}}^{n}\|^{2}\leq\big{(}\partial_{\tau}^{\alpha-1}(\|\partial_{\tau}^{-\alpha}\bar{\cal
E}_{\boldsymbol{\sigma}}\|)^{n}\big{)}^{2}\leq
Ct_{n}^{1-\alpha}\partial_{\tau}^{\alpha-1}(\|\partial_{\tau}^{-\alpha}\bar{\cal
E}_{\boldsymbol{\sigma}}\|^{2})^{n},$

and thus,

$\partial_{\tau}^{-\alpha}(\|\partial_{\tau}^{-1}\bar{\cal
E}_{\boldsymbol{\sigma}}\|^{2})^{n}\leq
Ct_{n}^{1-\alpha}\partial_{\tau}^{-\alpha}\big{(}\partial_{\tau}^{\alpha-1}(\|\partial_{\tau}^{-\alpha}\bar{\cal
E}_{\boldsymbol{\sigma}}\|^{2})\big{)}^{n}=Ct_{n}^{1-\alpha}\partial_{\tau}^{-1}(\|\partial_{\tau}^{-\alpha}\bar{\cal
E}_{\boldsymbol{\sigma}}\|^{2})^{n}.$ (41)

Operating with $\partial_{\tau}^{-\alpha}$ on (37) and using
$\partial_{\tau}^{\alpha}\partial_{\tau}^{-\alpha}\bar{\cal
E}_{u}^{n}=\bar{\cal E}_{u}^{n}$, we get

$(\bar{\cal E}_{u}^{n},v_{h})+(\nabla\cdot\partial_{\tau}^{-\alpha}\bar{\cal
E}_{\boldsymbol{\sigma}}^{n},v_{h})=0.$ (42)

Choosing $v_{h}=\partial_{\tau}^{-\alpha}\bar{\cal E}_{u}^{n}$ in (42), and
${\bf w}_{h}=\partial_{\tau}^{-\alpha}\bar{\cal E}_{\boldsymbol{\sigma}}^{n}$
in (38) after operating with $\partial_{\tau}^{-\alpha}$ both sides, then
adding them up and applying the Cauchy-Schwarz inequality to obtain

$(\bar{\cal E}_{u}^{n},\partial_{\tau}^{-\alpha}\bar{\cal
E}_{u}^{n})+\|\sqrt{\kappa^{-1}}\partial_{\tau}^{-\alpha}\bar{\cal
E}_{\boldsymbol{\sigma}}^{n}\|^{2}=(\beta\partial_{\tau}^{-\alpha}{\cal
E}_{u}^{n},\partial_{\tau}^{-\alpha}\bar{\cal
E}_{\boldsymbol{\sigma}}^{n})\leq C\|\partial_{\tau}^{-\alpha}{\cal
E}_{u}^{n}\|^{2}+\frac{1}{2}\|\sqrt{\kappa^{-1}}\partial_{\tau}^{-\alpha}\bar{\cal
E}_{\boldsymbol{\sigma}}^{n}\|^{2}.$

Canceling the common terms, summing up over $n$, multiplying by $\tau,$ and
noting (39), we deduce

$\partial_{\tau}^{-1}(\|\sqrt{\kappa^{-1}}\partial_{\tau}^{-\alpha}\bar{\cal
E}_{\boldsymbol{\sigma}}\|^{2})^{n}\leq
C\partial_{\tau}^{-1}(\|\partial_{\tau}^{-\alpha}{\cal E}_{u}\|^{2})^{n},$

and thus, with
$g_{u}^{n}=\partial_{\tau}^{-1}(\|\partial_{\tau}^{-\alpha}\bar{\cal
E}_{u}\|^{2})^{n}$ and
$q_{u}^{n}=\partial_{\tau}^{-1}(\|\partial_{\tau}^{-\alpha}\bar{\Phi}_{u}\|^{2})^{n}$,
we have

$g_{\boldsymbol{\sigma}}^{n}\leq g_{u}^{n}+q_{u}^{n}\,.$ (43)

Operating with $\partial_{\tau}^{-\alpha}$ both sides of (38) and (42), then
choosing $v_{h}=\partial_{\tau}^{-\alpha}\bar{\cal E}_{u}^{n}$ and ${\bf
w}_{h}=\partial_{\tau}^{-2\alpha}\bar{\cal E}_{\boldsymbol{\sigma}}^{n}$.
Adding them up to obtain

$\|\partial_{\tau}^{-\alpha}\bar{\cal
E}_{u}^{n}\|^{2}+(\kappa^{-1}\partial_{\tau}^{-\alpha}(\partial_{\tau}^{-\alpha}\bar{\cal
E}_{\boldsymbol{\sigma}})^{n},\partial_{\tau}^{-\alpha}\bar{\cal
E}_{\boldsymbol{\sigma}}^{n})=(\beta\partial_{\tau}^{-\alpha}{\cal
E}_{u}^{n},\partial_{\tau}^{-2\alpha}\bar{\cal
E}_{\boldsymbol{\sigma}}^{n})\\\
\leq\frac{1}{2}\|\partial_{\tau}^{-\alpha}\bar{\cal
E}_{u}^{n}\|^{2}+\|\partial_{\tau}^{-\alpha}\bar{\Phi}_{u}^{n}\|^{2}+C\|\partial_{\tau}^{-2\alpha}\bar{\cal
E}_{\boldsymbol{\sigma}}^{n}\|^{2}\,.$

Canceling the common terms, summing over $n,$ and using (39), we reach

$g_{u}^{n}\leq
q_{u}^{n}+C\partial_{\tau}^{-1}(\|\partial_{\tau}^{-2\alpha}\bar{\cal
E}_{\boldsymbol{\sigma}}\|^{2})^{n}\,.$ (44)

Using

$\partial_{\tau}^{-1}(\|\partial_{\tau}^{-2\alpha}\bar{\cal
E}_{\boldsymbol{\sigma}}\|^{2})^{n}=\partial_{\tau}^{-1}(\|\partial_{\tau}^{-\alpha}(\partial_{\tau}^{-\alpha}\bar{\cal
E}_{\boldsymbol{\sigma}})\|^{2})^{n}\leq
Ct_{n}^{\alpha}\partial_{\tau}^{-1}(\partial_{\tau}^{-\alpha}(\|\partial_{\tau}^{-\alpha}\bar{\cal
E}_{\boldsymbol{\sigma}}\|^{2}))^{n}=Ct_{n}^{\alpha}\partial_{\tau}^{-\alpha}g_{\boldsymbol{\sigma}}^{n},$

and the estimates in (44) and (43), we deduce that
$g_{\boldsymbol{\sigma}}^{n}\leq
Cq_{u}^{n}+Ct_{n}^{\alpha}\partial_{\tau}^{-\alpha}g_{\boldsymbol{\sigma}}^{n}.$
Then, if the time step $\tau$ is such that
$Ca_{0}^{(-\alpha)}\tau^{\alpha}<1,$ an application of a weakly singular
discrete Gronwall inequality yields $g_{\boldsymbol{\sigma}}^{n}\leq
C\,q_{u}^{n}\,.$ Inserting this in (41) and the proof is completed. ∎

In the next lemma, we derive an upper bound for $\bar{\cal E}_{u}^{n}$ that
will be used later to show the $L^{2}(\Omega)$-norm convergence rate for
${\cal E}_{u}^{n}$. For brevity, for a given function $f,$ we set
$(tf)^{n}=t_{n}f(t_{n})=t_{n}f^{n}.$

######  Lemma 9.

For $n\geq 1$, we have

$\|t_{n}\bar{\cal E}_{u}^{n}\|^{2}\leq C\max_{1\leq j\leq
n}\Big{(}\partial_{\tau}^{-\alpha}(\|(t\bar{\Phi}_{u})^{j}\|^{2})+t_{j-1}^{1-\alpha}\partial_{\tau}^{-1}(\|(\partial_{\tau}^{-\alpha}\bar{\Phi}_{u})\|^{2})^{j-1}\Big{)}\,.$

###### Proof.

Multiply both sides of (37) and (38) by $t_{n}$ and use the identity

$t_{n}\partial_{\tau}^{\alpha}U_{h}^{n}=\partial_{\tau}^{\alpha}((tU_{h})^{n})-\alpha\partial_{\tau}^{\alpha-1}U_{h}^{n-1},\quad
n\geq 1,$ (45)

see [22, (5.15)], we get

$\displaystyle(\partial_{\tau}^{\alpha}(t\bar{\cal
E}_{u})^{n},v_{h})+(\nabla\cdot(t\bar{\cal
E}_{\boldsymbol{\sigma}})^{n},v_{h})$ $\displaystyle=$
$\displaystyle\alpha(\partial_{\tau}^{\alpha-1}\bar{\cal
E}_{u}^{{n-1}},v_{h}),$ (46) $\displaystyle(\kappa^{-1}(t\bar{\cal
E}_{\boldsymbol{\sigma}})^{n},{\bf w}_{h})-(t_{n}\bar{\cal
E}_{u}^{n},\nabla\cdot{\bf w}_{h})$ $\displaystyle=$
$\displaystyle(\beta(t_{n}{\cal E}_{u}^{n}),{\bf w}_{h}).$ (47)

Choose $v_{h}=t_{n}\bar{\cal E}_{u}^{n}$ and ${\bf w}_{h}=t_{n}\bar{\cal
E}_{\boldsymbol{\sigma}}^{n}$ to get

$(\partial_{\tau}^{\alpha}(t\bar{\cal E}_{u})^{n},(t\bar{\cal
E}_{u})^{n})+\|\sqrt{\kappa^{-1}}(t_{n}\bar{\cal
E}_{\boldsymbol{\sigma}}^{n})\|^{2}=(\beta(t_{n}{\cal
E}_{u}^{n}),t_{n}\bar{\cal
E}_{\boldsymbol{\sigma}}^{n})+\alpha(\partial_{\tau}^{\alpha-1}\bar{\cal
E}_{u}^{{n-1}},t_{n}\bar{\cal E}_{u}^{n}).$ (48)

Now, consider (37) at the time level $t_{n-1}$, apply $\partial_{\tau}^{-1}$
to both sides of the equation followed by choosing $v_{h}=t_{n}\bar{\cal
E}_{u}^{n}$ then choosing ${\bf w}_{h}=t_{n}\partial_{\tau}^{-1}\bar{\cal
E}_{\boldsymbol{\sigma}}^{{n-1}}$ in (38) and adding up yields

$(\partial_{\tau}^{\alpha-1}\bar{\cal E}_{u}^{{n-1}},t_{n}\bar{\cal
E}_{u}^{n})=-(\kappa^{-1}t_{n}\bar{\cal
E}_{\boldsymbol{\sigma}}^{n},\partial_{\tau}^{-1}\bar{\cal
E}_{\boldsymbol{\sigma}}^{{n-1}})+\alpha(\beta t_{n}{\cal
E}_{u}^{n},\partial_{\tau}^{-1}\bar{\cal E}_{\boldsymbol{\sigma}}^{{n-1}}).$

Substituting the new expression for $(\partial_{\tau}^{\alpha-1}\bar{\cal
E}_{u}^{n-1},t_{n}\bar{\cal E}_{u}^{n})$ in (48) and using the property (40),
we deduce that

$\partial_{\tau}^{\alpha}(\|(t\bar{\cal
E}_{u})^{n}\|^{2})+2\|\sqrt{\kappa^{-1}}(t_{n}\bar{\cal
E}_{\boldsymbol{\sigma}}^{n})\|^{2}\leq C\|t_{n}{\cal
E}_{u}^{n}\|^{2}+\|\sqrt{\kappa^{-1}}(t_{n}\bar{\cal
E}_{\boldsymbol{\sigma}}^{n})\|^{2}+C{\|\partial_{\tau}^{-1}\bar{\cal
E}_{\boldsymbol{\sigma}}^{n-1}\|^{2}}.$

Cancelling similar terms, applying $\partial_{\tau}^{-\alpha}$ and using the
estimate in Lemma 8,

$\|t_{n}\bar{\cal E}_{u}^{n}\|^{2}\leq
C\partial_{\tau}^{-\alpha}(\|(t\bar{\cal E}_{u})^{n}\|^{2})+C\max_{1\leq j\leq
n}\Big{(}\partial_{\tau}^{-\alpha}(\|(t\bar{\Phi}_{u})^{j}\|^{2})+t_{j-1}^{1-\alpha}\partial_{\tau}^{-1}(\|(\partial_{\tau}^{-\alpha}\bar{\Phi}_{u})\|^{2})^{j-1}\Big{)}.$

If the time step $\tau$ satisfies $Ca_{0}^{(-\alpha)}\tau^{\alpha}<1$, then an
application of a weakly singular discrete Gronwall inequality yields the
desired estimate in the lemma. ∎

We are ready to show optimal error estimates for the fully discrete solution.

######  Theorem 10.

Let $(u,{\boldsymbol{\sigma}})$ be the solutions of (9)-(10). Let
$(U_{h}^{n},\Sigma_{h}^{n})$ be the solution of (33)-(34) with
$u_{0h}=P_{h}u_{0}$. Then, for a small time step $\tau$, we have

$\|U_{h}^{n}-u(t_{n})\|\leq C(\tau
t_{n}^{\alpha\delta/2-1}+h^{2}t_{n}^{-\alpha(2-\delta)/2})\|u_{0}\|_{\dot{H}^{\delta}(\Omega)},\quad\delta\in[0,2].$
(49)

###### Proof.

Using the first estimate in Theorem 7 and the fact that
$\partial_{\tau}^{-\alpha}(t_{n}^{\beta})\leq Ct_{n}^{\alpha+\beta}$ for
$\beta>-1$, we find that

$\partial_{\tau}^{-\alpha}(\|(t\bar{\Phi}_{u})^{n}\|^{2})\leq
C\left(\tau^{2}t_{n}^{\alpha(\delta+1)}+h^{4}t_{n}^{-\alpha(1-\delta)+2}\right),$

and

$\partial_{\tau}^{-1}(\|\partial_{\tau}^{-\alpha}\bar{\Phi}_{u}\|^{2})^{n}\leq
C\left(\tau^{2}t_{n}^{\alpha(\delta+2)-1}+h^{4}t_{n}^{\alpha\delta+1}\right),$

for $\delta\in[0,2]$. Then, by Lemma 9, we see that

$\|t_{n}\bar{\cal E}_{u}^{n}\|^{2}\leq
C\left(\tau^{2}t_{n}^{\alpha(\delta+1)}+h^{4}t_{n}^{-\alpha(1-\delta)+2}\right).$
(50)

The desired estimate (49) follows now by using the triangle inequality and the
first estimate in Theorem 7. ∎

######  Remark 11.

To show the estimate for $\bar{\cal E}_{\boldsymbol{\sigma}}^{n}$, we operate
$\partial_{\tau}^{\alpha}$ on (47) so that

$(\kappa^{-1}\partial_{\tau}^{\alpha}(t\bar{\cal
E}_{\boldsymbol{\sigma}})^{n},{\bf
w}_{h})-(\partial_{\tau}^{\alpha}(t\bar{\cal E}_{u})^{n},\nabla\cdot{\bf
w}_{h})=(\beta\partial_{\tau}^{\alpha}(t{\cal E}_{u})^{n},{\bf w}_{h}).$

Next, we choose ${\bf w}_{h}=t_{n}\bar{\cal E}_{\boldsymbol{\sigma}}^{n}$ and
$v_{h}=\partial_{\tau}^{\alpha}(t\bar{\cal E}_{u})^{n}$ in (46) to arrive at

$\|\partial_{\tau}^{\alpha}(t\bar{\cal
E}_{u})^{n}\|^{2}+\frac{1}{2}\partial_{\tau}^{\alpha}\|\sqrt{\kappa^{-1}}t_{n}\bar{\cal
E}_{\boldsymbol{\sigma}}^{n}\|^{2}\\\ \leq
C\left(\|\partial_{\tau}^{\alpha}(t\bar{\Phi}_{u})^{n}\|^{2}+\|\partial_{\tau}^{\alpha-1}\bar{\cal
E}_{u}^{{n-1}}\|^{2}+\|\sqrt{\kappa^{-1}}t_{n}\bar{\cal
E}_{\boldsymbol{\sigma}}^{n}\|^{2}\right)+\frac{1}{2}\|\partial_{\tau}^{\alpha}(t\bar{\cal
E}_{u})^{n})\|^{2}.$

Canceling the similar terms, then applying $\partial_{\tau}^{-\alpha}$ to both
sides,

$\|\sqrt{\kappa^{-1}}t_{n}\bar{\cal E}_{\boldsymbol{\sigma}}^{n}\|^{2}\leq
C\max_{1\leq j\leq
n}\Big{(}\partial_{\tau}^{-\alpha}\Big{(}\|\partial_{\tau}^{\alpha}(t\bar{\Phi}_{u})\|^{2}\Big{)}^{j}+\partial_{\tau}^{-\alpha}\Big{(}\|\partial_{\tau}^{\alpha-1}\bar{\cal
E}_{u}\|^{2}\Big{)}^{j-1}\Big{)}\\\
+C\partial_{\tau}^{-\alpha}\left(\|\sqrt{\kappa^{-1}}(t\bar{\cal
E}_{\boldsymbol{\sigma}})^{n}\|^{2}\right)\,.$

Thus, if the time step $\tau$ is such that
$Ca_{0}^{(-\alpha)}\tau^{\alpha}<1$, then an application of a weakly singular
discrete Gronwall inequality yields

$t_{n}^{2}\|\sqrt{\kappa}\bar{\cal E}_{\boldsymbol{\sigma}}^{n}\|^{2}\leq
C\,\max_{1\leq j\leq
n}\Big{(}\partial_{\tau}^{-\alpha}\Big{(}\|\partial_{\tau}^{\alpha}(t\bar{\Phi}_{u})\|^{2}\Big{)}^{j}+\partial_{\tau}^{-\alpha}\Big{(}\|\partial_{\tau}^{\alpha-1}\bar{\cal
E}_{u}\|^{2}\Big{)}^{j-1}\Big{)},\quad n\geq 1.$

The second term $\|\partial_{\tau}^{\alpha-1}\bar{\cal E}_{u}^{j-1}\|$ on the
right-hand side can be readily bounded using the estimate (50). However,
deriving an appropriate bound for the first term
$\partial_{\tau}^{-\alpha}\left(\|\partial_{\tau}^{\alpha}(t\bar{\Phi}_{u})\|\right)^{j}$
which preserves the convergence rate seems to be a challenging task.

##  5 Numerical experiments

In this section, we present numerical tests to validate the theoretical
predictions discussed in the convergence theorems (Theorems 5 and 10) for
different choices of the fractional order $\alpha$. In the model problem (1),
we choose $T=0.5$, the spatial domain $\Omega$ is a unit square, $\kappa=1$,
and ${\bf F}(x,y)=(x,y)$. We consider three numerical tests with smooth and
nonsmooth initial data $u_{0}$.

The calculations are performed on a uniform (spatial) triangular grid. The
construction of the grid involves dividing the unit square into two triangles
by drawing a diagonal from (0, 0) to (1, 1). Then, a sequence of red
refinements is applied to generate a sequence of shape-regular grids. For the
mixed FEM, we perform the computation using the Raviart-Thomas finite elements
$RT0/P0$ and $RT1/P1^{\rm disc}$, where $P0$ and $P1^{\rm disc}$ denote the
sets of piecewise constant and linear functions, respectively. All
computations are performed using the finite element package written in Julia,
which provides the necessary tools for finite element analysis.

Example 1. We consider a test problem with the following initial data

$\displaystyle u_{0}(x,y)=x(1-x)y(1-y)\,.$

Since the exact solution is not available, we compute a reference solution
$u^{*}$ on the refinement level 6, representing a $2^{7}\times 2^{7}=128\times
128$ spatial grid. To examine the spatial discretization error, we choose a
(fixed) small $\tau$ so that the error incurred by the spatial discretization
is dominant. We compute the $L^{2}(\Omega)$-norm of the errors at the final
time $T$. The numerical results are given in Table 1 for the case of the
lowest-order Raviart-Thomas finite elements $RT_{0}/P_{0}$ for two different
values of $\alpha$; $\alpha=0.3$ and $\alpha=0.7$. The table shows a
convergence rate $O(h)$ for both the scalar variable $u$ and the flux variable
${\boldsymbol{\sigma}}$, which is consistent with our observation in Remark 6.
For the finite elements $RT_{1}/P_{1}^{\rm disc}$, the results are presented
in Table 2. We observe $O(h^{2})$-rates for both variables $u$ and
${\boldsymbol{\sigma}}$. As the initial data
$u_{0}\in\dot{H}^{5/2-\epsilon}(\Omega)$ for $\epsilon>0$, the numerical
results agree with the achieved theoretical estimates in Theorem 5.

| $\alpha=0.3$ | $\alpha=0.7$
---|---|---
level | $\|u^{*}-u_{h}\|$ | OC | $\|{\boldsymbol{\sigma}}^{*}-{\boldsymbol{\sigma}}_{h}\|$ | OC | $\|u^{*}-u_{h}\|$ | OC | $\|{\boldsymbol{\sigma}}^{*}-{\boldsymbol{\sigma}}_{h}\|$ | OC
1 | 1.074e-02 |  | 4.389e-02 |  | 1.724e-03 |  | 5.984e-03 |
2 | 4.315e-03 | 1.32 | 1.684e-02 | 1.38 | 8.772e-04 | 0.97 | 2.311e-03 | 1.37
3 | 2.144e-03 | 1.01 | 8.325e-03 | 1.07 | 4.388e-04 | 0.99 | 1.146e-03 | 1.01
4 | 1.044e-03 | 1.04 | 4.105e-03 | 1.02 | 2.230e-04 | 0.98 | 5.659e-04 | 1.02
5 | 5.299e-04 | 0.98 | 2.011e-03 | 1.03 | 1.101e-04 | 1.02 | 2.781e-04 | 1.02
Table 1: $L^{2}$-errors and spatial order of convergence (OC) for Example 1 using $RT_{0}/P_{0}$.  | $\alpha=0.3$ | $\alpha=0.7$
---|---|---
level | $\|u^{*}-u_{h}\|$ | OC | $\|{\boldsymbol{\sigma}}^{*}-{\boldsymbol{\sigma}}_{h}\|$ | OC | $\|u^{*}-u_{h}\|$ | OC | $\|{\boldsymbol{\sigma}}^{*}-{\boldsymbol{\sigma}}_{h}\|$ | OC
1 | 1.059e-03 |  | 4.906e-03 |  | 6.977e-04 |  | 5.437e-04 |
2 | 2.118e-04 | 2.32 | 1.222e-03 | 2.01 | 1.138e-04 | 2.62 | 1.085e-04 | 2.32
3 | 5.378e-05 | 1.98 | 2.298e-04 | 2.41 | 2.617e-05 | 2.12 | 2.765e-05 | 1.97
4 | 1.323e-05 | 2.02 | 5.705e-05 | 2.01 | 6.602e-06 | 1.99 | 6.828e-06 | 2.02
Table 2: $L^{2}$-errors and spatial order of convergence (OC) for Example 1
using $RT_{1}/P_{1}^{\rm dc}$ .

Example 2 (non-homogeneous). In this example, the initial data $u_{0}$ is a
two-dimensional continuous piecewise linear function (having a hat shape) such
that $u_{0}(1/2,1/2)=1/2$ and $u_{0}=0$ at the four corners of $\Omega.$ So,
$u_{0}\in\dot{H}^{3/2-\epsilon}(\Omega)$ for $\epsilon>0$. We choose a source
term $g$, that is the right-hand side of (1) is $g$ instead of $0$ in this
case, so that the exact solution is

$u(x,y,t)=16\sum_{m,n=0}^{\infty}(-1)^{mn}(\lambda_{m}\lambda_{n})^{-2}\sin(\lambda_{m}x)\sin(\lambda_{n}y)E_{\alpha,1}(-\lambda_{mn}t^{\alpha}),$

where $E_{\mu,\beta}(t):=\sum_{p=0}^{\infty}\frac{t^{p}}{\Gamma(\mu p+\beta)}$
is the Mittag-Leffler function, with parameters $\mu,\beta>0,$
$\lambda_{m}=(2m+1)\pi$, and $\lambda_{mn}=[(2m+1)^{2}+(2n+1)^{2}]\pi^{2}.$ In
this case,

$g(x,y,t)=16t^{\alpha-1}\sum_{m,n=0}^{\infty}(-1)^{mn}(\lambda_{m}\lambda_{n})^{-2}\varphi_{mn}(x,y)E_{\alpha,\alpha}(-\lambda_{mn}t^{\alpha}),$

where

$\varphi_{mn}(x,y)=\Big{(}\lambda_{m}x\cos(\lambda_{m}x)+2\sin(\lambda_{m}x)\Big{)}\sin(\lambda_{n}y)+\lambda_{n}y\,\sin(\lambda_{m}x)\cos(\lambda_{n}y)\,.$

To find the errors in space and time, we use the standard strategy which is to
consider the time step size small enough to find the convergence order in
space and vice versa. In order to keep the error contribution from time small,
we set the time step length $\tau=1/1200$. Table 3 shows the errors and
convergence rates in space for $RT_{1}/P_{1}^{\rm disc}$ elements. The
expected order of convergence $\mathcal{O}(h^{2})$ is obtained for the scalar
variable $u$ and the corresponding flux ${\boldsymbol{\sigma}}$. For
completeness, we present in Figure 1 the errors computed using the
$RT_{0}/P_{0}$ and $RT_{1}/P_{1}^{\rm disc}$ elements for a fixed $\alpha=0.5$
and different refinement levels.

To illustrate the convergence rate in time, we excluded the error in space.
The numerical tests are performed with $RT_{1}/P_{1}^{\rm disc}$ and a spatial
mesh consists of $128\times 128$ cells. The calculations were done with time
step lengths $\tau=0.5/N$ for different choices of $N$. In Table 4, the errors
and convergence orders are listed for the scalar variable $u$. A first-order
convergence rate is observed, which aligns with theoretical results in Theorem
10.

| $\alpha=0.3$ | $\alpha=0.7$
---|---|---
level | $\|u-u_{h}\|$ | OC | $\|{\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}_{h}\|$ | OC | $\|u-u_{h}\|$ | OC | $\|{\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}_{h}\|$ | OC
1 | 3.704e-03 |  | 1.565e-02 |  | 2.256e-03 |  | 8.893e-03 |
2 | 9.715e-04 | 1.931 | 3.865e-03 | 2.018 | 5.255e-04 | 2.101 | 2.192e-03 | 2.020
3 | 2.147e-04 | 2.178 | 9.667e-04 | 2.000 | 1.286e-04 | 2.031 | 4.853e-04 | 2.175
4 | 5.317e-05 | 2.013 | 2.394e-04 | 2.014 | 3.140e-05 | 2.033 | 1.194e-04 | 2.023
5 | 1.437e-05 | 1.887 | 6.651e-05 | 1.848 | 8.263e-06 | 1.926 | 3.061e-05 | 1.963
Table 3: $L^{2}$-errors and spatial order of convergence (OC) for Example 2
using $RT_{1}/P_{1}^{\rm dc}$.

Figure 1: $L^{2}$-errors for Example 2 using $RT_{0}/P_{0}$ (left) and $RT_{1}/P_{1}^{\rm disc}(right)$ for fixed $\alpha=0.5$. | $\alpha=0.3$ | $\alpha=0.7$
---|---|---
$N$ | $\|u-u_{h}\|$ | OC | $\|u-u_{h}\|$ | OC
10 | 3.1744041e-02 |  | 1.679786e-02 |
20 | 1.7003601e-02 | 0.901 | 7.812162e-03 | 1.105
40 | 8.0542173e-03 | 1.078 | 3.287313e-03 | 1.249
80 | 4.3481248e-03 | 0.889 | 1.595068e-03 | 1.043
160 | 2.1347824e-03 | 1.026 | 8.042011e-04 | 0.988
Table 4: $L^{2}$-errors and temporal order of convergence (OC) for Example 2.

Example 3 (non-homogeneous less smooth). We choose the initial data $u_{0}$ to
be a two-dimensional piecewise constant (discontinuous) function such that
$u_{0}(x,y)=1$ for $1/4\leq x,y\leq 3/4$, and zero elsewhere. So,
$u_{0}\in\dot{H}^{1/2-\varepsilon}(\Omega)$, for $\varepsilon>0$. We choose a
source term $g$ and boundary conditions such that the exact solution is

$u(x,y,t)=8\sum_{m,n=0}^{\infty}(-1)^{\lceil m/2\rceil\lceil
n/2\rceil}(\lambda_{m}\lambda_{n})^{-1}\sin(\lambda_{m}x)\sin(\lambda_{n}y)E_{\alpha,1}(-\lambda_{mn}t^{\alpha})\,.$

The source term in that case is

$g(x,y,t)=8t^{\alpha-1}\sum_{m,n=0}^{\infty}(-1)^{\lceil m/2\rceil\lceil
n/2\rceil}(\lambda_{m}\lambda_{n})^{-1}\varphi_{mn}(x,y)E_{\alpha,\alpha}(-\lambda_{mn}t^{\alpha}).$

| $\alpha=0.3$ | $\alpha=0.7$
---|---|---
level | $\|u-u_{h}\|$ | OC | $\|{\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}_{h}\|$ | OC | $\|u-u_{h}\|$ | OC | $\|{\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}_{h}\|$ | OC
1 | 5.587e-03 |  | 4.003e-02 |  | 2.493e-03 |  | 1.784e-02 |
2 | 1.311e-03 | 2.091 | 1.084e-02 | 1.885 | 6.549e-04 | 1.928 | 5.182e-03 | 1.783
3 | 3.019e-04 | 2.119 | 3.219e-03 | 1.752 | 1.570e-04 | 2.061 | 1.626e-03 | 1.672
4 | 7.956e-05 | 1.924 | 9.967e-04 | 1.691 | 3.932e-05 | 1.998 | 5.356e-04 | 1.603
Table 5: $L^{2}$-errors and temporal order of convergence (OC) for Example 3
using $RT_{1}/P_{1}^{\rm dc}$ .

This example studies the convergence order with respect to space. We fix the
time step length $\tau=6.25\times 10^{-4}$ and perform the simulations for the
finite element $RT_{1}/P_{1}^{\rm disc}$ on different refinement levels. The
errors and the convergence orders using $\alpha=0.3$ and $\alpha=0.7$ are
presented in Table 5. One can observe the predicted convergence order
$\mathcal{O}(h^{2})$ for the errors estimated in Theorems 5 and 10 for the
scalar variable $u$. By contrast, the convergence rate is of order
$\mathcal{O}(h^{1.6})$ for the flux variable ${\boldsymbol{\sigma}}$, which
may be seen as expected results. However, as the initial data
$u_{0}\in\dot{H}^{1/2-\varepsilon}(\Omega)$ for $\varepsilon>0$, the numerical
results show a positive effect on the convergence rate. Indeed, by (31)-(32),
and interpolation between $\dot{H}^{1}(\Omega)$ and $L^{2}(\Omega)$, it
follows that

$\|{\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}_{h}\|\leq
Ch^{3/2-\varepsilon}t^{-\alpha(2+\varepsilon)/2}\|u_{0}\|_{\dot{H}^{1/2-\varepsilon}(\Omega)},$

confirming our theoretical estimates.

## References

  * [1] A. A. Alikhanov, Boundary value problems for the diﬀusion equation of the variable order in diﬀerential and diﬀerence settings, Appl. Math. Comput., 219 (2012), 3938–3946.
  * [2] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Second edition, Springer-Verlag, New York, 1991.
  * [3] X. N. Cao, J.-L. Fu, and H. Huang, Numerical method for the time fractional Fokker–Planck equation, Adv. Appl. Math. Mech., 4 (2012), 848–863.
  * [4] S. Chen, F. Liu, P. Zhuang, and V. Anh, Finite difference approximations for the fractional Fokker–Planck equation, Appl. Math. Model., 33 (2009), 256–273.
  * [5] M. Cui, Compact exponential scheme for the time fractional convection-diffusion reaction equation with variable coefficients, J. Comput. Phys., 280 (2015), 143–163.
  * [6] W. Deng, Numerical algorithm for the time fractional Fokker–Planck equation, J. Comput. Phys., 227 (2007), 1510–1522.
  * [7] J. Dixon and S. McKee, Weakly singular Gronwall inequalities, ZAMM Z. Angew. Math.637 Mech., 66 (1986), 535–544.
  * [8] G. Fairweather, H. Zhang, X. Yang, and D. Xu, A backward Euler orthogonal spline collocation method for the time-fractional Fokker–Planck equation, Numer. Methods PDEs, 31 (2015), 1534–1550.
  * [9] D. Goswami and A. K. Pani, An alternate approach to optimal L2-error analysis of semidiscrete Galerkin methods for linear parabolic problems with nonsmooth initial data, Numer. Funct. Anal. Optim., 32 (2011), 946–982.
  * [10] B. I. Henry, T. A. M. Langlands, and P. Straka, Fractional Fokker–Planck equations for subdiffusion with space- and time-dependent forces, Phys. Rev. Lett., 105 (2010), 170–602.
  * [11] C. Huang, K. N. Le, and M. Stynes, A new analysis of a numerical method for the time-fractional Fokker–Planck equation with general forcing, IMA J. Numer. Anal., 40 (2020), 1217–1240.
  * [12] Y. Jiang, A new analysis of stability and convergence for finite difference schemes solving the time fractional Fokker–Planck equation, Appl. Math. Model., 39 (2015), 1163–1171.
  * [13] Y. Jiang and X. Xu, A monotone finite volume method for time fractional Fokker–Planck equations, Sci. China Math., 62 (2019), 783–794 (2019).
  * [14] S. Karaa, K. Mustapha, and A. K. Pani, Optimal error analysis of a FEM for fractional diffusion problems by energy arguments, J. Sci. Comput., 74 (2018), 519–-535.
  * [15] S. Karaa, Semidiscrete finite element analysis time fractional parabolic problems: a unified approach, SIAM J. Numer., 56 (2018), 1673–-1692.
  * [16] S. Karaa, Galerkin type methods for semilinear time-fractional diffusion problems, J. Sci. Comput., 83 (2020).
  * [17] S. Karaa and A. K. Pani, Mixed FEM for time-fractional diffusion problems with time-dependent coefficients, J. Sci. Comput., 83 (2020).
  * [18] A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo, Theory and Applications of Fractional Differential Equations, Elsevier, Amsterdam, 2006.
  * [19] N. K. Le, W. McLean and K. Mustapha, Numerical solution of the time-fractional Fokker–Planck equation with general forcing, SIAM J. Numer. Anal., 54 (2016), 1763–-1784.
  * [20] N. K. Le, W. McLean and K. Mustapha, A semidiscrete finite element approximation of a time-fractional Fokker–Planck equation with nonsmooth initial data, SIAM J. Sci. Comput., 40 (2018), A3831–-A3852.
  * [21] M. Luskin and R. Rannacher, On the smoothing property of the Galerkin method for parabolic equations, SIAM J. Numer. Anal. 19, 1982, 93–113.
  * [22] M. Al-Maskari and S. Karaa, The time-fractional Cahn–Hilliard equation: analysis and approximation, IMA J. Numer. Anal., 42, 1831–1865.
  * [23] W. McLean, Regularity of solutions to a time-fractional diffusion equation, ANZIAM J., 52 (2010), 123–138.
  * [24] W. McLean, K. Mustapha, R. Ali, and O. Knio, Regularity theory for time-fractional advection-diffusion-reaction equations, Comp. Math. Appl., 79 (2020), 947–961.
  * [25] W. McLean, K. Mustapha, R. Ali, and O. Knio, Well-posedness of time-fractional advection-diffusion-reaction equations, Fract. Calc. Appl. Anal., 22 (2019), 918–-944.
  * [26] K. Mustapha and W. McLean, Uniform stability for a spatially discrete, subdiffusive Fokker–Planck equation, Numer. Algo., 89 (2022), 1441–1463.
  * [27] R. Metzler, E. Barkai, and J. Klafter, Deriving fractional Fokker–Planck equations from a generalised master equation, Europhys. Lett., 46 (1999), 431–436.
  * [28] K. Mustapha, FEM for time-fractional diffusion equations, novel optimal error analyses, Math. Comp., 87 (2018), 2259–2272.
  * [29] K. Mustapha, O. M. Knio, O. P. Le Maître, A second–order accurate numerical scheme for a time–fractional Fokker–Planck equation, IMA J. Numer. Anal., 43 (2023), 2115–2136.
  * [30] P. Raviart and J. A. Thomas, Mixed finite element method for second order elliptic problems. In Mathematical Aspects of the Finite Element Method, Galligani I, Magenes E (eds). Lecture Notes in Mathematics, vol. 606. Springer: Berlin, 1977.
  * [31] A. Saadatmandi, M. Dehghan, and M.-R. Azizi, The sinc–Legendre collocation method for a class of fractional convection–diffusion equations with variable coefficients, Commun. Nonlinear Sci. Numer. Simul., 17 (2012), 4125–4136.
  * [32] V. Thomée, Galerkin finite element methods for parabolic problems, Springer-Verlag, Berlin, 2006.
  * [33] S. Vong and Z. Wang, A high order compact finite difference scheme for time fractional Fokker–Planck equations, Appl. Math. Lett., 43 (2015), 38–43.
  * [34] Y. Yang, Y. Huang, and Y. Zhou, Numerical solutions for solving time fractional Fokker–Planck equations based on spectral collocation methods. J. Comput. Math., 339 (2018), 389–404.
  * [35] Y. Zhao, P. Chen, W. Bu, X. Liu, and Y. Tang, Two mixed finite element methods for time-fractional diffusion equations, J. Sci. Comput., 70 (2017), 407-–428.
  * [36] D. Wang, A. Xiao, J. Zou, Long-time behavior of numerical solutions to nonlinear fractional ODEs, ESAIM Math. Model. Numer. Anal., 54 (2020), 335–-358.