# A posteriori error analysis for approximations of time-fractional
subdiffusion problems

Lehel Banjai  Maxwell Institute for Mathematical Sciences, School of
Mathematical & Computer Sciences, Heriot-Watt University, Edinburgh EH14 4AS,
UK<EMAIL_ADDRESS>and  Charalambos G. Makridakis  Institute for Applied
and Computational Mathematics-FORTH, Heraklion-Crete, GR 70013, Greece and
Department of Mathematics, University of Sussex, Brighton BN1 9QH, UK
<EMAIL_ADDRESS>
###### Abstract.

In this paper we consider a sub-diffusion problem where the fractional time
derivative is approximated either by the L1 scheme or by Convolution
Quadrature. We propose new interpretations of the numerical schemes which lead
to a posteriori error estimates. Our approach is based on appropriate
pointwise representations of the numerical schemes as perturbed evolution
equations and on stability estimates for the evolution equation. A posteriori
error estimates in $L^{2}(H)$ and $L^{\infty}(H)$ norms of optimal order are
derived. Extensive numerical experiments indicate the reliability and the
optimality of the estimators for the schemes considered, as well as their
efficiency as error indicators driving adaptive mesh selection locating
singularities of the problem.

###### 2020 Mathematics Subject Classification:

Primary 35R11, 65M06, 65M15

##  1\. Introduction

We consider the sub-diffusion problem: Find $u$ such that

(1.1) $\displaystyle D_{t}^{\beta}u(t)+Au(t)$ $\displaystyle=f(t),$
$\displaystyle 0<t<T,$ $\displaystyle u(0)$ $\displaystyle=u^{0},$

where $\beta\in(0,1)$, $f:\mathbb{R}_{\geq 0}\rightarrow H$ a given
inhomogeneity, and $u^{0}\in H$ the initial data. Here $A$ is a positive
definite, selfadjoint, linear operator on a Hilbert space $H$ with inner
product $\left\langle\cdot,\cdot\right\rangle$ and norm $\|\cdot\|$. Equations
(1.1) are understood as equalities in the Hilbert space $H$. We also denote
$|u|_{1}=\left\langle Au,u\right\rangle^{1/2}$ for $u\in V=D(A^{1/2})$ and
assume that $D(A)$ is dense in $H$. We further assume that there exists a
$\sigma>0$ such that $e^{-\sigma t}\|f(t)\|$ is bounded for $t\geq 0$. The
sub-diffusion equation has been used numerous times as a model in natural
sciences to describe anomalous diffusion processes, see e.g., [32, 14, 19, 31,
36]. As a simple example of the above setting we can take $H=L^{2}(\Omega)$,
$V=H_{0}^{1}(\Omega)$, and $A\cdot=-\Delta\cdot$ the Dirichlet Laplacian. Here
$\Omega$ denotes a bounded Lipshitz domain and $H_{0}^{1}(\Omega)$ the space
of $H^{1}(\Omega)$ functions with vanishing boundary trace. The setting also
applies to the corresponding symmetric conforming and non-conforming Galerkin
discretizations.

The fractional derivative $D_{t}^{\beta}$ is the Caputo fractional derivative
${}^{C}_{0}D_{t}^{\beta}$ of order $\beta\in(0,1)$ given by

(1.2)
$D_{t}^{\beta}g(t)=\frac{1}{\Gamma(1-\beta)}\int_{0}^{t}(t-\tau)^{-\beta}g^{\prime}(\tau)d\tau.$

Note that

$D_{t}^{\beta}g(t)=D_{t}^{\beta-1}g^{\prime}(t),$

where $D_{t}^{\beta-1}$ is the Riemann-Liouville integral of order $1-\beta$

(1.3)
$D_{t}^{\beta-1}h(t)=I_{t}^{1-\beta}h(t)=\frac{1}{\Gamma(1-\beta)}\int_{0}^{t}(t-\tau)^{-\beta}h(\tau)d\tau.$

A physical derivation of fractional derivatives via discrete random walks is
described in [12].

In the case $f\equiv 0$, semigroup techniques are used in [6, Theorem 4.1] to
show that the unique solution is analytic as a function of $t$ in an open
sector around the positive real axis and bounded there as

(1.4) $\|AD_{t}^{\beta}u\|\leq C|t|^{-\beta}{\|u^{0}\|}.$

For non-zero $f$ a representation formula for $u$ is given in (4.2); see Lemma
4.6.

Two popular discretization schemes for fractional derivatives and the sub-
diffusion equation (1.1) are convolution quadrature [24, 25] and the L1 scheme
[21, 40], both usually using a uniform time-step. The bound (1.4) indicates
that the solution is non-smooth at $t=0$. Indeed unless a correction term to
the standard schemes using uniform time-steps is added [16, 42] only low order
convergence can be expected even in the case of a smooth $f$. Furthermore, if
$f$ is not smooth for $t>0$, similar singularities can also occur elsewhere;
see Section 6.3.

It is important to be able to use non-uniform time-steps and to locate the
regions requiring refinement by the use of reliable a posteriori error
estimators. Such estimators are not available for fractional time dependent
problems. One of the reasons is the quite involved nature of time
dicretisations of fractional derivatives, and consequently the difficulty to
connect the numerical schemes to the exact evolution equation. In this paper
we are able to derive a posteriori error estimates for low order time discrete
schemes. In previous works for time-dependent partial differential equations
of parabolic or hyperbolic nature _reconstruction operators_ are introduced to
recover continuous objects from the approximate solutions, see the review [28]
and e.g., [33, 20, 1, 38, 27, 11, 5, 8, 18]. Then the derivation of the
estimates is reduced to (i) estimate of the reconstruction error and (ii) the
application of PDE stability estimates. This requires a _pointwise
representation_ of the numerical scheme as a perturbed evolution equation.
Appropriate forms of this kind are not obvious for fractional equations and
require new interpretations of the numerical schemes. In this paper we address
this problem and we derive new a posteriori error estimates for (1.1) valid
for a general time-discretization scheme. We apply it to both the L1 scheme
and low order convolution quadrature based on the backward Euler scheme.
Extensive numerical experiments indicate the reliability and the optimality of
the estimators for these schemes as well as their efficiency as error
indicators driving adaptive mesh selection locating singularities of the
problem.

The paper is organised as follows. L1 scheme and convolution quadrature are
described in Section 2 and Section 3. Both schemes are cast in a similar form
that allows direct comparison. The a posteriori estimates will rely on
stability estimates for the evolution problem (1.1) proved in Section 4, see
Theorem 4.2 and Theorem 4.6. We provide detailed proofs of both results
controlling the $L^{2}(H)$ and $L^{\infty}(H)$ norms of the solution
respectively. It is important for the subsequent a posteriori analysis to note
that both results allow the inclusion of a kernel in the forcing term. In
Section 5 we cast both L1 scheme and convolution quadrature methods in a
unified formulation and we derive corresponding appropriate pointwise forms
and error equations. Then we apply the stability estimates of Section 4 to
readily conclude the main a posteriori error estimates in $L^{2}(H)$ and
$L^{\infty}(H)$ norms, Theorem 5.1 and Theorem 5.2. Numerical results are
presented in Section 6 for all methods and estimators considered.

##  2\. Numerical method based on piecewise linear time-discretization

Let $0=t_{0}<t_{1}<\dots<t_{N}=T$ be a partition of $[0,T]$,
$I_{n}:=(t_{n},t_{n+1}]$ and $\kappa_{n}:=t_{n+1}-t_{n}$. To construct a
numerical scheme, given a sequence of values $U_{0},\dots,U_{N}$ we need to
define an approximation of the fractional derivative $D_{t}^{\beta}$ on this
grid. A standard approach is to apply the fractional derivative to a linear
interpolant of this data. This is a standard approach to constructing
numerical methods including Fredholm and Volterra integal equations [7, 13].
In the case of fractional derivatives this approach is often called the L1
scheme [21, 40]. We give the details next.

Denote by $\hat{U}(t)$ the piecewise linear interpolant defined by

(2.1) $\hat{U}(t_{n})=U_{n},\quad
n=0,\dots,N,\qquad\hat{U}|_{I_{n}}\in\mathbb{P}_{1}(I_{n}),$

where $\mathbb{P}_{1}$ is the space of linear functions. We also define a
projection operator $\Pi_{1}$ mapping continuous functions to piecewise linear
functions by interpolating in $t_{n}$. Namely

(2.2) $\Pi_{1}u(t_{n})=u(t_{n}),\quad
n=0,\dots,N,\qquad\Pi_{1}u|_{I_{n}}\in\mathbb{P}_{1}(I_{n}).$

Recalling that $D_{t}^{\beta}=D_{t}^{\beta-1}D_{t}$, we now define the
discrete fractional derivative by

(2.3)
$D_{t,\kappa}^{\beta}U(t):=D_{t}^{\beta}\hat{U}(t)=D_{t}^{\beta-1}\hat{U}^{\prime}(t).$

Note that for $t\in I_{n}$

$\begin{split}D_{t,\kappa}^{\beta}U(t)=&\frac{1}{\Gamma(1-\beta)}\sum_{j=0}^{n-1}\int_{t_{j}}^{t_{j+1}}(t-\tau)^{-\beta}\,d\tau\tfrac{U_{j+1}-U_{j}}{\kappa_{j}}\\\
&+\frac{1}{\Gamma(1-\beta)}\int_{t_{n}}^{t}(t-\tau)^{-\beta}\,d\tau\tfrac{U_{n+1}-U_{n}}{\kappa_{n}}\\\
=&\sum_{j=0}^{n}\omega_{j}(t)\tfrac{U_{j+1}-U_{j}}{\kappa_{j}},\end{split}$

where

(2.4)
$\begin{split}\omega_{j}(t)&=\frac{1}{\Gamma(1-\beta)}\int_{t_{j}}^{\min(t,t_{j+1})}(t-\tau)^{-\beta}\,d\tau\\\
&=\frac{1}{\Gamma(2-\beta)}\left((t-t_{j})^{1-\beta}-(t-\min(t,t_{j+1}))^{1-\beta}\right)\end{split}$

for $t\geq t_{j}$ and

(2.5) $\omega_{j}(t)=0\text{ for }t\leq t_{j}.$

This implies that $D_{t,\kappa}^{\beta}U(t)$ is continuous as

$\begin{split}D_{t,\kappa}^{\beta}U(t_{n}^{+})&=\sum_{j=0}^{n-1}\omega_{j}(t_{n})\tfrac{U_{j+1}-U_{j}}{\kappa_{j}}+\omega_{n}(t_{n})\tfrac{U_{n+1}-U_{n}}{\kappa_{n}}\\\
&=\sum_{j=0}^{n-1}\omega_{j}(t_{n})\tfrac{U_{j+1}-U_{j}}{\kappa_{j}}=D_{t,\kappa}^{\beta}U(t_{n}^{-}).\end{split}$

We denote the evaluation of the weights $\omega_{j}$ at $t_{n+1}$ by
$\omega_{n,j}$:

(2.6) $\begin{split}\omega_{n,j}&:=\omega_{j}(t_{n+1})\\\
&=\frac{1}{\Gamma(2-\beta)}\left((t_{n+1}-t_{j})^{1-\beta}-(t_{n+1}-t_{j+1})^{1-\beta}\right),\end{split}$

for $n\geq j$. By definition $\omega_{n,j}=0$ for $n<j$. With all the notation
introduced the discrete fractional derivative (2.3) evaluated at $t_{n+1}$ is
given by

(2.7)
$D_{t,\kappa}^{\beta}U(t_{n+1})=\sum_{j=0}^{n}\omega_{n,j}\frac{U_{j+1}-U_{j}}{\kappa_{j}}.$

For later it is also useful to define

$a_{j}(t)=\frac{1}{\Gamma(2-\beta)}(t-t_{j})^{1-\beta},\quad t\geq t_{j},$

and $a_{j}(t)=0$ for $t\leq t_{j}$. Using this definition we have for $t\leq
t_{n+1}$

$\displaystyle D_{t,\kappa}^{\beta}U(t)=$
$\displaystyle\sum_{j=0}^{n}\omega_{j}(t)\tfrac{U_{j+1}-U_{j}}{\kappa_{j}}$
$\displaystyle=$
$\displaystyle\sum_{j=0}^{n}(a_{j}(t)-a_{j+1}(t))\tfrac{U_{j+1}-U_{j}}{\kappa_{j}}$
(2.8) $\displaystyle=$ $\displaystyle
a_{0}(t)\tfrac{U_{1}-U_{0}}{\kappa_{0}}+\sum_{j=1}^{n}a_{j}(t)\left(\tfrac{U_{j+1}-U_{j}}{\kappa_{j}}-\tfrac{U_{j}-U_{j-1}}{\kappa_{j-1}}\right).$

We can now write down the fully discrete system: Find $U_{n+1}\in H$,
$n=0,1,\dots,N-1$ such that

(2.9) $D_{t,\kappa}^{\beta}U(t_{n+1})+AU_{n+1}=f_{n+1},$

where $f_{n+1}=f(t_{n+1})$ and $U_{0}=u^{0}$ or some approximation of the
initial data. Alternatively, recalling the definition of $\omega_{n,j}$ (2.6)
we can rewrite the system in a more familiar form as a finite difference
formula

(2.10)
$\sum_{j=0}^{n}\omega_{n,j}\frac{1}{\kappa_{j}}(U_{j+1}-U_{j})+AU_{n+1}=f_{n+1},\qquad
n=0,\dots,N-1.$

This can be seen as a classical collocation scheme for Volterra integral
equations [7, 13] and is equivalent to the L1 scheme of [40] and [21]. A fast
and memory efficient implementation of the solution of such a discretization
is developed in [22].

We summarize available a-priori results. Most of these results are for uniform
time-step $\kappa_{j}=\kappa$ for all $j$. Convergence order of
$O(\kappa^{2-\beta})$ is proved in [40] under the assumption that $u\in
C^{2}[0,T]$. In [15] the authors argue that in general even for smooth data
this smoothness of the solution does not hold and only linear $O(\kappa)$
convergence order is obtained. Namely, they prove that for $f\equiv 0$

$\|u(t_{n})-U_{n}\|\leq C\kappa t_{n}^{-1}\|u^{0}\|,\qquad n\geq 1$

and

$\|u(t_{n})-U_{n}\|\leq C\kappa t_{n}^{\beta-1}\|Au^{0}\|,\qquad n\geq 1,$

if $u^{0}\in V$. The main reason for this is a singularity at $t=0$, as even
for $f\equiv 0$, in general $\|D_{t}^{\beta}u(t)\|\leq Ct^{-\beta}\|u^{0}\|$
[37]. This suggests that grading towards $t=0$ would be advantageous. Indeed
in [39] it is proved that the optimal convergence is recovered when using a
graded mesh, where the operator $A$ is the differential operator
$Au=-\frac{\partial^{2}}{\partial x^{2}}u+c(x)u$ in one spatial dimension. The
authors prove that for optimal convergence uniformly for $t>0$ in the case of
a smooth $f$, it is necessary to choose the graded mesh

(2.11) $t_{j}=T(j/N)^{k}\text{ with }k\geq\frac{2-\beta}{\beta},\quad
j=0,\dots,N.$

A modified L1 scheme described in [42] recovers the $O(\kappa^{2-\beta})$
convergence order for sufficiently smooth $f$

$\|u(t_{n})-U_{n}\|\leq Ct_{n}^{\beta-2}\kappa^{2-\beta},$

with $C$ a constant depending on $f$ and $u_{0}$. The modified scheme reads

(2.12)
$\begin{split}D_{t,\kappa}^{\beta}U(t_{1})+AU_{1}&=f(t_{1})-\frac{1}{2}(Au_{0}-f(0)),\\\
D_{t,\kappa}^{\beta}U(t_{n+1})+AU_{n+1}&=f_{n+1},\qquad
n=1,\dots,N-1.\end{split}$

##  3\. Convolution quadrature

Another popular discretization method for fractional derivatives is
convolution quadrature [24, 25] with non-uniform time-step schemes
investigated in [23]. The low order convolution quadrature (CQ) based on
backward Euler discretization can be given in the following form

(3.1) $[^{\text{\tiny
CQ}}D_{t,\kappa}^{\beta}V]_{n+1}:=\sum_{j=0}^{n}\omega^{\text{\tiny
CQ}}_{n,j}\frac{V_{j+1}-V_{j}}{\kappa_{j}},$

where $\omega^{\text{\tiny CQ}}_{n,j}$ are convolution weights for the
fractional integral of order $1-\beta$ and is given by

(3.2) $\omega^{\text{\tiny CQ}}_{n,j}=\kappa_{j}\frac{1}{2\pi
i}\int_{\sigma-i\infty}^{\sigma+i\infty}z^{\beta-1}\prod_{k=j}^{n}\frac{1}{1-\kappa_{k}z}dz,$

for a fixed $\sigma\in(0,\min_{k}\kappa_{k}^{-1})$; due to the analyticity of
the integrand, the value of $\omega^{\text{\tiny CQ}}_{n,j}$ is independent of
$\sigma$. The expression (3.1) is of the same form as the L1 discrete
derivative (2.7). As this is not the standard way to present convolution
quadrature we give a detailed derivation of the scheme in the appendix.

######  Remark 3.1.

In order to understand better the formula (3.2) it is of interest to compare
$\omega^{\text{\tiny CQ}}_{n,j}$ with the weights $\omega_{n,j}$ of the L1
scheme. Using the approximation $(1-z)^{-1}=e^{z}+O(z^{2})$ we have

$\begin{split}\omega^{\text{\tiny CQ}}_{n,j}&=\kappa_{j}\frac{1}{2\pi
i}\int_{\sigma-i\infty}^{\sigma+i\infty}z^{\beta-1}\prod_{k=j}^{n}\frac{1}{1-\kappa_{k}z}dz\\\
&\approx\kappa_{j}\frac{1}{2\pi
i}\int_{\sigma-i\infty}^{\sigma+i\infty}z^{\beta-1}e^{\sum_{k=j}^{n}\kappa_{k}z}dz\\\
&=\kappa_{j}\frac{1}{2\pi
i}\int_{\sigma-i\infty}^{\sigma+i\infty}z^{\beta-1}e^{(t_{n+1}-t_{j+1})z}dz\\\
&=\kappa_{j}\frac{1}{\Gamma(1-\beta)}(t_{n+1}-t_{j})^{-\beta},\end{split}$

where in the last step we used that the inverse Laplace transform of
$z^{\beta-1}$ is $\frac{1}{\Gamma(1-\beta)}t^{-\beta}$. A rigorous proof in a
more general setting and with error estimates is given in [25, Theorem 4.1]
for the case of uniform time-steps and in [4, Proposition 2.2] for non-uniform
time-steps.

Returning to the L1 discrete derivative (2.6)

$\begin{split}\omega_{n,j}&=\frac{1}{\Gamma(2-\beta)}\left((t_{n+1}-t_{j})^{1-\beta}-(t_{n+1}-t_{j+1})^{1-\beta}\right)\\\
&=\frac{1}{\Gamma(2-\beta)}\left((t_{n+1}-t_{j+1}+\kappa_{j})^{1-\beta}-(t_{n+1}-t_{j+1})^{1-\beta}\right)\\\
&=\kappa_{j}\frac{1}{\Gamma(1-\beta)}(t_{n+1}-t_{j+1})^{-\beta}+O((t_{n+1}-t_{j+1})^{-\beta-1}\kappa_{j}^{2}),\end{split}$

for $t_{n+1}>t_{j+1}$. Hence, the weights $\omega_{n,j}$ and
$\omega^{\text{\tiny CQ}}_{n,j}$ have a similar behaviour for $n\gg j$.

To simplify the computation of $\omega^{\text{\tiny CQ}}_{n,j}$ we can
transform the integration contour to the negative real axis

(3.3) $\begin{split}\omega^{\text{\tiny CQ}}_{n,j}&=\kappa_{j}\frac{1}{2\pi
i}\left[\int_{-\infty}^{0}e^{-(\beta-1)\pi
i}x^{\beta-1}\prod_{k=j}^{n}\frac{1}{1+x\kappa_{k}}dx+\int_{0}^{\infty}e^{(\beta-1)\pi
i}x^{\beta-1}\prod_{k=j}^{n}\frac{1}{1+x\kappa_{k}}dx\right]\\\
&=\kappa_{j}\frac{\sin((1-\beta)\pi)}{\pi}\int_{0}^{\infty}x^{\beta-1}\prod_{k=j}^{n}\frac{1}{1+x\kappa_{k}}dx.\end{split}$

Note that

$\omega^{\text{\tiny
CQ}}_{n,n}=\kappa_{n}^{1-\beta}\frac{\sin((1-\beta)\pi)}{\pi}\int_{0}^{\infty}x^{\beta-1}\frac{1}{1+x}dx=\kappa_{n}^{1-\beta}.$

For $j<n$ we truncate the integral and then apply an $hp$-quadrature as done
for the uniform time-step version in [3]. This approach to computing the
weights can be used to develop a computationally and memory efficient
algorithm [3].

With this definition of the discrete derivative, the time-discrete system is
given by

(3.4) $[^{\text{\tiny
CQ}}D_{t,\kappa}^{\beta}U]_{n+1}+AU_{n+1}=f(t_{n+1}),\qquad n=0,\dots,N-1,$

or written as

(3.5) $\sum_{j=0}^{n}\omega^{\text{\tiny
CQ}}_{n,j}\frac{U_{j+1}-U_{j}}{\kappa_{j}}+AU_{n+1}=f(t_{n+1}),\qquad
n=0,\dots,N-1.$

For non-uniform steps, convolution quadrature has been analyzed in [23]. This
analysis is however for hyperbolic problems with assumptions that are too
stringent for the parabolic problem investigated here. For uniform steps [26,
Theorem 2.2] and the estimate
$\left\|\left(A+\lambda^{\beta}I\right)^{-1}\right\|\leq C|\lambda|^{-\beta}$,
see (4.1), imply that

$\|U_{n}-u(t_{n})\|\leq Ct_{n}^{\beta-1}\kappa,\qquad n=1,\dots,N=T/\kappa.$

This implies that

$\displaystyle\left(\kappa\sum_{n=1}^{N}\|U_{n}-u(t_{n})\|^{2}\right)^{1/2}$
$\displaystyle\leq
C\kappa\left(\kappa\sum_{n=1}^{N}t_{n}^{2\beta-2}\right)^{1/2}$
$\displaystyle\leq C{\kappa\max(T^{\beta-1/2},\kappa^{\beta-1/2})}$ (3.6)
$\displaystyle\leq C(T)\kappa^{\min(1,\beta+\frac{1}{2})},$

where in the penultimate step we used an upper bound derived from the
equivalent continuous integral; see [2, Lemma 4.1]. Therefore for $t_{n}$ away
from the origin linear order convergence is obtained, whereas in the $L^{2}$
norm we expect to have suboptimal convergence order if $\beta<1/2$ unless a
graded discretization is used.

######  Remark 3.2.

We notice an important difference between the L1 scheme and CQ. The optimal
convergence for the L1 scheme is $O(\kappa^{2-\beta})$, whereas for CQ the
optimal convergence order is linear $O(\kappa)$. For non-uniform time-step
schemes we would expect the same to hold with $\kappa=\max_{j}\kappa_{j}$.

##  4\. Stability of the continuous problem

In this section we derive stability estimates for the evolution problem (1.1)
in $L^{2}(H)$ and $L^{\infty}(H),$ Theorem 4.2 and Theorem 4.6. These bounds
will be then instrumental to derive the a posteriori estimates in the next
section. We provide detailed proofs of both results aiming to include a kernel
in the forcing term. This is crucial in order to obtain an optimal estimator
which can recognise the $O(\kappa^{2-\beta})$, i.e., higher than linear,
convergence order achievable by the L1 scheme.

Important for the analysis in $L^{2}(H)$ will be the following positivity
result. It can be deduced from Lemma 1.7.2 in [41]; see also Lemma 3.1 in
[34]. A similar result to the second inequality below can also be found in
[30, Theorem A.1] but with a slightly less optimal constant. We nevertheless
believe that the following proof is of interest. It is similar to the proof
for the discrete stability given in [40].

######  Lemma 4.1.

Let $g\in C^{1}([0,T];H)$ and $\beta\in(0,1)$, where $H$ is a Hilbert space
with inner product $\langle\cdot,\cdot\rangle$ and norm $\|\cdot\|$. Then

$\begin{split}\int_{0}^{T}\langle
D_{t}^{\beta}g(t),g(t)\rangle\,dt\geq&\frac{1/2}{\Gamma(1-\beta)}\int_{0}^{T}\left((T-t)^{-\beta}+t^{-\beta}\right)\|g(t)\|^{2}\,dt\\\
&\qquad-\frac{1}{\Gamma(1-\beta)}\int_{0}^{T}t^{-\beta}\langle
g(0),g(t)\rangle\,dt\\\
\geq&\frac{(T/2)^{-\beta}}{\Gamma(1-\beta)}\int_{0}^{T}\|g(t)\|^{2}\,dt-\frac{1}{\Gamma(1-\beta)}\int_{0}^{T}t^{-\beta}\langle
g(0),g(t)\rangle\,dt.\end{split}$

###### Proof.

To simplify notation, we set $H=\mathbb{R}$. Let $\varepsilon>0$ and consider

$\begin{split}\int_{0}^{T}\int_{0}^{t}(t-\tau+\varepsilon)^{-\beta}g^{\prime}(\tau)g(t)d\tau
dt=&-\beta\int_{0}^{T}\int_{0}^{t}(t-\tau+\varepsilon)^{-\beta-1}g(\tau)g(t)d\tau
dt\\\
&+\varepsilon^{-\beta}\int_{0}^{T}g^{2}(t)dt-g(0)\int_{0}^{T}(t+\varepsilon)^{-\beta}g(t)dt.\end{split}$

Next we bound the first term on the right-hand side as follows

$\begin{split}\beta\int_{0}^{T}\int_{0}^{t}(t-\tau+\varepsilon)^{-\beta-1}g(\tau)g(t)d\tau
dt\leq&\frac{1}{2}\beta\int_{0}^{T}\int_{0}^{t}(t-\tau+\varepsilon)^{-\beta-1}g^{2}(\tau)d\tau
dt\\\
&+\frac{1}{2}\beta\int_{0}^{T}\int_{0}^{t}(t-\tau+\varepsilon)^{-\beta-1}g^{2}(t)d\tau
dt\\\
=&\frac{1}{2}\beta\int_{0}^{T}\int_{\tau}^{T}(t-\tau+\varepsilon)^{-\beta-1}dt\,g^{2}(\tau)d\tau\\\
&+\frac{1}{2}\int_{0}^{T}(\varepsilon^{-\beta}-(t+\varepsilon)^{-\beta})g^{2}(t)dt\\\
=&-\frac{1}{2}\int_{0}^{T}((T-\tau+\varepsilon)^{-\beta}-\varepsilon^{-\beta})g^{2}(\tau)d\tau\\\
&+\frac{1}{2}\int_{0}^{T}(\varepsilon^{-\beta}-(t+\varepsilon)^{-\beta})g^{2}(t)dt.\end{split}$

Returning to the original calculation, this inequality implies

$\begin{split}\int_{0}^{T}\int_{0}^{t}(t-\tau+\varepsilon)^{-\beta}g^{\prime}(\tau)g(t)d\tau
dt\geq&\frac{1}{2}\int_{0}^{T}(T-\tau+\varepsilon)^{-\beta}g^{2}(\tau)d\tau+\frac{1}{2}\int_{0}^{T}(t+\varepsilon)^{-\beta}g^{2}(t)dt\\\
&-g(0)\int_{0}^{T}(t+\varepsilon)^{-\beta}g(t)dt.\end{split}$

Taking the limit $\varepsilon\rightarrow 0^{+}$ and including the constant
$\frac{1}{\Gamma(1-\beta)}$ gives the first inequality. Finding that
$\min_{\tau\in[0,T]}\frac{1}{2}((T-\tau)^{-\beta}+\tau^{-\beta})=(T/2)^{-\beta}$
gives the second. ∎

######  Theorem 4.2.

Let $u$ be the solution of (1.1). Then

$\begin{split}2^{-1}\int_{0}^{T}g_{\beta,T}(t)&\|u(t)\|^{2}\,dt+\int_{0}^{T}|u(t)|^{2}_{1}dt\\\
&\leq\int_{0}^{T}\langle
f(t),u(t)\rangle\,dt+\frac{1}{\Gamma(1-\beta)}\int_{0}^{T}t^{-\beta}\langle
u^{0},u(t)\rangle\,dt,\end{split}$

where

$g_{\beta,T}(t)=\frac{1}{\Gamma(1-\beta)}\left((T-t)^{-\beta}+t^{-\beta}\right).$

###### Proof.

Testing (1.1) with $u$ and using Lemma 4.1 gives the result. ∎

In the following corollary we estimate the forcing term by including a kernel
dictated by the appearance of $g_{\beta,T}$ in the lower bound of the estimate
in Theorem 4.2. This is in contrast to estimators of the forcing term in the
dual of the $V$ norm typically appearing in diffusion problems. The next bound
will be instrumental in the a posteriori analysis of the next section.

######  Corollary 4.3.

Let $u$ be a solution of (1.1). Then

$\begin{split}\int_{0}^{T}\|u(t)\|^{2}+|u(t)|_{1}^{2}\,dt\leq&C^{1}_{T,\beta}\int_{0}^{T}(g_{\beta,T}(t))^{-1}\|f(t)\|^{2}\,dt\\\
&+C^{2}_{T,\beta}\|u^{0}\|^{2}\end{split}$

where $g_{\beta,T}$ as in Theorem 4.2,

$C^{1}_{T,\beta}=2\max(2^{1-\beta}\Gamma(1-\beta)T^{\beta},1)$

and

$C^{2}_{T,\beta}=2\frac{T^{1-\beta}}{\Gamma(2-\beta)}\max(2^{1-\beta}\Gamma(1-\beta)T^{\beta},1).$

###### Proof.

Using the Cauchy-Schwarz and Young’s inequalities in Theorem 4.2 gives

$\begin{split}2^{-1}(1-\epsilon_{1}-\epsilon_{2})\int_{0}^{T}g_{\beta,T}(t)\|u(t)\|^{2}\,dt&+\int_{0}^{T}|u(t)|_{1}\,dt\\\
\leq&\frac{1}{2}\epsilon_{1}^{-1}\int_{0}^{T}(g_{\beta,T}(t))^{-1}\|f(t)\|^{2}\,dt\\\
&+\frac{1}{2}\epsilon_{2}^{-1}\frac{1}{\Gamma(1-\beta)}\int_{0}^{T}\frac{1}{t^{2\beta}\left((T-t)^{-\beta}+t^{-\beta}\right)}dt\|u^{0}\|^{2}\\\
\leq&\frac{1}{2}\epsilon_{1}^{-1}\int_{0}^{T}(g_{\beta,T}(t))^{-1}\|f(t)\|^{2}\,dt+\frac{1}{2}\epsilon_{2}^{-1}\frac{T^{1-\beta}}{\Gamma(2-\beta)}\|u^{0}\|^{2}.\end{split}$

Using
$g_{\beta,T}(t)=\frac{1}{\Gamma(1-\beta)}\left((T-t)^{-\beta}+t^{-\beta}\right)\geq\frac{2}{\Gamma(1-\beta)}(T/2)^{-\beta}$
we have

$\begin{split}\frac{(1-\epsilon_{1}-\epsilon_{2})(T/2)^{-\beta}}{\Gamma(1-\beta)}\int_{0}^{T}\|u(t)\|^{2}\,dt+\int_{0}^{T}|u(t)|_{1}\,dt\leq&\frac{1}{2}\epsilon_{1}^{-1}\int_{0}^{T}(g_{\beta,T}(t))^{-1}\|f(t)\|^{2}\,dt\\\
&+\frac{1}{2}\epsilon_{2}^{-1}\frac{T^{1-\beta}}{\Gamma(2-\beta)}\|u^{0}\|^{2}.\end{split}$

Setting $\epsilon_{1}=\epsilon_{2}=1/4$ gives the result. ∎

Our next task is the $L^{\infty}(H)$ stability and the proof of Theorem 4.6
below. Towards this goal, an Abelian-Tauberian theorem with the names of Hardy
and Littlewood, and Karamata associated with it [10, p. 445], connects the
asymptotic behaviour of a function $k(t)$ at $t=0$ with the behaviour of its
transform $K(s)$ at infinity. The more modest aim here is to give a bound on
$k$ with an explicit constant.

######  Lemma 4.4.

Assume that $K(s)$ is an analytic function such that

$|K(s)|\leq C_{K}|s|^{-\mu}\qquad|\arg s|<\pi-\varphi,$

for some $\mu>0$ and $\varphi\in[0,\pi/2)$. Then $k$, the inverse Laplace
transform of $K$, is analytic for $t>0$ and

$|k(t)|\leq C_{K}C_{t,\mu,\varphi}t^{\mu-1},$

where

$\begin{split}C_{t,\mu,\varphi}=\min_{r_{0}\geq
0}\left(e^{r_{0}t}(r_{0}t)^{1-\mu}+\frac{1}{\pi}(\cos\varphi)^{\mu-1}\Gamma(1-\mu,r_{0}\cos\varphi)\right).\end{split}$

For $\mu\in(0,1)$, choosing $r_{0}=0$, we have the explicit bound

$C_{t,\mu,\varphi}\leq\frac{1}{\pi}(\cos\varphi)^{\mu-1}\Gamma(1-\mu),$

whereas for $\mu=1$, choosing $r_{0}=t^{-1}$,

$C_{t,1,\varphi}\leq e+\frac{1}{\pi}E_{1}(t^{-1}\cos\varphi).$

###### Proof.

We begin by using the inverse Laplace transform

$k(t)=\frac{1}{2\pi i}\int_{\Gamma}e^{st}K(s)ds$

to represent $k$. Here $\Gamma$ consists of the two arcs $\Gamma_{-}\colon
s=re^{(-\pi+\varphi)i}$ and $\Gamma_{+}\colon s=re^{(\pi-\varphi)i}$, for
$r\in[r_{0},\infty)$ and the circular contour connecting these:
$\Gamma_{0}=\\{r_{0}e^{i\varphi}:\varphi\in(-\pi+\varphi,\pi-\varphi)\\}$
where $r_{0}\in(0,1]$. Starting with $\Gamma_{0}$ we have

$\left|\frac{1}{2\pi i}\int_{\Gamma_{0}}e^{st}K(s)ds\right|\leq
C_{K}e^{r_{0}t}r_{0}^{1-\mu}.$

The bound on $\Gamma_{+}\cup\Gamma_{-}$ can be done at the same time

$\begin{split}\left|\frac{1}{2\pi
i}\int_{\Gamma_{-}\cup\Gamma_{+}}e^{st}K(s)ds\right|&\leq\frac{C_{K}}{\pi}\int_{r_{0}}^{\infty}r^{-\mu}e^{-rt\cos\varphi}dr\\\
&=\frac{C_{K}}{\pi}t^{\mu-1}(\cos\varphi)^{\mu-1}\int_{r_{0}\cos\varphi}^{\infty}r^{-\mu}e^{-r}dr\\\
&=\frac{C_{K}}{\pi}t^{\mu-1}(\cos\varphi)^{\mu-1}\Gamma(1-\mu,r_{0}\cos\varphi),\end{split}$

where $\Gamma(\cdot,\cdot)$ is the incomplete Gamma function. For
$\mu\in(0,1)$ we set $r_{0}=0$

$|k(t)|\leq\frac{C_{K}}{\pi}t^{\mu-1}(\cos\varphi)^{\mu-1}\Gamma(1-\mu,0)=\frac{C_{K}}{\pi}t^{\mu-1}(\cos\varphi)^{\mu-1}\Gamma(1-\mu).$

For $\mu=1$, we choose $r_{0}=t^{-1}$. Using $\Gamma(0,r)=E_{1}(r)$ [35,
(6.11.1)], where $E_{1}(\cdot)$ is the exponential integral [35, (6.2.1)], we
obtain that

$|k(t)|\leq C_{K}(e+\frac{1}{\pi}E_{1}(t^{-1}\cos\varphi)).$

∎

######  Remark 4.5.

Note that $E_{1}(t^{-1}\cos\varphi)$ grows as $O(\log t^{-1})$ for
$t\rightarrow\infty$.

Combining the above lemma with a resolvent bound for $A$ will give another
stability bound. As $A$ is a positive definite, self-adjoint operator its
spectrum lies on the positive real axis, and hence from [17, V (3.16)] we have
the resolvent bound

(4.1) $\|(A+\lambda
I)^{-1}\|{=\frac{1}{\operatorname{dist}(\lambda,\sigma(A))}}\leq\frac{1}{\sin\theta|\lambda|}\quad\text{
for }|\arg\lambda|<\pi-\theta$

and any fixed $\theta\in(0,\pi/2)$.

######  Theorem 4.6.

Let $u$ be the solution of (1.1) with the additional assumption $u^{0}\in
D(A)$. Then the bound

$\|u(t)\|\leq\|u^{0}\|+\frac{1}{\sin\theta}C_{\beta,\varphi}\int_{0}^{t}(t-\tau)^{\beta-1}\|f(\tau)-Au^{0}\|d\tau$

holds for $t>0$, where $\varphi=\max(0,\pi-(\pi-\theta)/\beta)$ and

$C_{\beta,\varphi}=\frac{1}{\pi}(\cos\varphi)^{\beta-1}\Gamma(1-\beta).$

###### Proof.

Writing $w=u-u^{0}$ (1.1) becomes

$D_{t}^{\beta}w+Aw=f-Au^{0}.$

After taking the Laplace transform we obtain that

$\hat{w}(s)=(s^{\beta}{+}A)^{-1}\left(\hat{f}(s)-\frac{1}{s}Au_{0}\right),$

where $\hat{w}$ and $\hat{f}$ are the Laplace transforms of $w$ and $f$; note
that the growth condition on $f$ ensures the existence of the Laplace
transform.

Thus, we have the representation

(4.2) $u(t)=u^{0}+\int_{0}^{t}k_{1}(t-\tau)(f(\tau)-Au^{0})d\tau,$

with

$K_{1}(s)=(s^{\beta}{+}A)^{-1}.$

The resolvent bound (4.1) implies

$\|K_{1}(s)\|\leq\frac{1}{\sin\theta|s|^{\beta}},$

for $|\arg s^{\beta}|<\pi-\theta$ or equivalently for $|\arg
s|<(\pi-\theta)/\beta$. Combining this with Lemma 4.4 and
$\varphi=\max(0,\pi-(\pi-\theta)/\beta)$ gives that

$\|u(t)\|\leq\|u^{0}\|+\frac{1}{\sin\theta}C_{\beta,\varphi}\int_{0}^{t}(t-\tau)^{\beta-1}\|f(\tau)-Au^{0}\|d\tau.$

∎

##  5\. A posteriori error analysis: piecewise linear reconstruction

Recall that $\hat{U}$ denotes the continuous piecewise linear interpolant of
the data $U_{0},U_{1},\dots,U_{N}$. Due to (2.9), (3.4) both schemes can be
written in the form

(5.1)
$Q^{\beta}\hat{U}(t_{n+1})+AU_{n+1}=Q^{\beta}\hat{U}(t_{n+1})+A\hat{U}(t_{n+1})=f(t_{n+1}),$

where for $n=0,\dots,N-1,$

(5.2)
$Q^{\beta}\hat{U}(t_{n+1})=\sum_{j=0}^{n}\int_{I_{j}}K_{n,j}\hat{U}^{\prime}(\tau)d\tau\quad\text{where
}K_{n,j}=\begin{cases}\omega_{n,j}/{\kappa_{j}}&\quad\text{for L1 scheme},\\\
\omega^{\text{\tiny CQ}}_{n,j}/{\kappa_{j}}&\quad\text{for C-Q
scheme.}\end{cases}$

Furthermore, for consistency reasons we set $Q^{\beta}\hat{U}(0)=0.$ We need a
piecewise equation to be valid for all $t.$ In contrast to the parabolic
problems where $\hat{U}^{\prime}(t_{n+1})=\hat{U}^{\prime}(t),$ for all $t\in
I_{n},$ and thus the pointwise equation can be extended for all $t,$ in our
case we should proceed in a different way. In fact, we derive a pointwise
equation for $\hat{U}$ by applying piecewise linear interpolation to (5.1).
Then we conclude,

(5.3) $\Pi_{1}Q^{\beta}\hat{U}(t)+A\hat{U}(t)=\Pi_{1}f(t)+\hat{G}(t),\qquad
t\in[0,T],$

with initial data $\hat{U}(0)=U_{0}$, where $\hat{G}$ is a piecewise linear
correction term whose role we explain next. Here, $\Pi_{1}Q^{\beta}\hat{U}$ is
the piecewise linear interpolant of the function $Q^{\beta}\hat{U}$ defined in
(5.2). By definition, (5.3) is satisfied for $t=t_{j}$, $j=1,\dots,N$ if
$\hat{G}(t_{j})=0$. However, as we have set $Q^{\beta}\hat{U}(0)=0$, in order
that the equation is satisfied also in the first interval we need that
$\hat{G}(0)=AU_{0}-f(0)$. Hence, the piecewise linear correction function
$\hat{G}$ is defined by interpolating

(5.4) $\hat{G}(0)=AU_{0}-f(0),\quad\hat{G}(t_{j})=0,\;j=1,\dots,N.$

We conclude therefore that $\hat{U}$ is the solution of the original evolution
problem (1.1) with a modified right-hand side

(5.5)
$D_{t}^{\beta}\hat{U}(t)+A\hat{U}(t)=\Pi_{1}f(t)+\hat{G}(t)+D_{t}^{\beta}\hat{U}(t)-\Pi_{1}Q^{\beta}\hat{U}(t),\qquad
t\in(0,T].$

###  5.1. Error equation

The error $e=u-\hat{U}$ is then the solution of

$D_{t}^{\beta}e+Ae=f-D_{t}^{\beta}\hat{U}-A\hat{U}=:\hat{R},\qquad
e(0)=u^{0}-U_{0},$

where in view of (5.5) the residual $\hat{R}$ is given by

(5.6)
$\hat{R}=f(t)-\Pi_{1}f(t)-\hat{G}(t)-(D_{t}^{\beta}\hat{U}(t)-\Pi_{1}Q^{\beta}\hat{U}(t)),\qquad
t\in(0,T].$

The continuous stability of the fractional problem (1.1) implies the desired a
posteriori bounds. The stability estimates in Corollary 4.3 and Theorem 4.6
imply two different a posteriori error estimates.

######  Theorem 5.1.

Let $u$ be the solution of (1.1) and $U_{n}$, $n=0,\dots,N$, the solution of
the discrete system (2.9) (or (2.12)). The error $e=u-\hat{U}$ satisfies the
bound

$\begin{split}\int_{0}^{t}\|e(\tau)\|^{2}+|e(\tau)|_{1}^{2}\,d\tau\leq&C_{t,\beta}^{1}\int_{0}^{t}(g_{\beta,t}(\tau))^{-1}\left\|f(\tau)-D_{t}^{\beta}\hat{U}(\tau)-A\hat{U}(\tau)\right\|^{2}\,d\tau\\\
&+C^{2}_{t,\beta}\|u^{0}-U_{0}\|^{2},\end{split}$

for $t>0$, and $g_{\beta,t}$, $C^{1}_{t,\beta}$, $C^{2}_{t,\beta}$ as in
Theorem 4.2 and Corollary 4.3

The next estimator controls the error in $L^{\infty}(H)\,.$

######  Theorem 5.2.

Let $u$ be the solution of (1.1) and $U_{n}$, $n=0,\dots,N$, the solution of
the discrete system (2.9) (or (2.12)). The error $e=u-\hat{U}$ satisfies the
bound

$\|e(t)\|\leq\|u^{0}-U_{0}\|+\frac{1}{\sin\theta}C_{\beta,\varphi}\int_{0}^{t}(t-\tau)^{\beta-1}\left\|f(\tau)-D_{t}^{\beta}\hat{U}(\tau)-A\hat{U}(\tau)\right\|d\tau,$

for any $\theta\in(0,\pi/2)$, $t\in(0,T]$, with
$\varphi=\max(0,\pi-(\pi-\theta)/\beta)$, and

$C_{\beta,\varphi}=\frac{1}{\pi}(\cos\varphi)^{\beta-1}\Gamma(1-\beta).$

Next we investigate the optimality of the above a posteriori error estimators.

###  5.2. Asymptotic behaviour of the estimators

In the extended numerical experiments presented in the last section we
demonstrate that both estimators accurately capture the asymptotic behaviour
of the error in several cases. In the rest of this section we consider the
question of asymptotic behaviour of the a posteriori estimators under certain
assumptions on the solution. We are particularly interested whether the
estimator can converge at the optimal convergence order of the L1 scheme
$O(\kappa^{2-\beta})$, i.e., better than linear. As a proof of concept, we
just consider the estimator in Theorem 5.2 for the L1 scheme. Given the
presence of the kernels in the estimators and the complicated a priori
analysis required, the discussion of the other cases is left for a future
work. For the L1 scheme $\hat{U}$ is the solution of

(5.7)
$\Pi_{1}D_{t}^{\beta}\hat{U}(t)+A\hat{U}(t)=\Pi_{1}f(t)+\hat{G}(t),\qquad
t\in[0,T],$

with initial data $\hat{U}(0)=U_{0}$, where $\hat{G}$ is the piecewise linear
correction term introduced above. The approximations $\hat{U}$ satisfy the
original problem (1.1) with a modified right-hand side

$D_{t}^{\beta}\hat{U}(t)+A\hat{U}(t)=\Pi_{1}f(t)+\hat{G}(t)+D_{t}^{\beta}\hat{U}(t)-\Pi_{1}D_{t}^{\beta}\hat{U}(t),\qquad
t\in(0,T],$

and the error $e=u-\hat{U}$ solves

(5.8)
$D_{t}^{\beta}e(t)+Ae(t)=f(t)-\Pi_{1}f(t)-\hat{G}(t)+\Pi_{1}D_{t}^{\beta}\hat{U}(t)-D_{t}^{\beta}\hat{U}(t),\qquad
t\in(0,T],$

and $e(0)=u^{0}-U_{0}$.

First of all let us investigate the effect of $\hat{G}$.

######  Lemma 5.3.

Let $\hat{G}$ be the correction function defined by (5.4). Then for
$t>\kappa_{0}$

$\int_{0}^{t}(g_{\beta,t}(\tau))^{-1}\|\hat{G}(\tau)\|^{2}d\tau\leq
C\|AU_{0}-f(0)\|^{2}\kappa_{0}^{1+\beta},$

and

$\begin{split}\int_{0}^{t}(t-\tau)^{\beta-1}\|\hat{G}(\tau)\|d\tau&\leq\frac{1}{\beta}\|AU_{0}-f(0)\|\left(t^{\beta}-(t-\kappa_{0})^{\beta}\right)\\\
&\leq C\kappa_{0}\|AU_{0}-f(0)\|(t{-\kappa_{0}})^{\beta-1},\end{split}$

for some constant $C>0$ depending on $\beta$.

###### Proof.

The first inequality follows from

$\begin{split}\int_{0}^{t}(g_{\beta,t}(\tau))^{-1}\|\hat{G}(\tau)\|^{2}d\tau&\leq\Gamma(1-\beta)\|AU_{0}-f(0)\|^{2}\int_{0}^{\kappa_{0}}((t-\tau)^{-\beta}+\tau^{-\beta})^{-1}d\tau\\\
&\leq\frac{1}{\beta+1}\Gamma(1-\beta)\|AU_{0}-f(0)\|^{2}\kappa_{0}^{1+\beta},\end{split}$

where we used

$\int_{0}^{\kappa_{0}}((t-\tau)^{-\beta}+\tau^{-\beta})^{-1}d\tau\leq\int_{0}^{\kappa_{0}}\tau^{\beta}d\tau\leq\frac{1}{1+\beta}\kappa_{0}^{\beta+1}.$

The second inequality is obtained by direct computation and the estimate

$\frac{1}{\kappa_{0}}\left(t^{\beta}-(t-\kappa_{0})^{\beta}\right)\leq\beta\max_{\tau\in[t-\kappa_{0},t]}\tau^{\beta-1}=\beta(t-\kappa_{0})^{\beta-1}.$

∎

######  Remark 5.4.

This suggests that in general the $L^{2}$ norm of the error is no better than
$O(\kappa_{0}^{\frac{1+\beta}{2}})$ and the $L^{\infty}$ norm no better than
$O(t^{\beta-1}\kappa_{0})$. Hence, denoting by $\kappa=\max_{j}\kappa_{j}$, we
require at least that $\kappa_{0}=O(\kappa^{\frac{4-2\beta}{1+\beta}})$ in
order to obtain optimal assymptotic convergence order $O(\kappa^{2-\beta})$ in
the $L^{2}$ norm. In the $L^{\infty}$ norm, we require at least
$\kappa_{0}=O(\kappa^{2-\beta})$ to obtain optimal convergence for $t_{n}>c$
for any fixed constant $c>0$. Further,
$\kappa_{0}=O(\kappa^{\frac{2-\beta}{\beta}})$ is required if optimal
convergence is to be expected uniformly for $t_{n}>\kappa_{0}$; this is
compatible with the result from [39]; see (2.11).

If the data $f$ is smooth for $t>0$ then so is the solution $u$ with a
possible singularity at $t=0$; see (4.2).

######  Lemma 5.5.

Let $u\in{C[0,T]\cap}C^{2}(0,T]$ and denote by $\hat{u}$ the piecewise linear
interpolant of the data $u(t_{j})$, $j=0,\dots,n+1$ where
$0=t_{0}<t_{1}<\dots<t_{n+1}$. Then there exist $z_{0}\in(0,t_{1})$ and
$\tilde{z}_{j}\in(t_{j-1},t_{j+1})$, $j=1,\dots,n$, such that for $t\in
I_{n}$, $n\geq 1$,

$\begin{split}\left|\Pi_{1}D_{t}^{\beta}\hat{u}(t)-D_{t}^{\beta}\hat{u}(t)\right|\leq&C\Bigl{(}\kappa_{0}^{2}t^{-1-\beta}|u^{\prime}(z_{0})|\\\
&+\sum_{j=1}^{n-1}\kappa_{j}^{2}(\kappa_{j-1}+\kappa_{j})(t-t_{j})^{-\beta-1}|u^{\prime\prime}(\tilde{z}_{j})|+(\kappa_{n-1}+\kappa_{n})\kappa_{n}^{1-\beta}|u^{\prime\prime}(\tilde{z}_{n})|\Bigr{)}.\end{split}$

For $t\in I_{0}$

$\left|\Pi_{1}D_{t}^{\beta}\hat{u}(t)-D_{t}^{\beta}\hat{u}(t)\right|\leq
C\kappa_{0}^{1-\beta}|u^{\prime}(z_{0})|.$

###### Proof.

First note that there exist $z_{j}\in(t_{j},t_{j+1})$ and
$\tilde{z}_{j}\in(z_{j},z_{j+1})$ such that

$\frac{u(t_{1})-u(0)}{\kappa_{0}}=u^{\prime}(z_{0})$

and

$\begin{split}\frac{u(t_{j+1})-u(t_{j})}{\kappa_{j}}-\frac{u(t_{j})-u(t_{j-1})}{\kappa_{j}}&=u^{\prime}(z_{j})-u^{\prime}(z_{j-1})\\\
&=(z_{j}-z_{j-1})u^{\prime\prime}(\tilde{z}_{j}).\end{split}$

Using the representation (2.8) we hence have for $t\leq t_{n+1}$

$\begin{split}\Pi_{1}D_{t}^{\beta}\hat{u}(t)-D_{t}^{\beta}\hat{u}(t)=&(\Pi_{1}a_{0}(t)-a_{0}(t))u^{\prime}(z_{0})\\\
&+\sum_{j=1}^{n}(\Pi_{1}a_{j}(t)-a_{j}(t))(z_{j}-z_{j-1})u^{\prime\prime}(\tilde{z}_{j}).\end{split}$

For $t\in I_{n}$ and $j\leq n-1$ we have

$|\Pi_{1}a_{j}(t)-a_{j}(t)|\leq C\kappa_{n}^{2}(t-t_{j})^{-1-\beta}.$

Whereas for $t\in I_{n}$

$|\Pi_{1}a_{n}(t)-a_{n}(t)|\leq C\kappa_{n}^{1-\beta}.$

Combining the last three statements gives the first result, whereas the last
statement for $n=0$ gives the second required result. ∎

If the solution $u$ is smooth, namely $u\in C^{2}[0,T]$ it is shown in [21,
40] that for uniform time-steps $\kappa=\kappa_{j}$ optimal convergence order

$\|u(t_{n})-U_{n}\|\leq C\kappa^{2-\beta}$

is obtained. We investigate now if our estimator in Theorem 5.2 achieves this.

First of all we note that if $u$ is continuously differentiable, then
$D_{t}^{\beta}u(0)=0$ and consequently (1.1) implies $Au_{0}=f(0)$. Therefore
if we take as initial data $U_{0}=u^{0}$, the correction function vanishes
$\hat{G}\equiv 0$. Next, for $u$ to be smooth, in general $f$ is not smooth
but behaves as

$f(t)\sim f(0)+c_{1}t^{1-\beta}$

asymptotically as $t\rightarrow 0$. To understand the origin of this
singularity, simply substitute $u=u(0)+u^{\prime}(0)t$ in the fractional
equation (1.1).

Hence, in order to investigate the term due to $f$ in the a posteriori error
estimate we need the following lemma.

######  Lemma 5.6.

Let $\beta\in(0,1)$. Then for a uniform mesh and $n\geq 1$

$\int_{0}^{t_{n+1}}(t_{n+1}-\tau)^{\beta-1}|\Pi_{1}\tau^{1-\beta}-\tau^{1-\beta}|d\tau\leq
C\kappa^{2-\beta}(t_{n+1}-\kappa)^{\beta-1}.$

###### Proof.

By results on linear interpolation we have that

$|\Pi_{1}t^{1-\beta}-t^{1-\beta}|\leq\frac{1}{2}(1-\beta)\beta\kappa^{2}t_{n}^{-\beta-1}\qquad
t\in I_{n}$

and for $t\in I_{0}$

$\begin{split}|\Pi_{1}t^{1-\beta}-t^{1-\beta}|&\leq{\max_{t\in[0,\kappa]}|\kappa^{-\beta}t-t^{1-\beta}|}\\\
&{=}(1-\beta)^{\frac{1}{\beta}}((1-\beta)^{-1}-1)\kappa^{1-\beta}.\end{split}$

Hence for $t\in I_{n}=(t_{n},t_{n+1}]$

$\begin{split}\int_{0}^{t_{n+1}}(t_{n+1}-\tau)^{\beta-1}|\Pi_{1}\tau^{1-\beta}-\tau^{1-\beta}|d\tau\leq&C\kappa^{1-\beta}\int_{0}^{\kappa}(t_{n+1}-\tau)^{\beta-1}d\tau\\\
&+C\kappa^{2}\sum_{j=1}^{n}t_{j}^{-\beta-1}\int_{t_{j}}^{t_{j+1}}(t_{n+1}-\tau)^{\beta-1}d\tau\\\
\leq&C\kappa^{2-\beta}(t_{n+1}-\kappa)^{\beta-1}\\\
&+C\kappa^{3}\sum_{j=1}^{n-1}t_{j}^{-\beta-1}(t_{n+1}-t_{j+1})^{\beta-1}\\\
&+C\kappa^{2+\beta}t_{n}^{-\beta-1}\\\
=&C\kappa^{2-\beta}(t_{n+1}-\kappa)^{\beta-1}+C\kappa\sum_{j=1}^{n-1}j^{-\beta-1}(n-j)^{\beta-1}\\\
&+C\kappa^{2+\beta}t_{n}^{-\beta-1}\\\
\leq&C\kappa^{2-\beta}(t_{n+1}-\kappa)^{\beta-1}+C\kappa^{2-\beta}t_{n}^{\beta-1}\\\
&+C\kappa^{2+\beta}t_{n}^{-\beta-1}.\end{split}$

In the last step we used the following estimate from [25, Lemma 5.3]

$\sum_{j=1}^{n-1}j^{-\beta-1}(n-j)^{\beta-1}=O(n^{\beta-1}).$

As $\kappa^{2+\beta}t_{n}^{-\beta-1}\leq C\kappa^{2-\beta}t_{n}^{\beta-1}$ for
$t_{n}>\kappa$ we have the required result. ∎

Finally we investigate the term due to
$\Pi_{1}D_{t}^{\beta}\hat{U}(\tau)-D_{t}^{\beta}\hat{U}(\tau)$.

######  Lemma 5.7.

Let $u\in C^{2}[0,T]$, then for $t>\kappa$

$\int_{0}^{t}(t-\tau)^{\beta-1}\left|\Pi_{1}D_{t}^{\beta}\hat{u}(\tau)-D_{t}^{\beta}\hat{u}(\tau)\right|d\tau\leq
C\kappa^{2-\beta}(t-\kappa)^{\beta-1}.$

###### Proof.

For $t\in I_{n}$ and $n>0$, Lemma 5.5 implies

$\begin{split}\left|\Pi_{1}D_{t}^{\beta}\hat{u}(t)-D_{t}^{\beta}\hat{u}(t)\right|&\leq
C\left(\kappa^{2}t^{-1-\beta}+\kappa^{3}\sum_{j=1}^{n-1}(t-t_{j})^{-\beta-1}+\kappa^{2-\beta}\right)\\\
&\leq
C\left(\kappa^{2}t^{-1-\beta}+\kappa^{2}\int_{0}^{t_{n-1}}(t-\tau)^{-\beta-1}d\tau+\kappa^{2-\beta}\right)\\\
&\leq C\left(\kappa^{2}t^{-1-\beta}+\kappa^{2-\beta}\right).\end{split}$

For $t\in I_{0}$

$\left|\Pi_{1}D_{t}^{\beta}\hat{u}(t)-D_{t}^{\beta}\hat{u}(t)\right|\leq
C\kappa^{1-\beta}.$

Hence

$\begin{split}\int_{0}^{t}(t-\tau)^{\beta-1}\left|\Pi_{1}D_{t}^{\beta}\hat{u}(\tau)-D_{t}^{\beta}\hat{u}(\tau)\right|d\tau\leq&C\left(\kappa^{1-\beta}\int_{0}^{\kappa}(t-\tau)^{\beta-1}d\tau\right.\\\
&+\kappa^{2-\beta}\int_{\kappa}^{t}(t-\tau)^{\beta-1}d\tau\\\
&+\left.\kappa^{2}\int_{\kappa}^{t}\tau^{-1-\beta}(t-\tau)^{\beta-1}d\tau\right)\\\
\leq&C\left(\kappa^{1-\beta}(t^{\beta}-(t-\kappa)^{\beta})\right.\\\
&+\left.\kappa^{2}\int_{\kappa}^{t}\tau^{-1-\beta}(t-\tau)^{\beta-1}d\tau\right)\\\
\leq&C\kappa^{2-\beta}(t-\kappa)^{\beta-1}.\end{split}$

∎

##  6\. Numerical experiments

###  6.1. Fractional differential equation with $f\equiv 0$

First we consider the simple fractional ordinary differential equation

(6.1) $D_{t}^{\beta}u+\lambda u=0,\qquad u(0)=u^{0},\qquad t\in(0,T),$

with $\lambda>0$ a fixed constant.

The solution is given by $u(t)=E_{\beta}(-\lambda t^{\beta})u^{0}$ where

$E_{\beta}(z)=\sum_{m=0}^{\infty}\frac{z^{m}}{\Gamma(\beta m+1)}$

is the Mittag–Leffler function. Note that
$u\sim(1-\frac{1}{\Gamma(\beta+1)}t^{\beta})u^{0}$ for $t\rightarrow 0$, hence
already the first derivative of $u$ is unbounded at $t=0$.

We will compare the convergence of the L1 scheme and CQ with our estimators.
The meshes will be of the form

$t_{j}=T(j/N)^{k}$

with $k\geq 1$; the mesh is uniform if $k=1$ and graded towards 0 if $k>1$.
The exact initial data will be used, i.e., $U_{0}=u^{0}$.

The estimators in Theorem 5.1 and Theorem 5.2 require the numerical
computation of the outer integral, which we throughout this section compute
using a compound midpoint rule. The inner term $D_{t}^{\beta}\hat{U}$ is just
the L1 fractional derivative and is computed exactly.

We denote by

(6.2)
$E^{1}(t)\approx\left(\int_{0}^{t}|u(\tau)-\hat{U}(\tau)|^{2}d\tau\right)^{1/2}$

the exact error in $L^{2}$ norm approximated by the compound midpoint rule and
by

(6.3) $E^{2}(t)=|u(t)-\hat{U}(t)|$

the exact error at time $t$. The estimators, again approximated by the
compound midpoint rule, are denoted by $E^{1}_{\text{est}}(t)$ and
$E^{2}_{\text{est}}(t)$ corresponding to Theorem 5.1 and Theorem 5.2
respectively.

$N$ | $E^{1}$ | eoc | $E^{1}_{\text{est}}$ | eoc
---|---|---|---|---
10 | $0.01$ |  | $0.054$ |
20 | $0.0055$ | 0.89 | $0.029$ | 0.92
40 | $0.0029$ | 0.91 | $0.015$ | 0.93
80 | $0.0016$ | 0.92 | $0.0080$ | 0.93
160 | $0.00081$ | 0.94 | $0.0042$ | 0.92
320 | $0.00042$ | 0.95 | $0.0022$ | 0.92
(a) Uniform mesh with $\beta=0.8$.

$N$ | $E^{1}$ | eoc | $E^{1}_{\text{est}}$ | eoc
---|---|---|---|---
10 | $0.0069$ |  | $0.042$ |
20 | $0.0032$ | 1.13 | $0.019$ | 1.15
40 | $0.0014$ | 1.16 | $0.0083$ | 1.18
80 | $0.00064$ | 1.17 | $0.0036$ | 1.19
160 | $0.00028$ | 1.19 | $0.0016$ | 1.20
320 | $0.00012$ | 1.19 | $0.00069$ | 1.20
(b) Graded mesh with $k=2$ and $\beta=0.8$.

$N$ | $E^{1}$ | eoc | $E^{1}_{\text{est}}$ | eoc
---|---|---|---|---
$10$ | $0.01$ |  | $0.11$ |
$20$ | $0.0071$ | $0.56$ | $0.076$ | $0.57$
$40$ | $0.0048$ | $0.57$ | $0.051$ | $0.57$
$80$ | $0.0032$ | $0.59$ | $0.034$ | $0.58$
$160$ | $0.0021$ | $0.6$ | $0.023$ | $0.58$
(c) Uniform mesh with $\beta=0.2$.

$N$ | $E^{1}$ | eoc | $E^{1}_{\text{est}}$ | eoc
---|---|---|---|---
$10$ | $0.0011$ |  | $0.012$ |
$20$ | $0.00037$ | $1.6$ | $0.0038$ | $1.63$
$40$ | $0.00012$ | $1.64$ | $0.0012$ | $1.65$
$80$ | $3.8\times 10^{-5}$ | $1.67$ | $0.00038$ | $1.67$
$160$ | $1.2\times 10^{-5}$ | $1.7$ | $0.00012$ | $1.69$
$320$ | $3.5\times 10^{-6}$ | $1.71$ | $3.6\times 10^{-5}$ | $1.7$
(d) Graded mesh with $k=3$ and $\beta=0.2$.

Table 1. The numerical results for the L1 scheme including the error measured
in the $L^{2}$ norm at time $t=1/2$ and the estimated order of convergence for
differerent choices of $\beta$ and with uniform and graded meshes.

In Table 1 we show the results for the L1 scheme using the $L^{2}$ error
measure (6.2) evaluated at $t=1/2$ and for two values of the parameter
$\beta\in\\{0.2,0.8\\}$. We see that suboptimal convergence order is obtained
using the uniform mesh, but the optimal convergence order is recovered with
the graded mesh. In both cases the estimator converges at the correct order.
Remark 5.4 predicts the requirement for the graded mesh
$k=\frac{4-2\beta}{1+\beta}$, i.e., $k=4/3$ for $\beta=0.8$ and $k=3$ for
$\beta=0.2$. For $\beta=0.8$ we see that for $k=2>4/3$ optimal convergence is
obtained. For $k=4/3$ and $\beta=0.8$ numerical experiments suggest that
asymptotically the optimal convergence is achieved but only for quite small
time steps; see Table 2. For $\beta=0.2$, we see that the observed rate is
slowly approaching the optimal rate 1.8 when the borderline grading $k=3$ is
used. Numerical experiments not reported here for $k$ smaller than $k=3$ for
$\beta=0.2$ give less than optimal convergence order.

$N$ | 160 | 320 | 640 | 1280
---|---|---|---|---
eoc. | 1.114 | 1.126 | 1.136 | 1.144
eoc.(est) | 1.130 | 1.136 | 1.141 | 1.146
Table 2. Estimated convergence order in the $L^{2}$ norm for the L1 scheme.
Here $\beta=0.8$ and $k=4/3$, which is the grading suggested by Remark 5.4.

The corresponding results for CQ based on the backward Euler scheme are shown
in Table 3. Recall that the optimal convergence order for the L1 scheme and
the CQ is different; see Remark 3.2. It is better than linear
$O(\kappa^{2-\beta})$ for the L1 scheme and linear for CQ. Further, linear
convergence is obtained for the CQ scheme for $\beta>0.5$ even with the
uniform time-step, hence we only consider $\beta=0.2$ in the numerical
results. Estimate (3.6) predicts a convergence order
$O(\kappa^{\beta+1/2})=O(\kappa^{0.7})$. In Table 3, the exact estimated order
of convergence is slowly approaching $0.7$. However, our estimator converges
at a slighly lower rate or $0.6$. The latter fits with Remark 5.4 which
predicts $O(\kappa^{\frac{1+\beta}{2}})=O(\kappa^{0.6})$. Again, for $k=2$,
optimal convergence is obtained, i.e., in the case of CQ, linear convergence.

$N$ | $E^{1}$ | eoc | $E^{1}_{\text{est}}$ | eoc
---|---|---|---|---
$10$ | $0.0069$ |  | $0.11$ |
$20$ | $0.0047$ | $0.55$ | $0.072$ | $0.58$
$40$ | $0.0032$ | $0.57$ | $0.048$ | $0.58$
$80$ | $0.0021$ | $0.59$ | $0.032$ | $0.58$
$160$ | $0.0014$ | $0.60$ | $0.021$ | $0.59$
$320$ | $0.00091$ | $0.61$ | $0.014$ | $0.59$
$640$ | $0.00059$ | $0.62$ | $0.0094$ | $0.59$
$1280$ | $0.00038$ | $0.63$ | $0.0063$ | $0.59$
(a) Uniform mesh with $\beta=0.2$.

$N$ | $E^{1}$ | eoc | $E^{1}_{\text{est}}$ | eoc
---|---|---|---|---
$10$ | $0.0041$ |  | $0.031$ |
$20$ | $0.0023$ | $0.86$ | $0.014$ | $1.1$
$40$ | $0.0012$ | $0.9$ | $0.0067$ | $1.1$
$80$ | $0.00064$ | $0.93$ | $0.0032$ | $1.09$
$160$ | $0.00033$ | $0.95$ | $0.0015$ | $1.08$
(b) Graded mesh with $k=2$ and $\beta=0.2$.

Table 3. The numerical results for the backward Euler based CQ including the
error measured in the $L^{2}$ norm at time $t=1/2$ and the estimated order of
convergence with uniform and graded meshes. The fractional power is
$\beta=0.2$ throughout.

The corresponding results in the $L^{\infty}$ error measure (6.3) evaluated at
$t=1/2$ are shown in Table 4 and Table 5. Again the estimator converges at the
correct order. As indicated before, convolution quadrature does not require a
graded mesh to reach its optimal convergence order, i.e., linear, in this norm
at a fixed time away from the singularity.

$N$ | $E^{2}$ | eoc | $E^{2}_{\text{est}}$ | eoc
---|---|---|---|---
10 | $0.012$ |  | $0.039$ |
20 | $0.0060$ | 0.95 | $0.020$ | 0.96
40 | $0.0030$ | 1.0 | $0.010$ | 0.97
80 | $0.0015$ | 1.0 | $0.0053$ | 0.98
160 | $0.00076$ | 1.0 | $0.0027$ | 0.98
320 | $0.00038$ | 1.0 | $0.0014$ | 0.98
(a) Uniform mesh with $\beta=0.8$.

$N$ | $E^{2}$ | eoc | $E^{2}_{\text{est}}$ | eoc
---|---|---|---|---
10 | $0.011$ |  | $0.037$ |
20 | $0.0051$ | 1.14 | $0.017$ | 1.12
40 | $0.0023$ | 1.17 | $0.0076$ | 1.16
80 | $0.0010$ | 1.18 | $0.0034$ | 1.18
160 | $0.00044$ | 1.19 | $0.0015$ | 1.19
320 | $0.00019$ | 1.19 | $0.00065$ | 1.19
(b) Graded mesh with $k=2$ and $\beta=0.8$.

Table 4. The numerical results for the L1 scheme including the error measured
in the $L^{\infty}$ norm at time $t=1/2$ and the estimated order of
convergence for differerent choices of $\beta$ and with uniform and graded
meshes.

$N$ | $E^{2}$ | eoc | $E^{2}_{\text{est}}$ | eoc
---|---|---|---|---
$10$ | $0.0088$ |  | $0.027$ |
$20$ | $0.0045$ | $0.98$ | $0.014$ | $0.98$
$40$ | $0.0022$ | $0.99$ | $0.007$ | $0.99$
$80$ | $0.0011$ | $1$ | $0.0035$ | $0.99$
$160$ | $0.00056$ | $1$ | $0.0018$ | $0.99$
(a) Uniform mesh with $\beta=0.8$.

$N$ | $E^{2}$ | eoc | $E^{2}_{\text{est}}$ | eoc
---|---|---|---|---
$10$ | $0.0026$ |  | $0.016$ |
$20$ | $0.0013$ | $1.0$ | $0.0077$ | $1.0$
$40$ | $0.00064$ | $1.0$ | $0.0038$ | $1.0$
$80$ | $0.00032$ | $1.0$ | $0.0019$ | $1.0$
$160$ | $0.00016$ | $1.0$ | $0.00095$ | $1.0$
(b) Uniform mesh with $\beta=0.2$.

Table 5. The numerical results for the backward Euler based CQ including the
error measured in the $L^{\infty}$ norm at time $t=1/2$ and the estimated
order of convergence with uniform meshe for $\beta=0.8$ and $\beta=0.2$.

###  6.2. Fractional differential equation without an analytic solution

For completeness, we investigate a fractional differential equation where an
analytic expression for the solution does not exist:

(6.4) $D_{t}^{\beta}u(t)+u(t)=2\cos(t),\qquad u(0)=1,\qquad t\in(0,T).$

We repeat just a single experiment with the L1 scheme and $\beta=0.8$. Only
the error $E^{2}(t)$ for $t=1/2$ is investigated. The results for the uniform
and graded meshes are shown Table 6 with the results analogous to the
corresponding experiment in the previous section (see Table 4). As the exact
solution is not available, we make use of a graded, fine mesh with $k=2$ to
obtain that $u(1/2)=1.3877\dots$, correct to the digits shown.

$N$ | $E^{2}$ | eoc | $E^{2}_{\text{est}}$ | eoc
---|---|---|---|---
10 | $0.015$ |  | $0.047$ |
20 | $0.0072$ | 1.0 | $0.024$ | 0.99
40 | $0.0035$ | 1.0 | $0.012$ | 1.0
80 | $0.0017$ | 1.0 | $0.0060$ | 1.0
160 | $0.00084$ | 1.0 | $0.0030$ | 1.0
(a) Uniform mesh with $\beta=0.8$.

$N$ | $E^{2}$ | eoc | $E^{2}_{\text{est}}$ | eoc
---|---|---|---|---
10 | $0.016$ |  | $0.050$ |
20 | $0.0071$ | 1.14 | $0.023$ | 1.13
40 | $0.0032$ | 1.17 | $0.010$ | 1.16
80 | $0.0014$ | 1.19 | $0.0045$ | 1.18
160 | $0.00060$ | 1.21 | $0.0020$ | 1.19
(b) Graded mesh with $k=2$ and $\beta=0.8$.

Table 6. The numerical results for the L1 scheme applied to the problem (6.4)
for which we have no exact solution. We include the error measured in the
$E^{2}$ norm (6.3) at time $t=1/2$ and the estimated order of convergence for
$\beta=0.8$ with uniform and graded meshes.

###  6.3. Fractional differential equation with non-smooth $f$

Including a simple correction term can recover optimal convergence if $f$ is
smooth for both the L1 scheme [42] and convolution quadrature [16]. However,
if $f$ is not globally smooth, non-uniform time-steps are needed for both the
schemes.

To illustrate this we consider the fractional differential equation

(6.5) $D_{t}^{\beta}u+\lambda u=f,\qquad u(0)=u^{0},\qquad t\in(0,T),$

with $\lambda>0$ a fixed constant and $f$ chosen so that the exact solution is
given by

$u(t)=1+t^{\beta}+H(t-r)t^{\beta},$

where $H(\cdot)$ is the Heaviside function and $r\in(0,T)$ some constant. In
this case the right-hand side $f$ is piecewise smooth

$f(t)=\lambda u(t)+\Gamma(\beta+1)(1+H(t-r)).$

Instead of a priori defining a graded mesh towards $t=0$ and $t=r$, we will
adaptively construct a non-uniform mesh using the developed estimators. We
denote again by $E^{2}_{\text{est}}$ the $L^{\infty}$ estimate and given a
parameter $\theta\in(0,1)$ we mark the interval $I_{n}$ if

$E^{2}_{\text{est}}(t_{n+1})\geq\theta\max_{j=1,\dots,N}E^{2}_{\text{est}}(t_{j}).$

To obtain a refined mesh, each marked interval is split into two.

To investigate convergence, we use as the error measure the maximum error

$e^{\max}=\max_{j=1,\dots,N}E^{2}(t_{j})$

and the corresponding a posteriori error estimate

$e^{\max}_{\text{est}}=\max_{j=1,\dots,N}E^{2}_{\text{est}}(t_{j}).$

Figure 1. On the left we compare the convergence of the L1 adaptive scheme
with the uniform discretization. The fractional power is $\beta=0.6$. The
time-steps chosen adaptively are shown on the right. Refinement near the
singularities at $t=0$ and $t=0.28$ can clearly be seen. Figure 2. Comparison
of the convergence of the CQ adaptive scheme compared with the uniform
discretization. The fractional power is $\beta=0.6$.

We perform experiments with $\beta=0.6$, $\theta=0.75$, and the singularity at
$r=0.28$. The comparison of the convergence in the above norms of the uniform
schemes with the above described adaptive scheme are shown on the left of
Figure 1 for the L1 scheme and in Figure 2 for CQ. We see that the adaptive
scheme in both cases does considerably better than the scheme with a uniform
time step. The difference would be much more pronounced for a smaller $\beta$
and hence less smooth $u$. On the right of Figure 1 we plot the time-steps
chosen by the scheme. Clearly the adaptive scheme was able to locate the
singularities near $t=0$ and $t=r$.

###  6.4. Subdiffusion equation with FEM in space

In this section we consider a subdiffusion equation on a 1D domain
$\Omega=(0,\pi)$: Find $u(t)\in H^{1}_{0}(\Omega)$ such that

(6.6) $\displaystyle D_{t}^{\beta}u-\Delta u$ $\displaystyle=f$
$\displaystyle\text{for }(t,x)\in(0,T)\times\Omega,$ $\displaystyle u(t,x)$
$\displaystyle=0$ $\displaystyle\text{for }(t,x)\in(0,T)\times\partial\Omega,$
$\displaystyle u(0,x)$ $\displaystyle=u^{0}(x)$ $\displaystyle\text{for
}x\in\Omega.$

The right-hand side $f$ is chosen so that the exact solution is given by

$u(t,x)=(1+t^{\beta}+H(t-r)t^{\beta})\sin x$

with $r\in(0,T)$ some constant. In this case the right-hand side $f$ is
piecewise smooth in time

$f(t,x)=u(t,x)+\Gamma(\beta+1)(1+H(t-r))\sin x.$

Choosing a finite dimensional subspace $V_{h}\subset H^{1}_{0}(\Omega)$, the
Galerkin discretization of (6.6) results in the problem: Find
$\mathbf{u}(t)\in V_{h}$

(6.7) $\displaystyle D_{t}^{\beta}\mathbf{u}+A_{h}\mathbf{u}$
$\displaystyle=P_{h}f$ $\displaystyle\text{for }t\in(0,T)$
$\displaystyle\mathbf{u}(0)$ $\displaystyle=P_{h}u^{0}.$

$A_{h}:V_{h}\rightarrow V_{h}$ is the Galerkin discretization of the Dirichlet
Laplacian defined by

$(A_{h}u,v)_{L^{2}(\Omega)}=(\nabla u,\nabla v)_{L^{2}(\Omega)}\qquad\text{for
all }u,v\in V_{h}.$

The operator $P_{h}:L^{2}(\Omega)\rightarrow V_{h}$ is the $L^{2}$-projection
defined by

$(P_{h}u,v)_{L^{2}(\Omega)}=(u,v)_{L^{2}(\Omega)}\qquad\text{for all }v\in
V_{h}.$

The operator $A_{h}$ satisfies all the assumptions required so that this
problem fits into the setting (1.1) with $A=A_{h}$, $V=H=V_{h}$. The above
setting is used for convenience and in order to demonstrate the behaviour of
the time-estimators derived in this paper. A complete a posteriori and
adaptive treatment of (6.6) will require dynamical change of the finite
element spaces (i.e. variable with $n$ discrete spaces $V_{h}^{n}$) in the
scheme definition and an analysis taking into account the spatial
discretisation error. Given the “multistep nature” of the time-discretisation
of the fractional equation (6.6), even the definition of the scheme requires
special attention. Still, point-wise representations of the fully-discrete
methods are possible using the elliptic reconstruction, [29]. In this case the
application of the framework derived herein will be applicable. However, this
is not a straightforward task, compare with [28, 20] for parabolic problems,
and it will be the subject of a forthcoming work.

In the numerical experiments we fix $V_{h}$ to a piecewise quadratic finite
element space on a uniform spatial mesh with mesh-width $h>0$. The
discretization in time is achieved using the adaptive L1 scheme described in
the previous section up to time $T=1$. As the error measure we use

(6.8) $e^{\max}=\max_{j=1,\dots,N}\|P_{h}u(t_{j})-U_{j}\|_{L^{2}(\Omega)},$

whereas the a posteriori measure is

(6.9) $e^{\max}_{\text{est}}=\max_{j=1,\dots,N}E^{2}_{\text{est}}(t_{j}),$

where $E^{2}_{\text{est}}(t_{j})$ is obtained using Theorem 5.2.

The numerical results for $\beta=0.2$ and $r=0.28$ are shown in Figure 3. In
the mark-and-refine scheme we again set $\theta=0.75$. Similar results as in
the previous section are obtained.

Figure 3. Numerical results for the subdiffusion equation (6.6) with
$\beta=0.2$. On the left we plot the convergence of the adaptive scheme in the
maximum norm (6.8) and (6.9) and compare it with the scheme with uniform time-
step. On the right is a plot of the solution where in the visible grid
adaptive refinement towards the singularities at $t=0$ and $t=0.28$ can be
discerned.

##  Appendix A Non-uniform convolution quadrature

Let $K$ be a sectorial operator, i.e., for some $\theta\in(0,\pi/2)$ and
$\mu\in\mathbb{R}$, K is analytic in $\mathbb{C}\setminus\\{z\;;\;|\arg
z|<\pi-\theta\\}$ and satisfies

$|K(s)|\leq C|s|^{\mu}\qquad s\in\mathbb{C}\setminus\\{z\;;\;|\arg
z|<\pi-\theta\\}.$

Denoting by $k=\mathscr{L}^{-1}K$ the inverse Laplace transform of $K$, we
consider the convolution

$u(t)=\int_{0}^{t}k(t-\tau)g(\tau)d\tau,$

If $\mu<0$, $k$ is integrable otherwise the integral needs to be understood as
a Hadamard finite part integral or equivalently setting $K_{m}(s)=s^{-m}K(s)$
with $m=\lceil\mu\rceil$ we have

$u(t)=\frac{d^{m}}{dt^{m}}\int_{0}^{t}k_{m}(t-\tau)g(\tau)d\tau,$

where $k_{m}=\mathscr{L}^{-1}K_{m}$ .

Replacing $k$ with the inverse Laplace transform of $K$ and exchanging
integrals gives

(A.1) $u(t)=\frac{1}{2\pi i}\int_{\sigma\pm
i\infty}K(s)\int_{0}^{t}e^{s(t-\tau)}g(t)ds=\frac{1}{2\pi i}\int_{\sigma\pm
i\infty}K(s)y(t;s)ds,$

where $\sigma>0$ is a constant and $y$ solves the ODE

$y^{\prime}=sy+g\qquad y(0)=0.$

This calculation can be justified for $\mu<0$. For $\mu\geq 0$, it is a formal
argument that nevertheless leads to a well-defined numerical scheme.

Convolution quadrature is obtained by discretizing the ODE and substituting
the result in (A.1). Applying backward Euler discretization to the ODE gives

$\frac{1}{\kappa_{n-1}}(y_{n}-y_{n-1})=sy_{n}+g(t_{n})\qquad y_{0}=0.$

Solving the recursion gives

$\begin{split}y_{n+1}&=(1-s\kappa_{n})^{-1}\left[y_{n}+\kappa_{n-1}g(t_{n+1})\right]\\\
&=\sum_{j=0}^{n}\kappa_{j}\prod_{k=j}^{n}(1-s\kappa_{k})^{-1}g(t_{j+1}),\end{split}$

Hence

(A.2) $u_{n+1}=\sum_{j=0}^{n}\omega^{\text{\tiny CQ}}_{n,j}(K)g(t_{j+1})$

where

$\begin{split}\omega^{\text{\tiny CQ}}_{n,j}(K)&=\kappa_{j}\frac{1}{2\pi
i}\int_{\sigma-i\infty}^{\sigma+i\infty}K(s)\prod_{k=j}^{n}(1-s\kappa_{k})^{-1}ds\\\
&=\kappa_{j}\frac{1}{2\pi
i}\int_{\mathcal{C}}K(s)\prod_{k=j}^{n}(1-s\kappa_{k})^{-1}ds,\end{split}$

where $\mathcal{C}$ is a negatively oriented contour contained in the right-
hand complex plane surrounding the poles at $s=\kappa_{k}^{-1}$. In this form,
we see that no condition on the growth parameter $\mu$ of $K$ is needed for
the weights to be well-defined.

An important property of convolution quadrature is the composition rule. Let
$K(s)=K_{1}(s)K_{2}(s)$, then

$u_{n+1}=\sum_{j=0}^{n}\omega^{\text{\tiny
CQ}}_{n,j}(K)g(t_{j+1})=\sum_{j=0}^{n}\omega^{\text{\tiny
CQ}}_{n,j}(K_{2})\sum_{\ell=0}^{j}\omega^{\text{\tiny
CQ}}_{j,\ell}(K_{1})g(t_{\ell+1}).$

This property was shown in [23] for $K_{1}(s)=K_{2}^{-1}(s)$ by using
properties of divided differences. We explain the steps for general $K_{1}$
and $K_{2}$.

$\begin{split}\omega^{\text{\tiny
CQ}}_{n,j}(K)&=\kappa_{j}\prod_{k=j}^{n}(-\kappa_{k})^{-1}\frac{1}{2\pi
i}\int_{\mathcal{C}}K(s)\frac{1}{\prod_{k=j}^{n}(\kappa_{k}^{-1}-s)}ds\\\
&=-\left(\prod_{k=j+1}^{n}(-\kappa_{k})^{-1}\right)\frac{1}{2\pi
i}\int_{\mathcal{C}}K(s)\frac{1}{\prod_{k=j}^{n}(s-\kappa_{k}^{-1})}ds\\\
&=\left(\prod_{k=j+1}^{n}(-\kappa_{k})^{-1}\right)[\kappa_{j}^{-1},\dots,\kappa_{n}^{-1}]K,\end{split}$

where $[x_{0},\dots,x_{n}]f$ are Newton’s divided differences. For the
equivalence of the above contour integral and divided differences see [9,
Equation (51)]; note that the negative sign is absorbed due to the contour
$\mathcal{C}$ being negatively oriented.

To show the composition rule we need the Leibniz rule for divided differences
[9, Corollary 28]

$[x_{0},\dots,x_{n}](fg)=\sum_{k=0}^{n}[x_{0},\dots,x_{k}]f[x_{k},\dots,x_{n}]g.$

Hence

$\begin{split}\omega^{\text{\tiny
CQ}}_{n,j}(K)&=\left(\prod_{k=j+1}^{n}(-\kappa_{k})^{-1}\right)[\kappa_{j}^{-1},\dots,\kappa_{n}^{-1}]K_{1}K_{2}\\\
&=\sum_{\ell=j+1}^{n}\left(\prod_{k=j+1}^{n}(-\kappa_{k})^{-1}\right)[\kappa_{j}^{-1},\dots,\kappa_{\ell}^{-1}]K_{1}[\kappa_{\ell}^{-1},\dots,\kappa_{n}^{-1}]K_{2}\\\
&=\sum_{\ell=j+1}^{n}\left(\prod_{k=j+1}^{\ell}(-\kappa_{k})^{-1}\right)[\kappa_{j}^{-1},\dots,\kappa_{\ell}^{-1}]\left(\prod_{k=\ell+1}^{n}(-\kappa_{k})^{-1}\right)K_{1}[\kappa_{\ell}^{-1},\dots,\kappa_{n}^{-1}]K_{2}\\\
&=\sum_{\ell=j+1}^{n}\omega^{\text{\tiny
CQ}}_{\ell,j}(K_{1})\omega^{\text{\tiny CQ}}_{n,\ell}(K_{2}).\end{split}$

Substituting this into (A.2) gives the composition rule.

Let us now apply all this to the fractional derivative where
$K(s)=s^{\beta}=K_{1}(s)K_{2}(s)$ with $K_{1}(s)=s^{\beta-1}$, $K_{2}(s)=s$.
Using the divided difference definition we see that

$\omega^{\text{\tiny CQ}}_{n,n}(K_{2})=\kappa_{n}^{-1},\;\omega^{\text{\tiny
CQ}}_{n,n-1}(K_{2})=-\kappa_{n}^{-1},\omega^{\text{\tiny CQ}}_{n,j}=0,\qquad
j<n-1$

whereas $\omega^{\text{\tiny CQ}}_{n,j}(K_{1})$ are, by definition, as in
(3.2).

The definition of convolution used in this section gives rise to the Riemann-
Liouville derivative

${}^{\text{RL}}_{\;\;\;0}D_{t}^{\beta}u=\frac{d}{dt}\frac{1}{\Gamma(1-\beta)}\int_{0}^{t}(t-\tau)^{-\beta}u(\tau)d\tau.$

In order to apply convolution quadrature to the Caputo derivative we rewrite
it as a Riemann-Liouville derivative using the identity

${}^{\text{C}}_{0}D_{t}^{\beta}u=^{\text{RL}}_{\;\;\;0}D_{t}^{\beta}(u-u(0))$

valid for sufficiently smooth $u$. Hence given a sequence of values
$U_{0},\dots,U_{N}$ we apply convolution quadrature to the sequence
$U_{j}-U_{0}$, $j=0,\dots,N$ giving

$\begin{split}\sum_{j=0}^{n}\omega^{\text{\tiny
CQ}}_{n,j}(K)(U_{j+1}-U_{0})&=\sum_{j=0}^{n}\omega^{\text{\tiny
CQ}}_{n,j}(K_{1})\sum_{\ell=0}^{j}\omega^{\text{\tiny
CQ}}_{j,\ell}(K_{2})(U_{\ell+1}-U_{0})\\\ &=\sum_{j=0}^{n}\omega^{\text{\tiny
CQ}}_{n,j}(K_{1})\frac{1}{\kappa_{j}}\left(U_{j+1}-U_{j}\right),\end{split}$

using the composition rule. This is the representation (3.1) used in this
paper.

## References

  * [1] G. Akrivis, C. Makridakis, and R. H. Nochetto.  Optimal order a posteriori error estimates for a class of Runge-Kutta and Galerkin methods.  Numer. Math., 114(1):133–160, 2009.
  * [2] K. Baker and L. Banjai.  Numerical analysis of a wave equation for lossy media obeying a frequency power law.  IMA Journal of Numerical Analysis, 05 2021.  drab028.
  * [3] L. Banjai and M. López-Fernández.  Efficient high order algorithms for fractional integrals and fractional differential equations.  Numer. Math., 141(2):289–317, 2019.
  * [4] L. Banjai and F.-J. Sayas.  Integral equation methods for evolutionary PDE: A convolution quadrature approach.  Springer Series in Computational Mathematics. Springer, To be published in 2022.
  * [5] E. Bänsch, F. Karakatsani, and C. G. Makridakis.  A posteriori error estimates for fully discrete schemes for the time dependent Stokes problem.  Calcolo, 55(2):Art. 19, 32, 2018.
  * [6] E. Bazhlekova.  The abstract Cauchy problem for the fractional evolution equation.  Fract. Calc. Appl. Anal., 1(3):255–270, 1998.
  * [7] H. Brunner.  Collocation methods for Volterra integral and related functional differential equations, volume 15 of Cambridge Monographs on Applied and Computational Mathematics.  Cambridge University Press, Cambridge, 2004.
  * [8] A. Cangiani, E. H. Georgoulis, I. Kyza, and S. Metcalfe.  Adaptivity and blow-up detection for nonlinear evolution problems.  SIAM J. Sci. Comput., 38(6):A3833–A3856, 2016.
  * [9] C. de Boor.  Divided differences.  Surv. Approx. Theory, 1:46–69, 2005.
  * [10] W. Feller.  An introduction to probability theory and its applications. Vol. II.  John Wiley & Sons, Inc., New York-London-Sydney, 1966.
  * [11] E. H. Georgoulis, O. Lakkis, C. G. Makridakis, and J. M. Virtanen.  A posteriori error estimates for leap-frog and cosine methods for second order evolution problems.  SIAM J. Numer. Anal., 54(1):120–136, 2016.
  * [12] R. Gorenflo, F. Mainardi, D. Moretti, and P. Paradisi.  Time fractional diffusion: a discrete random walk approach.  Nonlinear Dynam., 29(1-4):129–143, 2002.  Fractional order calculus and its applications.
  * [13] W. Hackbusch.  Integral equations, volume 120 of International Series of Numerical Mathematics.  Birkhäuser Verlag, Basel, 1995.  Theory and numerical treatment, Translated and revised by the author from the 1989 German original.
  * [14] Y. Hatano and N. Hatano.  Dispersive transport of ions in column experiments: An explanation of long-tailed profiles.  Water Resources Research, 34(5):1027–1033, 1998.
  * [15] B. Jin, R. Lazarov, and Z. Zhou.  An analysis of the L1 scheme for the subdiffusion equation with nonsmooth data.  IMA J. Numer. Anal., 36(1):197–221, 2016.
  * [16] B. Jin, B. Li, and Z. Zhou.  Correction of high-order BDF convolution quadrature for fractional evolution equations.  SIAM J. Sci. Comput., 39(6):A3129–A3152, Jan. 2017.
  * [17] T. Kato.  Perturbation theory for linear operators.  Classics in Mathematics. Springer-Verlag, Berlin, 1995.  Reprint of the 1980 edition.
  * [18] T. Katsaounis and I. Kyza.  A posteriori error analysis for evolution nonlinear Schrödinger equations up to the critical exponent.  SIAM J. Numer. Anal., 56(3):1405–1434, 2018.
  * [19] A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo.  Theory and applications of fractional differential equations, volume 204 of North-Holland Mathematics Studies.  Elsevier Science B.V., Amsterdam, 2006.
  * [20] O. Lakkis and C. Makridakis.  Elliptic reconstruction and a posteriori error estimates for fully discrete linear parabolic problems.  Math. Comp., 75(256):1627–1658, 2006.
  * [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] M. López-Fernández, C. Lubich, and A. Schädle.  Adaptive, fast, and oblivious convolution in evolution equations with memory.  SIAM J. Sci. Comput., 30(2):1015–1037, 2008.
  * [23] M. López-Fernández and S. Sauter.  Generalized convolution quadrature based on Runge-Kutta methods.  Numer. Math., 133(4):743–779, 2016.
  * [24] C. Lubich.  Discretized fractional calculus.  SIAM J. Math. Anal., 17(3):704–719, 1986.
  * [25] C. Lubich.  Convolution quadrature and discretized operational calculus. I.  Numer. Math., 52(2):129–145, 1988.
  * [26] C. Lubich.  Convolution quadrature revisited.  BIT, 44(3):503–514, 2004.
  * [27] C. Lubich and C. Makridakis.  Interior a posteriori error estimates for time discrete approximations of parabolic problems.  Numer. Math., 124(3):541–557, 2013.
  * [28] C. Makridakis.  Space and time reconstructions in a posteriori analysis of evolution problems.  In ESAIM Proceedings. Vol. 21 (2007) [Journées d’Analyse Fonctionnelle et Numérique en l’honneur de Michel Crouzeix], volume 21 of ESAIM Proc., pages 31–44. EDP Sci., Les Ulis, 2007.
  * [29] C. Makridakis and R. H. Nochetto.  Elliptic reconstruction and a posteriori error estimates for parabolic problems.  SIAM J. Numer. Anal., 41(4):1585–1594, 2003.
  * [30] W. McLean.  Fast summation by interval clustering for an evolution equation with memory.  SIAM J. Sci. Comput., 34(6):A3039–A3056, 2012.
  * [31] R. Metzler and J. Klafter.  Boundary value problems for fractional diffusion equations.  Phys. A, 278(1-2):107–125, 2000.
  * [32] R. R. Nigmatullin.  The realization of the generalized transfer equation in a medium with fractal geometry.  physica status solidi (b), 133(1):425–430, Jan. 1986.
  * [33] R. H. Nochetto, G. Savaré, and C. Verdi.  A posteriori error estimates for variable time-step discretizations of nonlinear evolution equations.  Comm. Pure Appl. Math., 53(5):525–589, 2000.
  * [34] L. Oparnica and E. Süli.  Well-posedness of the fractional Zener wave equation for heterogeneous viscoelastic materials.  Fract. Calc. Appl. Anal., 23(1):126–166, 2020.
  * [35] NIST Digital Library of Mathematical Functions.  http://dlmf.nist.gov/, Release 1.1.1 of 2021-03-15.  F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
  * [36] H. E. Roman and P. A. Alemany.  Continuous-time random walks and the fractional diffusion equation.  J. Phys. A, 27(10):3407–3410, 1994.
  * [37] 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.
  * [38] D. Schötzau and T. P. Wihler.  A posteriori error estimation for $hp$-version time-stepping methods for parabolic partial differential equations.  Numer. Math., 115(3):475–509, 2010.
  * [39] 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.
  * [40] Z.-z. Sun and X. Wu.  A fully discrete difference scheme for a diffusion-wave system.  Appl. Numer. Math., 56(2):193–209, 2006.
  * [41] K. Šiškova.  Inverse Source Problems in Evolutionary PDE’s.  PhD thesis, Ghent University, 2018.
  * [42] 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, Jan. 2018.