# An $\alpha$-robust second-order accurate scheme for a subdiffusion
equation††thanks: School of Mathematics and Statistics, University of New
South Wales, Sydney, Australia
This work was supported by the Australian Research Council grant DP220101811.

Kassem Mustapha, William McLean and Josef Dick

###### Abstract

We investigate a second-order accurate time-stepping scheme for solving a
time-fractional diffusion equation with a Caputo derivative of order
$\alpha\in(0,1)$. The basic idea of our scheme is based on local integration
followed by linear interpolation. It reduces to the standard Crank–Nicolson
scheme in the classical diffusion case, that is, as $\alpha\to 1$. Using a
novel approach, we show that the proposed scheme is $\alpha$-robust and
second-order accurate in the $L^{2}(L^{2})$-norm, assuming a suitable time-
graded mesh. For completeness, we use the Galerkin finite element method for
the spatial discretization and discuss the error analysis under reasonable
regularity assumptions on the given data. Some numerical results are presented
at the end.

##  1 Introduction

We shall approximate the solution of the time-fractional diffusion equation

$\partial_{t}^{\alpha}u({\boldsymbol{x}},t)+\mathcal{A}u({\boldsymbol{x}},t)=f({\boldsymbol{x}},t)\quad\text{for
$({\boldsymbol{x}},t)\in\Omega\times(0,T]$,}$ (1.1)

subject to homogeneous Dirichlet boundary conditions, that is,
$u({\boldsymbol{x}},t)=0$ on $\partial\Omega\times(0,T]$, with
$u({\boldsymbol{x}},0)=u_{0}({\boldsymbol{x}})$ at the initial time level
$t=0.$ The spatial domain $\Omega\subset\mathbb{R}^{d}$ (with $d=1$, $2$, $3$)
is a convex polyhedron, $0<\alpha<1$, the time fractional Caputo derivative

$\partial_{t}^{\alpha}v(t):=\mathcal{I}^{1-\alpha}v^{\prime}(t)=\int_{0}^{t}\omega_{1-\alpha}(t-s)v^{\prime}(s)\,ds,\quad{\rm
with}\quad\omega_{1-\alpha}(t):=\frac{t^{-\alpha}}{\Gamma(1-\alpha)},$

where $v^{\prime}=\partial v/\partial t$ and $\Gamma$ denotes the gamma
function. We use the notation $\mathcal{I}v(t)$ for the standard time integral
of $v$ from $0$ to $t.$ In (1.1), ${\mathcal{A}}$ is an elliptic operator in
the spatial variables, defined by
${\mathcal{A}}w({\boldsymbol{x}})=-\nabla\cdot(\kappa({\boldsymbol{x}})\nabla
w)({\boldsymbol{x}})$. The diffusivity $\kappa\in L^{\infty}(\Omega)$
satisfies $0<\kappa_{\min}\leq\kappa$ on $\Omega,$ for some constant
$\kappa_{\min}$. For the error analysis, we also require that $\kappa\in
W^{1,\infty}(\Omega)$.

The presence of the nonlocal time fractional (Caputo) derivative in (1.1) and
the fact that the solution $u$ suffers from a weak singularity near $t=0$ have
a direct impact on the accuracy, and consequently the convergence rates, of
numerical methods. To overcome this difficulty, different approaches have been
applied including corrections, graded meshes, and convolution quadrature [6,
8, 14, 25, 30, 31]. Indeed, the numerical solutions for model problems of the
form (1.1), including _a priori_ and _a posteriori_ error analyses and fast
algorithms, were studied by various authors over the past fifteen years using
multiple approaches [1, 3, 5, 7, 9, 10, 12, 13, 18, 28, 32, 33, 34, 36, 37].
For more references and details, see the recent monograph by Jin and Zhou
[11].

In this work, we investigate rigorously the error from approximating the
solution of the initial-boundary value problem (1.1) using a uniform second-
order accurate time-stepping method. The latter is defined via a local time-
integration of problem (1.1) on each subinterval of the time mesh combined
with continuous piecewise linear interpolation. The proposed scheme is
identical to the piecewise-linear case of a discontinuous Petrov–Galerkin
method proposed in [22]. Therein, with $\tau$ being the maximum time mesh step
size, a suboptimal convergence rate of order $O(\tau^{(3-\alpha)/2})$ was
proved. A time-graded mesh (2.1) was employed to compensate for the singular
behaviour of the continuous solution at $t=0$. In the limiting case as
$\alpha\to 1$, the problem (1.1) reduces to the classical diffusion equation,
and the considered numerical scheme reduces to the classical Crank–Nicolson
method. In this case, $O(\tau^{(3-\alpha)/2})$ reduces to $O(\tau)$ which is
far from the optimal $O(\tau^{2})$ rate achieved in practice.

By using an innovative approach that relies on interesting implicit polynomial
interpolations and duality arguments, we show $O(\tau^{2})$ convergence,
whilst at the same time relaxing the imposed regularity assumptions from the
earlier analysis [22]. This convergence rate is $\alpha$-robust in the sense
that the constant in the error bound remains bounded as $\alpha\to 1$.
Implementation wise, although the proposed scheme is uniformly second-order
accurate, the computational cost is comparable to the well-known backward
Euler or $L1$ [16, 25, 29] methods, which are not even first-order accurate.

For completeness, we discretize the problem (1.1) over the spatial domain
$\Omega$ using the standard Galerkin finite element method (FEM), thereby
defining a fully discrete approximation to $u$. An additional error of order
$O(h^{2})$ is anticipated under certain regularity assumptions on the
continuous solution, where $h$ is the mximum spatial finite element mesh size.
This is proved via a concise approach that relies on the discrete version of
the earlier error analysis. To make this feasible, the solution of the
semidiscrete Galerkin finite element solution of problem (1.1) plays the role
of the comparison function.

Outline of the paper. In the next section, we define our time-stepping scheme,
introduce some notations and technical lemmas, and summarize the convergence
results in Theorem 2.1. The required regularity properties are also
highlighted. Section 3 proves some error bounds for the implicit piecewise-
linear interpolant $\widehat{u}$ defined in (2.6). Section 4 is devoted to
showing the second order of accuracy of the proposed time-stepping scheme via
a duality argument. In Section 5, we discretize in space via the Galerkin
finite element method and discuss the convergence of the fully discrete
solution. To support our theoretical findings, we present some numerical
results in Section 6. Finally, a short technical appendix derives an
$\alpha$-robust interpolation estimate.

##  2 Time-stepping scheme

This section is devoted to discretizing the model problem (1.1) over the time
interval $[0,T]$ through a second-order accurate method, and stating our main
convergence results. We begin by introducing some notations.

For $\ell\geq 0,$ the norm on $H^{\ell}(\Omega)$ is denoted by
$\|\cdot\|_{\ell}$. The Sobolev spaces $H^{\ell}(\Omega)$ and
$H^{1}_{0}(\Omega)$ are defined as usual, and 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 $\mathcal{A}$ on
$\Omega$. The inner product in $L^{2}(\Omega)$ is denoted by
$\langle\cdot,\cdot\rangle$, and the associated norm by $\|\cdot\|$. The
generic constant $C$ remains bounded for $0<\alpha\leq 1,$ and is independent
of the time mesh and the finite element mesh, but may depend on $\Omega$, $T$,
and other quantities, including $\kappa$, $u_{0}$ and $f$.

Define the time mesh $0=t_{0}<t_{1}<t_{2}<\ldots<t_{N}=T$ by

$t_{n}=(n\,\tau)^{\gamma},\quad\text{with $\tau=T^{1/\gamma}/N$ and
$\gamma\geq 1$,}\quad\text{for $0\leq n\leq N$,}$ (2.1)

and let $\tau_{n}=t_{n}-t_{n-1}$. Such a time-graded mesh has the properties
[19]

$t_{n}\leq 2^{\gamma}t_{n-1}\quad\text{and}\quad\gamma\tau
t_{n-1}^{1-1/\gamma}\leq\tau_{n}\leq\gamma\tau
t_{n}^{1-1/\gamma},\quad\text{for $n\geq 2$.}$ (2.2)

Integrating problem (1.1) over $I_{n}$ and then dividing by $\tau_{n}$ yields

$\frac{1}{\tau_{n}}\int_{I_{n}}\partial_{t}^{\alpha}u\,dt+\mathcal{A}\bar{u}_{n}=\bar{f}_{n},\quad{\rm
for}~{}~{}1\leq n\leq N,$ (2.3)

where $\bar{f}_{n}=\tau_{n}^{-1}\int_{I_{n}}f(t)\,dt$ denotes the average
value of a function $f$ over the time interval $I_{n}$. Motivated by (2.3),
for $t\in I_{n}$ and for $1\leq n\leq N,$ our semidiscrete approximate
solution $U(t)\approx u(t$ is defined by requiring that

$U(t)=\frac{t_{n}-t}{\tau_{n}}U^{n-1}+\frac{t-t_{n-1}}{\tau_{n}}U^{n},\quad
U^{n}:=U(t_{n}),$

with

$\frac{1}{\tau_{n}}\int_{I_{n}}\partial_{t}^{\alpha}U\,dt+\mathcal{A}U^{n-1/2}=\bar{f}_{n},\quad{\rm
with}~{}~{}U^{0}=U(0)=u_{0},$ (2.4)

where $U^{n-1/2}=\bar{U}_{n}=\tfrac{1}{2}(U^{n}+U^{n-1}).$ At the $n$th time
step, with $\tau_{n,\alpha}=\Gamma(3-\alpha)\tau_{n}^{\alpha}$, we have to
solve an elliptic problem for $W^{n}=U^{n}-U^{n-1}$, namely,

$(I+\tfrac{1}{2}\tau_{n,\alpha}\mathcal{A})W^{n}=\tau_{n,\alpha}\bar{f}_{n}-\tau_{n,\alpha}\mathcal{A}U^{n-1}-\tau_{n,\alpha}\sum_{j=1}^{n-1}(\tau_{n}\tau_{j})^{-1}b_{n,j}W^{j},$

where $b_{n,j}=a_{n,j}-a_{n,j+1},$ with
$a_{n,j}=\omega_{3-\alpha}(t_{n}-t_{j-1})-\omega_{3-\alpha}(t_{n-1}-t_{j-1}).$

If $\alpha\to 1$, then $\partial_{t}^{\alpha}u\to u^{\prime}$ and
$\partial_{t}^{\alpha}U\to U^{\prime}$, implying that our scheme reduces to
the Crank–Nicolson scheme for the classical diffusion equation.

Our convergence analysis relies on decomposing the error as

$\eta=u-U=\psi-\theta\quad\text{with}\quad\psi=u-\widehat{u}\quad\text{and}\quad\theta=U-\widehat{u},$
(2.5)

where $\widehat{u}$ is a continuous piecewise-linear function in time
satisfying

$\int_{I_{n}}\widehat{u}(t)\,dt=\int_{I_{n}}u(t)\,dt\quad\text{for $1\leq
n\leq N$,}\quad\text{with $\widehat{u}(0)=\widehat{u}^{0}=u_{0}$.}$ (2.6)

Alternatively, $\widehat{u}$ can be defined via
$\mathcal{I}\widehat{u}(t_{n})=\mathcal{I}u(t_{n})$ for $1\leq n\leq N$, with
$\widehat{u}(0)=u_{0}$, and we say that $\widehat{u}$ interpolates $u$
_implicitly_. The decomposition (2.5) of the error $\eta$ follows a well-known
pattern, but the novel choice of the piecewise linear function $\widehat{u}$
makes possible our improved error analysis under reasonable regularity
assumption. The continuous average of $u$ equals both the continuous and the
discrete average of $\widehat{u}$ on each time subinterval $I_{n}$. For
comparison, let $u_{I}$ denote the usual continuous piecewise-linear
interpolant to $u$, that is,

$u_{I}(t)=\frac{t_{n}-t}{\tau_{n}}\,u(t_{n-1})+\frac{t-t_{n-1}}{\tau_{n}}u(t_{n})\quad\text{for
$t\in I_{n}$,}$ (2.7)

and observe that $u_{I}$ and $u$ have the same _discrete_ average
$\tfrac{1}{2}(u(t_{n})+u(t_{n-1}))$ on each $I_{n}$, but their continuous
averages will differ unless $u$ is linear on $I_{n}$.

Subtracting (2.4) from (2.3) and using (2.6), we obtain

$\int_{I_{n}}\partial_{t}^{\alpha}(\psi-\theta)\,dt-\int_{I_{n}}\mathcal{A}\theta\,dt=0.$

Taking the $L^{2}(\Omega)$-inner product with a test function $\varphi\in
H^{1}_{0}(\Omega)$, and applying the divergence theorem, it follows that

$\int_{I_{n}}\langle\mathcal{I}^{1-\alpha}\theta^{\prime},\varphi\rangle\,dt+\int_{I_{n}}\langle\kappa\nabla\theta,\nabla\varphi\rangle\,dt=\int_{I_{n}}\langle\mathcal{I}^{1-\alpha}\psi^{\prime},\varphi\rangle\,dt.$
(2.8)

Choosing $\varphi=\theta^{\prime}$ and summing over $n$ yields

$\mathcal{I}(\langle\mathcal{I}^{1-\alpha}\theta^{\prime},\theta^{\prime}\rangle)(t_{n})+\mathcal{I}(\langle\kappa\nabla\theta,\nabla\theta^{\prime}\rangle)(t_{n})=\mathcal{I}(\langle\mathcal{I}^{1-\alpha}\psi^{\prime},\theta^{\prime}\rangle)(t_{n}).$

Since
$\mathcal{I}(\langle\kappa\nabla\theta,\nabla\theta^{\prime}\rangle)(t_{n})=\tfrac{1}{2}(\|\sqrt{\kappa}\nabla\theta(t_{n})\|^{2}-\|\sqrt{\kappa}\nabla\theta(0)\|^{2})=\tfrac{1}{2}\|\sqrt{\kappa}\nabla\theta(t_{n})\|^{2}$,

$\mathcal{I}(\langle\mathcal{I}^{1-\alpha}\theta^{\prime},\theta^{\prime}\rangle)(t_{n})\leq\mathcal{I}(\langle\mathcal{I}^{1-\alpha}\psi^{\prime},\theta^{\prime}\rangle)(t_{n}).$
(2.9)

To proceed in our analysis, we make use of the following technical lemma. For
the proof, we refer to Mustapha and Schötzau [23, Lemma 3.1 (iii)].

###### Lemma 2.1

For $0<\alpha\leq 1$ and $\epsilon>0$,

$\displaystyle\mathcal{I}(\langle\mathcal{I}^{1-\alpha}v,w\rangle)(t)$
$\displaystyle\leq\frac{1}{\alpha}\Big{(}\mathcal{I}(\langle\mathcal{I}^{1-\alpha}v,v\rangle)(t)\Big{)}^{1/2}\Big{(}\mathcal{I}(\langle\mathcal{I}^{1-\alpha}w,w\rangle)(t)\Big{)}^{1/2}$
$\displaystyle\leq\epsilon\mathcal{I}(\langle\mathcal{I}^{1-\alpha}v,v\rangle)(t)+\frac{1}{4\epsilon\alpha^{2}}\mathcal{I}(\langle\mathcal{I}^{1-\alpha}w,w\rangle)(t).$

For later use, by expanding $\langle\mathcal{I}^{1-\alpha}(v+w),v+w)\rangle$
then applying Lemma 2.1 with $\epsilon=1/(2\alpha)$ we deduce the inequality
in the next lemma.

###### Lemma 2.2

For $0<\alpha\leq 1$,

$\mathcal{I}(\langle\mathcal{I}^{1-\alpha}(v+w),v+w\rangle)(t)\leq(1+\alpha^{-1})\Big{(}\mathcal{I}(\langle\mathcal{I}^{1-\alpha}v,v\rangle)(t)+\mathcal{I}(\langle\mathcal{I}^{1-\alpha}w,w\rangle)(t)\Big{)}.$

We now apply Lemma 2.1 to the right-hand side of (2.9) with
$\epsilon=1/(2\alpha^{2})$. Multiplying through by $2$, and then cancelling
the similar terms, leads to the estimate below that will be used later in our
convergence analysis.

$\mathcal{I}(\langle\mathcal{I}^{1-\alpha}\theta^{\prime},\theta^{\prime}\rangle)(t_{n})\leq\frac{1}{\alpha^{2}}\mathcal{I}(\langle\mathcal{I}^{1-\alpha}\psi^{\prime},\psi^{\prime}\rangle)(t_{n}),$
(2.10)

Under reasonable regularity assumptions, a novel error analysis involving
implicit interpolations and a duality argument leads to the convergence
results in the next theorem. With $J=(0,T)$, an optimal $O(\tau^{2})$-rate of
convergence is achieved in the $L^{2}(J;L^{2}(\Omega))$-norm. Our numerical
results illustrate this in the stronger $L^{\infty}(J;L^{2}(\Omega))$-norm.
Moreover, our numerical results suggest that the condition on the graded mesh
exponent can be further relaxed. More precisely, instead of
$\gamma>\max\\{2/\sigma,(3-\alpha)/(2\sigma-\alpha)\\}$ it suffices to impose
$\gamma>2/\sigma$.

The developed error analysis requires the following regularity property [10,
Theorems 2.1 and 2.2], [27, Theorems 1 and 2], and [20]: for some $\sigma>0$,

$t\|u^{\prime}(t)\|+t^{2}\|u^{\prime\prime}(t)\|+t^{3}\|u^{\prime\prime\prime}(t)\|\leq
Ct^{\sigma}\quad\text{for $t>0$.}$ (2.11)

For example, if $f\equiv 0$ and $u_{0}\in\dot{H}^{r}(\Omega)$ with $1<r\leq
2$, then (2.11) holds true for $\sigma=r\alpha/2>\alpha/2$.

For a given time interval $Q$, the norms in
$L^{\infty}\bigl{(}Q;L^{2}(\Omega)\bigr{)}$ and
$L^{2}\bigl{(}Q;L^{2}(\Omega)\bigr{)}$ are respectively denoted by
$\|w\|_{Q}=\sup_{t\in Q}\|w(t)\|$ and
$\|w\|_{L^{2}(Q)}=\bigl{(}\int_{Q}\|w(t)\|^{2}\,dt\bigr{)}^{1/2}$.

###### Theorem 2.1

Let $u$ and $U$ be the solutions of (1.1) and (2.4), respectively. If the
graded time mesh exponent
$\gamma>\max\\{2/\sigma,(3-\alpha)/(2\sigma-\alpha)\\}$ and if the regularity
assumption (2.11) holds true with $\sigma>\alpha/2$, then for $0<\alpha<1,$ we
have $\|u-U\|_{L^{2}(J)}\leq C\alpha^{-2}\tau^{2}$.

Proof. The desired estimate follows from Lemma 3.4 and Theorem 4.1 below.
$\quad\Box$

##  3 Errors from implicit interpolations

In preparation for our convergence analysis, we now study the error from
approximating $u$ by $\widehat{u}$, and proceed to estimate $\|\psi\|$ and
$\mathcal{I}(\langle\mathcal{I}^{1-\alpha}\psi^{\prime},\psi^{\prime}\rangle)$.
These estimates assume that the regularity property (2.11) holds. For ease of
reference, we here introduce the parameter $\delta=\sigma-\frac{2}{\gamma}$
which will subsequently appear repeatedly. We start this section with the
following representation of the implicit interpolation error in the
approximation $u\approx\widehat{u}$ at $t=t_{n}$.

###### Lemma 3.1

For $1\leq n\leq N$, $\psi^{n}=\sum_{j=1}^{n}(-1)^{n+j+1}\Delta_{j}$ where

$\Delta_{j}=\frac{2}{\tau_{j}}\int_{I_{j}}(u-u_{I})\,dt=-\frac{1}{\tau_{j}}\int_{I_{j}}(t_{j}-t)(t-t_{j-1})u^{\prime\prime}(t)\,dt.$

Proof. Since $\int_{I_{n}}u_{I}\,dt=\tfrac{1}{2}\tau_{n}(u^{n}+u^{n-1})$, we
find using (2.6) that

$\psi^{n}=-\psi^{n-1}-\frac{2}{\tau_{n}}\int_{I_{n}}(u-u_{I})\,dt.$

The formula for $\psi^{n}$ then follows by induction on $n$, after noting that
$\psi^{0}=0$. Recalling the Peano kernel for the trapezoidal rule, we see that

$\int_{I_{j}}(u-u_{I})\,dt=-\frac{1}{2}\int_{I_{j}}(t_{j}-s)(s-t_{j-1})u^{\prime\prime}(s)\,ds,$

implying the second expression for $\Delta_{j}$. $\quad\Box$

###### Lemma 3.2

For $n\geq 1$,

$\|\psi^{n}\|\leq\int_{I_{1}}t\|u^{\prime\prime}(t)\|\,dt+\frac{1}{12}\biggl{(}\tau_{2}^{2}\|u^{\prime\prime}(t_{1})\|+2\tau_{n}^{2}\|u^{\prime\prime}(t_{n})\|+3\sum_{j=2}^{n}\tau_{j}^{2}\int_{I_{j}}\|u^{\prime\prime\prime}\|\,dt\biggr{)}$

and
$\|\psi(t)\|\leq\|u(t)-u_{I}(t)\|+\max\bigl{(}\|\psi^{n}\|,\|\psi^{n-1}\|\bigr{)}$
for $t\in I_{n}$.

Proof. For $t\in I_{j}$, we have the identity

$-u^{\prime\prime}(t)=\frac{1}{2}\int_{t}^{t_{j}}u^{\prime\prime\prime}(s)\,ds-\frac{1}{2}\int_{t_{j-1}}^{t}u^{\prime\prime\prime}(s)\,ds-(u^{\prime\prime}(t_{j})+u^{\prime\prime}(t_{j-1}))/2.$

Multiply both sides by $(t_{j}-t)(t-t_{j-1})$ and integrate to obtain

$\Delta_{j}=-\frac{\tau_{j}^{2}}{12}[u^{\prime\prime}(t_{j})+u^{\prime\prime}(t_{j-1})]+R_{j},$

where

$R_{j}=\frac{1}{2\tau_{j}}\int_{I_{j}}(t_{j}-t)(t-t_{j-1})\biggl{(}\int_{t}^{t_{j}}u^{\prime\prime\prime}(s)\,ds-\int_{t_{j-1}}^{t}u^{\prime\prime\prime}(s)\,ds\biggr{)}\,dt.$

Thus, by Lemma 3.1,

$(-1)^{n}\psi^{n}=\Delta_{1}+\frac{1}{12}\sum_{j=2}^{n}(-1)^{j}\tau_{j}^{2}[u^{\prime\prime}(t_{j})+u^{\prime\prime}(t_{j-1})]-\sum_{j=2}^{n}(-1)^{j}R_{j}.$

Shifting the summation index, so

$\sum_{j=2}^{n}(-1)^{j}\tau_{j}^{2}u^{\prime\prime}(t_{j-1})=\tau_{2}^{2}u^{\prime\prime}(t_{1})-\sum_{j=2}^{n-1}(-1)^{j}\tau_{j+1}^{2}u^{\prime\prime}(t_{j}).$

Since $\Delta_{1}=-\tau_{1}^{-1}\int_{I_{1}}(t_{1}-s)su^{\prime\prime}(s)\,ds$
and since
$\|R_{j}\|\leq\frac{\tau_{j}^{2}}{6}\int_{I_{j}}\|u^{\prime\prime\prime}(t)\|\,dt,$

$\displaystyle\|\psi^{n}\|$
$\displaystyle\leq\int_{I_{1}}t\|u^{\prime\prime}(t)\|\,dt+\frac{1}{12}\biggl{(}\tau_{2}^{2}\|u^{\prime\prime}(t_{1})\|+\tau_{n}^{2}\|u^{\prime\prime}(t_{n})\|$
$\displaystyle\qquad{}+\sum_{j=2}^{n-1}(\tau_{j+1}^{2}-\tau_{j}^{2})\|u^{\prime\prime}(t_{j})\|+2\sum_{j=2}^{n}\tau_{j}^{2}\int_{I_{j}}\|u^{\prime\prime\prime}(t)\|\,dt\biggr{)}.$

Using

$\sum_{j=2}^{n-1}\tau_{j+1}^{2}\|u^{\prime\prime}(t_{j})\|=\sum_{j=3}^{n}\tau_{j}^{2}\|u^{\prime\prime}(t_{j-1})\|\leq\sum_{j=3}^{n}\tau_{j}^{2}\Big{(}\|u^{\prime\prime}(t_{j})\|+\int_{I_{j}}\|u^{\prime\prime\prime}(t)\|\,dt\Big{)},$

and the bound for $\|\psi^{n}\|$ follows after canceling the common terms.

The interpolant $\psi_{I}$, defined as in (2.7), satisfies
$\psi_{I}=u_{I}-\widehat{u}$, leading to the representation

$\psi=u-\widehat{u}=\psi_{I}+u-u_{I},$ (3.1)

which implies the bound for $\|\psi\|$ because
$\|\psi_{I}\|\leq\max\bigl{(}\|\psi^{n}\|,\|\psi^{n-1}\|\bigr{)}$ on $I_{n}$.
$\quad\Box$

###### Corollary 3.1

Under the regularity assumption in (2.11) and for a time mesh of the form
(2.1) with grading parameter $\gamma\geq 1$ we have, for $n\geq 1$,

$\|\psi^{n}\|\leq\|\psi\|_{I_{n}}\leq
C\times\begin{cases}\tau^{2}\log(t_{n}/t_{1}),&\text{if $\gamma=2/\sigma$,}\\\
\tau^{\min(\gamma\sigma,2)}t_{n}^{\max(0,\delta)},&\text{if $\gamma\neq
2/\sigma$.}\end{cases}$

Proof. First we show that for $\gamma\geq 1,$

$\|u-u_{I}\|_{I_{n}}\leq C\tau^{\min(\gamma\sigma,2)}t_{n}^{\max(0,\delta)}.$
(3.2)

Since the interpolation error $u-u_{I}$ vanishes if $u$ is a polynomial of
degree $1$, by computing the Peano kernel one finds that for $t\in I_{n}$,

$u(t)-u_{I}(t)=-\frac{t_{n}-t}{\tau_{n}}\int_{t_{n-1}}^{t}(s-t_{n-1})u^{\prime\prime}(s)\,ds-\frac{t-t_{n-1}}{\tau_{n}}\int_{t}^{t_{n}}(t_{n}-s)u^{\prime\prime}(s)\,ds.$
(3.3)

The right-side is bounded by

$\int_{t_{n-1}}^{t}(s-t_{n-1})\|u^{\prime\prime}(s)\|\,ds+\int_{t}^{t_{n}}(t-t_{n-1})\|u^{\prime\prime}(s)\|\,ds\leq\int_{I_{n}}(s-t_{n-1})\|u^{\prime\prime}(s)\|\,ds$

and so, by using the time mesh properties in (2.2), we get

$\|u-u_{I}\|_{I_{n}}\leq
C\tau_{n}t_{n}^{-1}\int_{I_{n}}t\|u^{\prime\prime}(t)\|\,dt\leq
C\tau_{n}t_{n}^{-1}\int_{I_{n}}t^{\sigma-1}\,dt\leq
C\tau_{n}^{2}t_{n}^{\sigma-2},~{}~{}n\geq 1.$

Since $\tau_{n}\leq C\tau t_{n}^{1-1/\gamma}$, the proof of (3.2) is completed
after noting that

$\tau_{n}^{2}t_{n}^{\sigma-2}\leq C\tau^{2}t_{n}^{\delta}\leq
C\tau^{2}\max(t_{n}^{\delta},t_{1}^{\delta})\leq
C\tau^{\min(\gamma\sigma,2)}t_{n}^{\max(0,\delta)}.$ (3.4)

Turning to the estimate for $\psi^{n}$, Lemma 3.2 and (2.11) imply that

$\|\psi^{n}\|\leq
C\int_{I_{1}}t^{\sigma-1}\,dt+C\tau_{2}^{2}t_{1}^{\sigma-2}+C\tau_{n}^{2}t_{n}^{\sigma-2}+C\sum_{j=2}^{n}\tau_{j}^{2}\int_{I_{j}}t^{\sigma-3}\,dt,$

for $1\leq n\leq N.$ Since $t_{1}=\tau^{\gamma}$ and $\tau_{2}\leq
2^{\gamma}\tau^{\gamma},$
$\int_{I_{1}}t^{\sigma-1}\,dt+\tau_{2}^{2}t_{1}^{\sigma-2}\leq
C\tau^{\gamma\sigma},$ and we again bound $\tau_{n}^{2}t_{n}^{\sigma-2}$ using
(3.4). For the sum over $j$,

$\sum_{j=2}^{n}\tau_{j}^{2}\int_{I_{j}}t^{\sigma-3}\,dt\leq
C\sum_{j=2}^{n}\tau^{2}t_{j}^{2-2/\gamma}\int_{I_{j}}t^{\sigma-3}\,dt\leq
C\tau^{2}\int_{t_{1}}^{t_{n}}t^{\delta-1}\,dt,$

and the estimate for $\|\psi^{n}\|$ follows. $\quad\Box$

The next target is to estimate
$\mathcal{I}(\langle\mathcal{I}^{1-\alpha}\psi^{\prime},\psi^{\prime}\rangle)(t_{n})$.
Preceding this, we need to bound $\|\psi^{\prime}\|$ in the next lemma.

###### Lemma 3.3

We have $\|\psi^{\prime}(t)\|\leq t^{\sigma-1}$ for $t\in I_{1}$. Moreover,
$\|\psi^{\prime}(t)\|\leq C\tau^{2}\tau_{n}^{-1}t_{n}^{\delta}$ for $t\in
I_{n}$ with $n\geq 2,$ and for $\delta>0$.

Proof. Differentiating (3.1) and (3.3), we see for $t\in I_{n}$ that

$\psi^{\prime}(t)=\frac{1}{\tau_{n}}\int_{t_{n-1}}^{t}(s-t_{n-1})u^{\prime\prime}(s)\,ds-\frac{1}{\tau_{n}}\int_{t}^{t_{n}}(t_{n}-s)u^{\prime\prime}(s)\,ds+\tau_{n}^{-1}(\psi^{n}-\psi^{n-1}).$

Thus, if $t\in I_{1}$ then, by (2.11), Corollary 3.1, and the fact that
$\psi^{0}=0,$

$\|\psi^{\prime}(t)\|\leq
C\tau^{-1}\int_{0}^{t}s^{\sigma-1}\,ds+C\tau_{1}^{-1}\int_{t}^{t_{1}}(t_{1}-s)s^{\sigma-2}\,ds+C\tau_{1}^{\sigma-1}\leq
Ct^{\sigma-1}\,.$

If $\delta>0$, $n\geq 2$ and $t\in I_{n}$ then, recalling (3.4) and using
again (2.11),

$\tau_{n}\|\psi^{\prime}(t)\|\leq
C\tau_{n}\int_{I_{n}}t^{\sigma-2}\,dt+C\tau^{2}t_{n}^{\delta}\leq
C\tau_{n}^{2}t_{n}^{\sigma-2}+C\tau^{2}t_{n}^{\delta}\leq
C\tau^{2}t_{n}^{\delta},$

showing that $\|\psi^{\prime}(t)\|\leq C\tau^{2}\tau_{n}^{-1}t_{n}^{\delta}$.
$\quad\Box$

###### Lemma 3.4

Assume $\sigma>\alpha/2$. Then
$\big{|}\mathcal{I}(\langle\mathcal{I}^{1-\alpha}\psi^{\prime},\psi^{\prime}\rangle)(t_{n})\big{|}\leq
C\tau^{3-\alpha}$ for $n\geq 1$, provided that
$\gamma>\max\\{2/\sigma,(3-\alpha)/(2\sigma-\alpha)\\}$.

Proof. For $n=1$, the Cauchy-Schwarz inequality, Lemma 3.3, and the assumption
$\sigma>\alpha/2$ give

$\Big{|}\mathcal{I}(\langle\mathcal{I}^{1-\alpha}\psi^{\prime},\psi^{\prime}\rangle)(t_{1})\Big{|}\leq
C\int_{0}^{t_{1}}t^{\sigma-1}\int_{0}^{t}(t-s)^{-\alpha}s^{\sigma-1}\,ds\,dt\leq
C\tau_{1}^{2\sigma-\alpha}.$ (3.5)

To deal with the case $n\geq 2$, we make the splitting
$\mathcal{I}^{1-\alpha}\psi^{\prime}=T_{1}+T_{2}$ where

$T_{1}(t)=\int_{0}^{t_{j-1}}\omega_{1-\alpha}(t-s)\psi^{\prime}(s)\,ds\quad\text{and}\quad
T_{2}(t)=\int_{t_{j-1}}^{t}\omega_{1-\alpha}(t-s)\psi^{\prime}(s)\,ds$

for $t\in I_{j}$ and $j\geq 2$. Using Lemma 3.3 and (2.2), we observe that

$\|T_{2}(t)\|\leq
C\int_{t_{j-1}}^{t}\omega_{1-\alpha}(t-s)\,\frac{\tau^{2}s^{\delta}}{\tau_{j}}\,ds\leq
C\frac{\tau^{2}t^{\delta}}{\tau_{j}}\omega_{2-\alpha}(t-t_{j-1})\leq
C\tau^{2}t^{\delta}\tau_{j}^{-\alpha}.$

For estimating $T_{1}(t)$, integrate by parts recalling $\psi^{0}=0$,

$T_{1}(t)=\omega_{1-\alpha}(t-t_{j-1})\psi^{j-1}+\int_{0}^{t_{j-1}}\omega_{-\alpha}(t-s)\psi(s)\,ds,$

where $\omega_{-\alpha}(t)=\omega_{1-\alpha}^{\prime}(t)=-\alpha
t^{-\alpha-1}/\Gamma(1-\alpha)$. We apply Corollary 3.1 to conclude that

$\|T_{1}(t)\|\leq
C\tau^{2}t_{j-1}^{\delta}\biggl{(}\omega_{1-\alpha}(t-t_{j-1})-\int_{0}^{t_{j-1}}\omega_{-\alpha}(t-s)\,ds\biggr{)}\leq
C\tau^{2}t_{j}^{\delta}\omega_{1-\alpha}(t-t_{j-1}).$

Lemma 3.3 and above estimates for $\|T_{1}\|$ and $\|T_{2}\|$ yield

$\|\psi^{\prime}\|_{I_{j}}\int_{I_{j}}\bigl{(}\|T_{1}\|+\|T_{2}\|\bigr{)}\,dt\leq
C(\tau^{2}t_{j}^{\delta}\tau_{j}^{-1})\tau^{2}t_{j}^{\delta}(\tau_{j}^{1-\alpha})\leq\tau^{4}t_{j}^{2\delta}\tau_{j}^{-\alpha},~{}~{}{\rm
for}~{}~{}j\geq 2\,.$

By using this and (3.5), and noting that $\gamma(2\sigma-\alpha)>3-\alpha$, we
reach

$\displaystyle\Big{|}\mathcal{I}($
$\displaystyle\langle\mathcal{I}^{1-\alpha}\psi^{\prime},\psi^{\prime}\rangle)(t_{n})\Big{|}\leq\biggl{|}\int_{0}^{t_{1}}\langle\mathcal{I}^{1-\alpha}\psi^{\prime},\psi^{\prime}\rangle\,dt\biggr{|}+\biggl{|}\int_{t_{1}}^{t_{n}}\langle\mathcal{I}^{1-\alpha}\psi^{\prime},\psi^{\prime}\rangle\,dt\biggr{|}$
$\displaystyle\leq
C\tau^{3-\alpha}+\sum_{j=2}^{n}\|\psi^{\prime}\|_{I_{j}}\int_{I_{j}}\bigl{(}\|T_{1}\|+\|T_{2}\|\bigr{)}\,dt\leq
C\tau^{3-\alpha}+C\sum_{j=2}^{n}\tau^{4}t_{j}^{2\delta}\tau_{j}^{-\alpha},$

for $j\geq 2.$ By (2.2), if $j\geq 2$ then $\tau^{1+\alpha}\leq
C\tau_{j}^{1+\alpha}t_{j}^{-(1+\alpha)(1-1/\gamma)}$ so

$\tau^{1+\alpha}\sum_{j=2}^{n}t_{j}^{2\delta}\tau_{j}^{-\alpha}\leq
C\sum_{j=2}^{n}t_{j}^{2\delta-(1+\alpha)(1-1/\gamma)}\tau_{j}\leq
C\int_{t_{1}}^{t_{n}}t^{2\sigma-\alpha-(3-\alpha)/\gamma-1}\,dt\leq C,$

and therefore the desired bound holds. $\quad\Box$

##  4 Errors from the time discretizations

This section is devoted to estimating the error $\eta=u-U$ from the time
discretization in the norm of $L^{2}(;L^{2}(\Omega))$. To achieve an optimal
convergence rate, we employ a duality argument in addition to the usage of the
time graded meshes. By reversing the order of integration, we find that

$\mathcal{I}(\langle\mathcal{I}^{\alpha}v,w\rangle)(T)=\mathcal{I}(\langle
v,\mathcal{J}_{T}^{\alpha}w\rangle)(T)\quad\text{where}\quad(\mathcal{J}_{T}^{\alpha}w)(t)=\int_{t}^{T}\omega_{\alpha}(s-t)w(t)\,ds.$

Using
$\mathcal{I}(\langle\partial_{t}^{\alpha}v,w\rangle)(T)=\mathcal{I}(\langle
v^{\prime},\mathcal{J}_{T}^{1-\alpha}w\rangle)(T)$ and integrating by parts
yield

$\mathcal{I}(\langle\partial_{t}^{\alpha}v,w\rangle)(T)=-\langle
v(0),(\mathcal{J}_{T}^{1-\alpha}w)(0)\rangle-\mathcal{I}(\langle
v,(\mathcal{J}_{T}^{1-\alpha}w)^{\prime}\rangle)(T),$ (4.1)

and since
$\partial_{t}^{\alpha}v=\mathcal{I}^{1-\alpha}v^{\prime}=(\mathcal{I}^{1-\alpha}v)^{\prime}-v(0)\omega_{1-\alpha}$,

$\mathcal{I}(\langle(\mathcal{I}^{1-\alpha}v)^{\prime},w\rangle)(T)=-\mathcal{I}(\langle
v,(\mathcal{J}_{T}^{1-\alpha}w)^{\prime}\rangle)(T).$ (4.2)

We remark that Zhang et al. [35, Equation (89)] have recently exploited this
dual operator $w\mapsto-(\mathcal{J}_{T}^{1-\alpha}w)^{\prime}$ in the error
analysis of a discontinuous Galerkin scheme for (1.1).

Suppose that $\varphi$ satisfies the final-value problem

$-(\mathcal{J}_{T}^{1-\alpha}\varphi)^{\prime}+\mathcal{A}\varphi=\eta(t)\quad\text{for
$0<t<T$,}\quad\text{with $\varphi(T)=0$,}$ (4.3)

subject to homogeneous Dirichlet boundary conditions. Let
$y(t)=\varphi(0)+\int_{0}^{t}\varphi(s)\,ds$ so that $y$ solves the initial-
value problem

$y^{\prime}=\varphi\quad\text{for~{}$0<t<T$,}\quad\text{with
$y(0)=\varphi(0)$,}$ (4.4)

and with $y_{I}$ denotes the continuous piecewise-linear function that
interpolates $y$ at the time levels $t_{j}$, put

$Y=y-y_{I}.$ (4.5)

###### Lemma 4.1

With the notation above, $\|\eta\|^{2}_{L^{2}(J)}\leq\mathcal{I}(\langle
Y^{\prime},\mathcal{I}^{1-\alpha}\eta^{\prime}+\mathcal{A}\eta\rangle)(T).$

Proof. Using (4.1), (4.3), $\eta(0)=0$ and (4.4),

$\|\eta\|^{2}_{L^{2}(J)}=\mathcal{I}(\langle-(\mathcal{J}_{T}^{1-\alpha}\varphi)^{\prime}+\mathcal{A}\varphi,\eta\rangle)(T)=\mathcal{I}(\langle
y^{\prime},\mathcal{I}^{1-\alpha}\eta^{\prime}+\mathcal{A}\eta\rangle)(T).$

At the same time, (2.3) and (2.4) imply that

$\frac{1}{\tau_{n}}\int_{I_{n}}\bigl{(}\partial_{t}^{\alpha}\eta+\mathcal{A}\eta\bigr{)}\,dt=0\quad\text{so}\quad\int_{I_{n}}\langle
y_{I}^{\prime},\partial_{t}^{\alpha}\eta+\mathcal{A}\eta\rangle\,dt=0,$

because $y_{I}^{\prime}$ is constant on $I_{n}$. Since
$Y^{\prime}=y^{\prime}-y_{I}^{\prime}$, the inequality follows at once.
$\quad\Box$

We will show below in Theorem 4.2 that the interpolation error $Y$ satisfies

$\mathcal{I}(\langle
Y^{\prime},\mathcal{I}^{1-\alpha}Y^{\prime}\rangle)(T)+\mathcal{I}(\langle\mathcal{A}Y,\mathcal{I}^{\alpha}\mathcal{A}Y^{\prime}\rangle)(T)\leq
C\tau^{1+\alpha}\|\eta\|^{2}_{L^{2}(J)}.$ (4.6)

Assuming this fact for now, we can derive an estimate for $\eta$ in terms of
$\psi^{\prime}$. For convenience, we use the notations:
$F(g)=\bigl{(}\mathcal{I}(\langle\mathcal{I}^{\alpha}g^{\prime},g\rangle)(T)\bigr{)}^{1/2}$
and
$G(g)=\bigl{(}\mathcal{I}(\langle\mathcal{I}^{1-\alpha}g^{\prime},g^{\prime}\rangle)(T)\bigr{)}^{1/2}$.

###### Theorem 4.1

We have $\alpha^{2}\|\eta\|^{2}_{L^{2}(J)}\leq
C\,\tau^{\alpha+1}\mathcal{I}(\langle\mathcal{I}^{1-\alpha}\psi^{\prime},\psi^{\prime}\rangle)(T).$

Proof. It suffices to estimate the right-hand side of the inequality in Lemma
4.1. By Lemma 2.1,

$\alpha\mathcal{I}(\langle
Y^{\prime},\mathcal{I}^{1-\alpha}\psi^{\prime}\rangle)(T)\leq
G(Y)G(\psi)~{}~{}{\rm and}~{}~{}\alpha\mathcal{I}(\langle
Y^{\prime},\mathcal{I}^{1-\alpha}\theta^{\prime}\rangle)(T)\leq
G(Y)G(\theta)\,.$

After using (2.10) and $\eta^{\prime}=\psi^{\prime}-\theta^{\prime}$, we
conclude that

$\alpha^{2}\mathcal{I}(\langle
Y^{\prime},\mathcal{I}^{1-\alpha}\eta^{\prime}\rangle)(T)\leq 2G(Y)G(\psi).$
(4.7)

Since $Y(t_{j})=0$ for $0\leq j\leq N$, integrating by parts, using
$\mathcal{A}Y=\mathcal{I}\mathcal{A}Y^{\prime}=\mathcal{I}^{1-\alpha}(\mathcal{I}^{\alpha}\mathcal{A}Y^{\prime})$,
and applying Lemma 2.1,

$\displaystyle\alpha\mathcal{I}(\langle
Y^{\prime},\mathcal{A}\theta\rangle)(T)$
$\displaystyle=-\alpha\mathcal{I}(\langle\theta^{\prime},\mathcal{I}^{1-\alpha}(\mathcal{I}^{\alpha}\mathcal{A}Y^{\prime})\rangle)(T)\leq
G(\theta)F({\mathcal{A}}Y).$

The same estimate holds with $\theta$ replaced by $\psi$, so because
$\eta=\psi-\theta$, and using (2.10) again,

$\alpha^{2}\mathcal{I}(\langle Y^{\prime},\mathcal{A}\eta\rangle)(T)\leq
2F({\mathcal{A}}Y)G(\psi).$ (4.8)

Adding (4.7) and (4.8), we see $\alpha^{2}\|\eta\|^{2}_{L^{2}(J)}\leq
2G(\psi)(G(Y)+F({\mathcal{A}}Y))$ by Lemma 4.1. Squaring both sides, we have

$\alpha^{4}\|\eta\|^{4}_{L^{2}(J)}\leq
4(G(\psi))^{2}(G(Y)+F({\mathcal{A}}Y))^{2}\leq
8(G(\psi))^{2}((G(Y))^{2}+(F({\mathcal{A}}Y))^{2}).$

Since $(G(Y))^{2}+(F({\mathcal{A}}Y))^{2}$ is just the left-hand side of
(4.6), the desired inequality followed after cancelling the common factor
$\|\eta\|^{2}_{L^{2}(J)}$. $\quad\Box$

It remains to prove (4.6). We start by showing preliminary bounds for
$\big{|}\mathcal{I}(\langle
Y^{\prime},\mathcal{I}^{1-\alpha}Y^{\prime}\rangle)(T)\big{|}$ and
$\big{|}\mathcal{I}(\langle\mathcal{A}Y,\mathcal{I}^{1-\alpha}\mathcal{A}Y^{\prime}\rangle)(T)\big{|}$
in the next two lemmas.

###### Lemma 4.2

$\Big{|}\mathcal{I}(\langle
Y^{\prime},\mathcal{I}^{1-\alpha}Y^{\prime}\rangle)(T)\Big{|}\leq
C(1-\alpha)\sum_{j=1}^{N-2}\tau_{j}^{-\alpha}\|Y\|_{I_{j}}^{2}+C\tau^{1-\alpha}\|Y^{\prime}\|^{2}_{L^{2}(J)}.$

Proof. For $t\in I_{j}$ with $j\geq 2$, we write
$\mathcal{I}^{1-\alpha}Y^{\prime}(t)=S_{1}(t)+S_{2}(t)$ where

$S_{1}(t)=\int_{t_{j-2}}^{t}\omega_{1-\alpha}(t-s)Y^{\prime}(s)\,ds~{}\text{and}~{}S_{2}(t)=\int_{0}^{t_{j-2}}\omega_{1-\alpha}(t-s)Y^{\prime}(s)\,ds.$

Applying the Cauchy-Schwarz inequality, integrating, changing the order of
integration, and integrating again, yields

$\displaystyle\|S_{1}(t)$
$\displaystyle\|^{2}_{L^{2}(I_{j})}\leq\int_{I_{j}}\biggl{(}\int_{t_{j-2}}^{t}\omega_{1-\alpha}(t-s)\,ds\biggr{)}\biggl{(}\int_{t_{j-2}}^{t}\omega_{1-\alpha}(t-s)\|Y^{\prime}(s)\|^{2}\,ds\biggr{)}\,dt$
$\displaystyle\leq
C(\tau_{j}+\tau_{j-1})^{1-\alpha}\biggl{(}\int_{I_{j}}\int_{s}^{t_{j}}+\int_{I_{j-1}}\int_{I_{j}}\biggr{)}\omega_{1-\alpha}(t-s)\,dt\,\|Y^{\prime}(s)\|^{2}\,ds$
$\displaystyle\leq
C(\tau_{j}+\tau_{j-1})^{2(1-\alpha)}\int_{t_{j-2}}^{t_{j}}\|Y^{\prime}(s)\|^{2}\,ds\leq
C\tau_{j}^{2(1-\alpha)}\int_{t_{j-2}}^{t_{j}}\|Y^{\prime}(s)\|^{2}\,ds.$

For $t\in I_{1}$, let
$S_{1}(t)=\int_{0}^{t}\omega_{1-\alpha}(t-s)Y^{\prime}(s)\,ds$ and
$S_{2}(t)=0$. Following the steps above, we find that
$\|S_{1}\|^{2}_{L^{2}(I_{1})}\leq
C\tau_{1}^{2(1-\alpha)}\|Y^{\prime}\|^{2}_{L^{2}(I_{1})}$. Thus,
$\|S_{1}\|^{2}_{L^{2}(J)}\leq
C\tau^{2(1-\alpha)}\|Y^{\prime}\|^{2}_{L^{2}(J)},$ and consequently

$\mathcal{I}(|\langle
Y^{\prime},S_{1}\rangle|)(T)\leq\|Y^{\prime}\|_{L^{2}(J)}\|S_{1}\|_{L^{2}(J)}\leq
C\tau^{1-\alpha}\|Y^{\prime}\|^{2}_{L^{2}(J)}.$ (4.9)

Turning to the second term $S_{2}(t)$, we integrate by parts (noting that
$Y^{j-2}=0=Y^{0}$) and obtain

$\|S_{2}(t)\|=\biggl{\|}\int_{0}^{t_{j-2}}\omega_{-\alpha}(t-s)Y(s)\,ds\biggr{\|}\leq\sum_{i=1}^{j-2}\|Y\|_{I_{i}}\int_{I_{i}}|\omega_{-\alpha}(t-s)|\,ds.$

Since $|\omega_{-\alpha}(t-s)|\leq|\omega_{-\alpha}(t_{j-1}-s)|$ for $t\in
I_{j}$,

$\int_{I_{j}}|\langle
Y^{\prime},S_{2}\rangle|\,dt\leq\sum_{i=1}^{j-2}\int_{I_{j}}\|Y^{\prime}(t)\|\|Y\|_{I_{i}}\,dt\int_{I_{i}}|\omega_{-\alpha}(t_{j-1}-s)|\,ds,$

with

$\int_{I_{j}}\|Y^{\prime}(t)\|\|Y\|_{I_{i}}\,dt\leq\frac{1}{2}\biggl{(}\|Y\|_{I_{i}}^{2}+\tau_{j}\|Y^{\prime}\|^{2}_{L^{2}(I_{j})}\biggr{)},$

and, remembering that
$\omega_{-\alpha}(t)=-\alpha\,t^{-1}\omega_{1-\alpha}(t)$,

$\int_{I_{i}}|\omega_{-\alpha}(t_{j-1}-s)|\,ds=\omega_{1-\alpha}(t_{j-1}-t_{i})-\omega_{1-\alpha}(t_{j-1}-t_{i-1}).$

Since $t_{j}-t_{i}\geq t_{j-1}-t_{i-1}$,
$\omega_{1-\alpha}(t_{j-1}-t_{i-1})\geq\omega_{1-\alpha}(t_{j}-t_{i})$, and
thus,

$\int_{I_{j}}|\langle
Y^{\prime},S_{2}\rangle|\,dt\leq\frac{1}{2}\sum_{i=1}^{j-2}\|Y\|_{i}^{2}[\omega_{1-\alpha}(t_{j-1}-t_{i})-\omega_{1-\alpha}(t_{j}-t_{i})]\\\
+\frac{\tau_{j}}{2}\|Y^{\prime}\|^{2}_{L^{2}(I_{j})}\int_{0}^{t_{j-2}}|\omega_{-\alpha}(t_{j-1}-s)|\,ds,$

and in the second term, noting that
$1/\Gamma(1-\alpha)=(1-\alpha)/\Gamma(2-\alpha)$,

$\int_{0}^{t_{j-2}}|\omega_{-\alpha}(t_{j-1}-s)|\,ds=\omega_{1-\alpha}(t_{j-1}-t_{j-2})-\omega_{1-\alpha}(t_{j})\leq
C(1-\alpha)\tau_{j-1}^{-\alpha}.$

By interchanging the order of the double sum,

$\displaystyle\mathcal{I}(|\langle
Y^{\prime},S_{2}\rangle|)(T)=\sum_{j=2}^{N}\int_{I_{j}}|\langle
Y^{\prime},S_{2}\rangle|\,dt$ $\displaystyle\leq
C(1-\alpha)\Big{(}\sum_{i=1}^{N-2}\|Y\|_{I_{i}}^{2}\sum_{j=i+2}^{N}[(t_{j-1}-t_{i})^{-\alpha}-(t_{j}-t_{i})^{-\alpha}]+\Bigl{(}\max_{2\leq
j\leq
N}\tau_{j}\tau_{j-1}^{-\alpha}\Bigr{)}\|Y^{\prime}\|^{2}_{L^{2}(J)}\Big{)}$
$\displaystyle\leq
C(1-\alpha)\biggl{(}\sum_{i=1}^{N-2}\tau_{i+1}^{-\alpha}\|Y\|_{I_{i}}^{2}+\tau^{1-\alpha}\|Y^{\prime}\|^{2}_{L^{2}(J)}\biggr{)},$

which, combined with (4.9), yields the desired estimate. $\quad\Box$

###### Lemma 4.3

$\Big{|}\mathcal{I}(\langle\mathcal{A}Y,\mathcal{I}^{\alpha}\mathcal{A}Y^{\prime}\rangle)(T)\Big{|}\leq
C\sum_{j=1}^{N}\tau_{j}^{\alpha}\biggr{[}\biggl{(}\int_{I_{j}}\|\mathcal{A}Y^{\prime}\|\,dt\biggr{)}^{2}+\|\mathcal{A}Y\|_{I_{j}}^{2}\biggr{]}.$

Proof. Integrating by parts,

$\mathcal{I}(\langle\mathcal{A}Y,\mathcal{I}^{\alpha}\mathcal{A}Y^{\prime}\rangle)(T)=-\int_{0}^{T}\langle\mathcal{A}Y^{\prime},S_{3}\rangle\,dt-\int_{0}^{T}\langle\mathcal{A}Y^{\prime},S_{4}\rangle\,dt,$
(4.10)

where we used the splitting $\mathcal{I}^{\alpha}{\mathcal{A}}Y=S_{3}+S_{4}$
with

$S_{3}(t)=\int_{t_{j-1}}^{t}\omega_{\alpha}(t-s)\mathcal{A}Y(s)\,ds~{}~{}\text{and}~{}~{}S_{4}(t)=\int_{0}^{t_{j-1}}\omega_{\alpha}(t-s)\mathcal{A}Y(s)\,ds$

for $t\in I_{j}$. Since
$\|S_{3}(t)\|\leq\|\mathcal{A}Y\|_{I_{j}}\int_{t_{j-1}}^{t}\omega_{\alpha}(t-s)\,ds\leq
C\tau_{j}^{\alpha}\|\mathcal{A}Y\|_{I_{j}}$,

$\biggl{|}\int_{0}^{T}\langle\mathcal{A}Y^{\prime},S_{3}\rangle\,dt\biggr{|}\leq
C\sum_{j=1}^{N}\tau_{j}^{\alpha}\|\mathcal{A}Y\|_{I_{j}}\int_{I_{j}}\|\mathcal{A}Y^{\prime}(t)\|\,dt.$
(4.11)

For the estimate involving $S_{4}$, we reverse the order of integration and
then integrate by parts, to obtain

$\displaystyle\int_{I_{j}}\langle\mathcal{A}Y^{\prime},S_{4}\rangle\,dt$
$\displaystyle=-\int_{0}^{t_{j-1}}\Big{\langle}\mathcal{A}Y(s),\int_{I_{j}}\omega_{\alpha-1}(t-s)\mathcal{A}Y(t)\,dt\Big{\rangle}\,ds,$

and thus, applying the Cauchy-Schwarz inequality and using
$\|\mathcal{A}Y\|_{I_{i}}\|\mathcal{A}Y\|_{I_{j}}\leq\|\mathcal{A}Y\|_{I_{i}}^{2}+\|\mathcal{A}Y\|_{I_{j}}^{2},$
we get

$\biggl{|}\int_{0}^{T}\langle\mathcal{A}Y^{\prime},S_{4}\rangle\,dt\biggr{|}\leq\sum_{j=1}^{N}\sum_{i=1}^{j-1}\bigl{(}\|{\mathcal{A}}Y\|_{I_{i}}^{2}+\|\mathcal{A}Y\|_{I_{j}}^{2}\bigr{)}\int_{I_{i}}\int_{I_{j}}|\omega_{\alpha-1}(t-s)|\,dt\,ds\\\
=\sum_{i=1}^{N-1}\|\mathcal{A}Y\|_{I_{i}}^{2}\int_{I_{i}}\int_{t_{i}}^{t_{n}}|\omega_{\alpha-1}(t-s)|\,dt\,ds+\sum_{j=1}^{N}\|\mathcal{A}Y\|_{I_{j}}^{2}\int_{I_{j}}\int_{0}^{t_{j-1}}|\omega_{\alpha-1}(t-s)|\,ds\,dt\\\
\leq\sum_{i=1}^{N-1}\|{\mathcal{A}}Y\|_{I_{i}}^{2}\int_{I_{i}}\omega_{\alpha}(t_{i}-s)\,ds+\sum_{j=1}^{N}\|{\mathcal{A}}Y\|_{I_{j}}^{2}\int_{I_{j}}\omega_{\alpha}(t-t_{j-1})\,dt\leq
C\sum_{j=1}^{N}\tau_{j}^{\alpha}\|{\mathcal{A}}Y\|_{I_{j}}^{2}.$

The proof is concluded by inserting this and (4.11) into the splitting (4.10).
$\quad\Box$

By using the achieved estimates in the previous two lemmas, we are now able in
the next theorem to provide the missing part in the proof of Theorem 4.1.

###### Theorem 4.2

The inequality (4.6) is satisfied by the function $Y$ defined via (4.3), (4.4)
and (4.5).

Proof. Recall that $Y=y-y_{I}$ where $y_{I}$ is the piecewise linear
polynomial that interpolates $y$ at the time nodes, and $y^{\prime}=\varphi$.
Thus, if $t\in I_{j}$ then

$Y(t)=y(t)-y_{I}(t)=\frac{t_{j}-t}{\tau_{j}}\int_{t_{j-1}}^{t}\varphi(s)\,ds-\frac{t-t_{j-1}}{\tau_{j}}\int_{t}^{t_{j}}\varphi(s)\,ds$
(4.12)

so $\|Y(t)\|\leq\int_{I_{j}}\|\varphi\|\,ds$. Similarly, replacing $u$ with
$y$ in (3.3), we have
$\|Y(t)\|\leq\tau_{j}\int_{I_{j}}\|\varphi^{\prime}\|\,ds$. Therefore

$\|Y\|_{I_{j}}^{2}\leq\tau_{j}\|\varphi\|^{2}_{L^{2}(I_{j})}\quad\text{and}\quad\|Y\|_{I_{j}}^{2}\leq\tau_{j}^{3}\|\varphi^{\prime}\|^{2}_{L^{2}(I_{j})}.$
(4.13)

Consider the linear operator $B$ defined by
$(B\varphi)(t)=\tau_{j}^{-\alpha/2}\|Y\|_{I_{j}}$ for $t\in I_{j}$ and $1\leq
j\leq N$. The estimates (4.13) give

$\|B\varphi\|^{2}_{L^{2}(J)}\leq
C\tau^{1-\alpha}\|\varphi\|^{2}_{L^{2}(J)}\quad\text{and}\quad\|B\varphi\|^{2}_{L^{2}(J)}\leq
C\tau^{3-\alpha}\|\varphi^{\prime}\|^{2}_{L^{2}(J)},$

and, since $\varphi(T)=0$ and
$(\tau^{1-\alpha})^{1-\alpha}(\tau^{3-\alpha})^{\alpha}=\tau^{1+\alpha}$, we
may apply Corollary 7.1 from the Appendix to deduce that

$\sum_{j=1}^{N}\tau_{j}^{-\alpha}\|Y\|_{I_{j}}^{2}=\|B\varphi\|^{2}_{L^{2}(J)}\leq
C\tau^{1+\alpha}\|(\mathcal{J}_{T}^{1-\alpha}\varphi)^{\prime}\|^{2}_{L^{2}(J)}.$
(4.14)

Furthermore, by differentiating (4.12) and (3.3),

$Y^{\prime}(t)=\varphi(t)-\frac{1}{\tau_{j}}\int_{I_{j}}\varphi(s)\,ds=\frac{1}{\tau_{j}}\int_{t_{j-1}}^{t}(s-t_{j-1})\varphi^{\prime}(s)\,ds-\frac{1}{\tau_{j}}\int_{t}^{t_{j}}(t_{j}-s)\varphi^{\prime}(s)\,ds$

for $t\in I_{j}$, so

$\|Y^{\prime}(t)\|\leq\|\varphi(t)\|+\frac{1}{\tau_{j}}\int_{I_{j}}\|\varphi(s)\|\,ds\quad\text{and}\quad\|Y^{\prime}(t)\|\leq\int_{I_{j}}\|\varphi^{\prime}(s)\|\,ds,$

implying that

$\|Y^{\prime}\|^{2}_{L^{2}(I_{j})}\leq
4\|\varphi\|^{2}_{L^{2}(I_{j})}\quad\text{and}\quad\|Y^{\prime}\|^{2}_{L^{2}(I_{j})}\leq\tau_{j}^{2}\|\varphi^{\prime}\|^{2}_{L^{2}(I_{j})}.$
(4.15)

After summing over $j$ and once again applying Corollary 7.1, we arrive at

$\|Y^{\prime}\|^{2}_{L^{2}(J)}\leq
C\tau^{2\alpha}\|(\mathcal{J}_{T}^{1-\alpha}\varphi)^{\prime}\|^{2}_{L^{2}(J)}.$
(4.16)

Now take the inner product of (4.3) with
$-(\mathcal{J}_{T}^{1-\alpha}\varphi)^{\prime}$ in $L_{2}(\Omega)$, and then
integrate in time to obtain

$\displaystyle\|(\mathcal{J}_{T}^{1-\alpha}\varphi)^{\prime}\|^{2}_{L^{2}(J)}-\mathcal{I}(\langle\mathcal{A}\varphi,(\mathcal{J}_{T}^{1-\alpha}\varphi)^{\prime}\rangle)(T)$
$\displaystyle=-\mathcal{I}(\langle\eta,(\mathcal{J}_{T}^{1-\alpha}\varphi)^{\prime}\rangle)(T)$
$\displaystyle\leq\frac{1}{2}\|\eta\|^{2}_{L^{2}(J)}+\frac{1}{2}\|(\mathcal{J}_{T}^{1-\alpha}\varphi)^{\prime}\|^{2}_{L^{2}(J)}.$

By (4.2),

$-\mathcal{I}(\langle\mathcal{A}\varphi,(\mathcal{J}_{T}^{1-\alpha}\varphi)^{\prime}\rangle)(T)=\mathcal{I}(\langle(\mathcal{I}^{1-\alpha}\mathcal{A}^{1/2}\varphi)^{\prime},\mathcal{A}^{1/2}\varphi\rangle)(T)\geq
0$ (4.17)

and therefore
$\|(\mathcal{J}_{T}^{1-\alpha}\varphi)^{\prime}\|^{2}_{L^{2}(J)}\leq\|\eta\|^{2}_{L^{2}(J)}.$
Combining this with (4.14), (4.16), we conclude that

$\sum_{j=1}^{N}\tau_{j}^{-\alpha}\|Y\|_{I_{j}}^{2}\leq
C\tau^{1+\alpha}\|\eta\|^{2}_{L^{2}(J)}\quad\text{and}\quad\|Y^{\prime}\|^{2}_{L^{2}(J)}\leq
C\tau^{2\alpha}\|\eta\|^{2}_{L^{2}(J)}$

and so, applying Lemma 4.2,

$\mathcal{I}(\langle
Y^{\prime},\mathcal{I}^{1-\alpha}Y^{\prime}\rangle)(T)\leq
C\tau^{1+\alpha}\|\eta\|^{2}_{L^{2}(J)}.$ (4.18)

By taking the inner product of (4.3) with $\mathcal{A}\varphi$ and proceeding
as above, we deduce that

$\|\mathcal{A}\varphi\|^{2}_{L^{2}(J)}\leq\|\eta\|^{2}_{L^{2}(J)}.$ (4.19)

Repeating the arguments leading to the first estimate in (4.15) but with
$Y^{\prime}$ replaced by $\mathcal{A}Y^{\prime}$, we see that
$\|\mathcal{A}Y^{\prime}\|^{2}_{L^{2}(J)}\leq\|\mathcal{A}\varphi\|^{2}_{L^{2}(J)}$
and so

$\sum_{j=1}^{N}\tau_{j}^{\alpha}\biggl{(}\int_{I_{j}}\|\mathcal{A}Y^{\prime}\|\,dt\biggr{)}^{2}\leq\sum_{j=1}^{N}\tau_{j}^{1+\alpha}\|\mathcal{A}Y^{\prime}\|^{2}_{L^{2}(I_{j})}\leq
C\tau^{1+\alpha}\|\mathcal{A}\varphi\|^{2}_{L^{2}(J)}.$

Likewise,
$\|\mathcal{A}Y\|_{I_{j}}^{2}\leq\tau_{j}\|\mathcal{A}\varphi\|^{2}_{L^{2}(I_{j})}$
by the arguments leading to (4.12), so

$\sum_{j=1}^{N}\tau_{j}^{\alpha}\|\mathcal{A}Y\|_{I_{j}}^{2}\leq\sum_{j=1}^{N}\tau_{j}^{1+\alpha}\|\mathcal{A}\varphi\|^{2}_{L^{2}(I_{j})}\leq
C\tau^{1+\alpha}\|\mathcal{A}\varphi\|^{2}_{L^{2}(J)}.$

Hence, by Lemma 4.3 and (4.19),

$\mathcal{I}(\langle\mathcal{A}Y,\mathcal{I}^{\alpha}\mathcal{A}Y^{\prime}\rangle)(T)\leq
C\tau^{1+\alpha}\|\eta\|^{2}_{L^{2}(J)}.$ (4.20)

Together, (4.18) and (4.20) imply the desired estimate (4.6). $\quad\Box$

##  5 A fully discrete scheme and error analysis

In this section, we discretize the time-stepping scheme (2.4) in space using
the continuous piecewise-linear Galerkin FEM and hence define a fully-discrete
method. Thus, we introduce a family of regular (conforming) triangulations
$\mathcal{T}_{h}$ of the domain $\overline{\Omega}$ indexed by
$h=\max_{K\in\mathcal{T}_{h}}(h_{K})$, where $h_{K}$ denotes the diameter of
the element $K$. Let $V_{h}$ denote the space of continuous, piecewise-linear
functions with respect to $\mathcal{T}_{h}$ that vanish on $\partial\Omega$.
Let ${\mathcal{W}}(V_{h})\subset{\mathcal{C}}([0,T];V_{h})$ denote the space
of linear polynomials on $\overline{I}_{n}$ for $1\leq n\leq N$, with
coefficients in $V_{h}$. Motivated by the weak formulation of (2.4), our
fully-discrete solution $U_{h}\in{\mathcal{W}}(V_{h})$ is defined by requiring

$\Big{\langle}\int_{I_{n}}\partial_{t}^{\alpha}U_{h}\,dt,v_{h}\Big{\rangle}+\tau_{n}\langle\kappa\nabla
U_{h}^{n-1/2},\nabla
v_{h}\rangle=\tau_{n}\langle\bar{f}_{n},v_{h}\rangle\quad\text{for all
$v_{h}\in V_{h}$,}$ (5.1)

and for $1\leq n\leq N$, where $U_{h}^{n}:=U_{h}(t_{n})$ and
$U_{h}^{n-1/2}=\tfrac{1}{2}(U_{h}^{n}+U_{h}^{n-1})$. For the discrete initial
data, we choose $U_{h}^{0}=R_{h}u_{0}\in V_{h}$, where
$R_{h}:H^{1}_{0}(\Omega)\mapsto V_{h}$ is the Ritz projection defined by
$\langle\kappa\nabla(R_{h}w-w),\nabla v_{h}\rangle=0$ for all $v_{h}\in
V_{h}$.

In the next theorem, we prove that the numerical solution defined by (5.1) is
second order accurate in both time and space, provided that the time mesh
exponent $\gamma$ is chosen appropriately. In comparison, under heavier
regularity assumptions and stronger graded meshes, convergence of order
$h^{2}+\tau^{(3-\alpha)/2}$ was proved by Mustapha et al. [22]. Furthermore,
the proof therein is more technical and lengthy. Use of the piecewise-linear
polynomial function $\widehat{u}$, see (2.6), and a duality argument allowed
us to improve the convergence rate, simplify the proof and also relax the
regularity assumptions. In addition to the regularity assumption in (2.11),
for $t>0,$ we impose

$\|u^{\prime}(t)\|_{2}+t\|u^{\prime\prime}(t)\|_{2}\leq
C\,t^{\upsilon-1},\quad\text{for some~{}$\upsilon>0.$}$ (5.2)

###### Theorem 5.1

Let $u$ be the solution of (1.1) and let $U_{h}^{n}$ be the approximate
solution defined by (5.1). Assume that the regularity assumptions in (2.11)
and (5.2) are satisfied for $\sigma,\upsilon>\alpha/2$, and choose the mesh
grading exponent
$\gamma>\max\\{2/\sigma,1/\upsilon,(3-\alpha)/(2\sigma-\alpha)\\}$. Then,

$\|u-U_{h}\|_{L^{2}(J)}\leq C(\tau^{2}+h^{2}).$

Proof. Decompose the error as $u-U_{h}=(u-u_{h})+(u_{h}-U_{h})$, where $u_{h}$
is the Galerkin finite element solution of (1.1) defined by

$\langle\partial_{t}^{\alpha}u_{h},v_{h}\rangle+\langle\kappa\nabla
u_{h},\nabla v_{h}\rangle=\langle f,v_{h}\rangle\quad\text{for all $v_{h}\in
V_{h}$,}$ (5.3)

for each fixed $t>0,$ with $u_{h}(0)=U_{h}^{0}=R_{h}u_{0}.$ From this, the
weak formulation of (1.1), and the orthogonality property of the Ritz
projection, we have

$\langle\partial_{t}^{\alpha}(u_{h}-R_{h}u),v_{h}\rangle+\langle\kappa\nabla(u_{h}-R_{h}u),\nabla
v_{h}\rangle=\langle\partial_{t}^{\alpha}(u-R_{h}u),v_{h}\rangle\quad\text{for
$v_{h}\in V_{h}$,}$

Choose $v_{h}=(u_{h}-R_{h}u)^{\prime}$, integrate over $(0,t)$ and apply Lemma
2.1 to the right-hand side with $\epsilon=1/(4\alpha^{2})$. After canceling
the common terms,

$4\alpha^{2}\|\sqrt{\kappa}\nabla(u_{h}-R_{h}u)(t)\|^{2}\leq\mathcal{I}(\langle\partial_{t}^{\alpha}e_{h},e_{h}^{\prime}\rangle)(t),~{}~{}{\rm
with}~{}~{}e_{h}=u-R_{h}u\,.$

The error bound for the Ritz projection and the regularity assumption in (5.2)
yield $\|e_{h}^{\prime}(t)\|\leq Ch^{2}\|u^{\prime}(t)\|_{2}\leq
Ch^{2}t^{\upsilon-1}$. Hence, $\|\partial_{t}^{\alpha}e_{h}(t)\|\leq
Ch^{2}t^{\upsilon-\alpha}$ and consequently,
$\mathcal{I}(\langle\partial_{t}^{\alpha}e_{h},e_{h}^{\prime}\rangle)(t)\leq
Ch^{4}$ for $\upsilon>\alpha/2$. Inserting this estimate in the above
equation, we obtain $\|\nabla(u_{h}-R_{h}u)(t)\|\leq Ch^{2}$ for
$\upsilon>\alpha/2$, and thus, by applying the Poincaré and triangle
inequalities, we get $\|u(t)-u_{h}(t)\|\leq Ch^{2}.$

The remaining target now is to estimate $U_{h}-u_{h}.$ By analogy with our
earlier splitting (2.5), we let

$\eta_{h}=u_{h}-U_{h},\quad\psi_{h}=u-R_{h}\widehat{u}\quad{\rm
and}\quad\theta_{h}=U_{h}-R_{h}\widehat{u}.$

From (5.1) and (5.3), and with $\chi_{h}=u_{h}-R_{h}\widehat{u},$ we have

$\int_{I_{n}}[\langle\partial_{t}^{\alpha}\theta_{h}(t),v_{h}\rangle+\langle\kappa\nabla\theta_{h}(t),\nabla
v_{h}\rangle]\,dt=\int_{I_{n}}[\langle\partial_{t}^{\alpha}\chi_{h}(t),v_{h}\rangle+\langle\kappa\nabla\chi_{h}(t),\nabla
v_{h}\rangle]\,dt,$

for $v_{h}\in V_{h}$. Using the orthogonality property of the Ritz projection
and the definition of $\widehat{u}$, in addition to (5.3),

$\int_{I_{n}}\langle\kappa\nabla\chi_{h}(t),\nabla
v_{h}\rangle\,dt=\int_{I_{n}}\langle\kappa\nabla(u_{h}-u)(t),\nabla
v_{h}\rangle\,dt=\int_{I_{n}}\langle\partial_{t}^{\alpha}(u-u_{h})(t),v_{h}\rangle\,dt,$

and hence

$\int_{I_{n}}[\langle\partial_{t}^{\alpha}\theta_{h}(t),v_{h}\rangle+\langle\kappa\nabla\theta_{h}(t),\nabla
v_{h}\rangle]\,dt=\int_{I_{n}}\langle\partial_{t}^{\alpha}\psi_{h}(t),v_{h}\rangle\,dt,\quad\text{for
$v_{h}\in V_{h}$.}$

By repeating the steps from (2.8) to (2.10), we deduce that

$\alpha^{2}\mathcal{I}(\langle\partial_{t}^{\alpha}\theta_{h},\theta_{h}^{\prime}\rangle)(T)\leq\mathcal{I}(\langle\partial_{t}^{\alpha}\psi_{h},\psi_{h}^{\prime}\rangle)(T).$
(5.4)

Applying Lemma 2.2 with
$v=\widehat{e}_{h}^{\prime}:=(\widehat{u}-R_{h}\widehat{u})^{\prime}$ and
$w=\psi^{\prime}$ so that $\psi_{h}^{\prime}=v+w$,

$\mathcal{I}(\langle\partial_{t}^{\alpha}\psi_{h},\psi_{h}^{\prime}\rangle)(T)\leq(1+\alpha^{-1})\Big{(}\mathcal{I}(\langle\partial_{t}^{\alpha}\widehat{e}_{h},\widehat{e}_{h}^{\prime}\rangle)(T)+\mathcal{I}(\langle\partial_{t}^{\alpha}\psi,\psi^{\prime}\rangle)(T)\Big{)}.$

For $t\in I_{j}$, since
$\widehat{u}^{\prime}(t)=\tau_{j}^{-1}(\widehat{u}^{j}-\widehat{u}^{j-1}),$
the Ritz projection error bound gives

$\|\widehat{e}_{h}^{\prime}(t)\|\leq Ch^{2}\|\widehat{u}^{\prime}(t)\|_{2}\leq
Ch^{2}\tau_{j}^{-1}\|\psi^{j}-\psi^{j-1}\|_{2}+Ch^{2}\tau_{j}^{-1}\int_{I_{j}}\|u^{\prime}(s)\|_{2}\,ds\,.$

From Lemma 3.1, (5.2) and the time mesh property (2.2), we have

$\|\psi^{j}-\psi^{j-1}\|_{2}\leq\|\psi^{j}\|_{2}+\|\psi^{j-1}\|_{2}\leq
2\sum_{i=1}^{j}\int_{I_{j}}(t-t_{j-1})\|u^{\prime\prime}(t)\|_{2}\,dt\\\ \leq
C\tau_{1}^{\upsilon}+C\tau\int_{t_{1}}^{t_{j}}t^{\upsilon-1/\gamma-1}\,dt\leq
C\tau_{1}^{\upsilon}+C\tau t_{j}^{\upsilon-1/\gamma}\leq C\tau
t_{j}^{\upsilon-1/\gamma}\leq C\tau_{j}t_{j}^{\upsilon-1},$

for $\gamma>1/\upsilon$ with $j\geq 1.$ Combining the above two estimates, we
conclude that $\|\widehat{e}_{h}^{\prime}(t)\|\leq
Ch^{2}t_{j}^{\upsilon-1}\leq Ch^{2}t^{\upsilon-1}$ for $t\in I_{j}$. This
leads to $\|\widehat{e}_{h}^{\prime}(t)\|\leq Ch^{2}\omega_{\upsilon}(t)$ and
$\|\partial_{t}^{\alpha}\widehat{e}_{h}(t)\|\leq
Ch^{2}\omega_{1-\alpha+\upsilon}(t).$ Thence,

$\mathcal{I}(\langle\partial_{t}^{\alpha}\widehat{e}_{h},\widehat{e}_{h}^{\prime}\rangle)(T)\leq
Ch^{4}\int_{0}^{T}\omega_{1-\alpha+\upsilon}\,\omega_{\upsilon}\,dt\leq
Ch^{4}\int_{0}^{T}t^{2\upsilon-\alpha-1}\,dt\leq Ch^{4},$

for $\upsilon>\alpha/2$, so by Lemma 3.4,

$\mathcal{I}(\langle\partial_{t}^{\alpha}\psi_{h},\psi_{h}^{\prime}\rangle)(T)\leq
C\Big{(}h^{4}+\mathcal{I}(\langle\partial_{t}^{\alpha}\psi,\psi^{\prime}\rangle)(T)\Big{)}\leq
C(h^{4}+\tau^{3-\alpha}).$ (5.5)

Adapting (4.3), suppose that $\varphi_{h}(t)\in V_{h}$ satisfies the discrete
final-value problem

$-(\mathcal{J}_{T}^{1-\alpha}\varphi_{h})^{\prime}(t)+\mathcal{A}_{h}\varphi_{h}(t)=\eta_{h}(t)\quad\text{for
$0<t<T$,}\quad\text{with $\varphi_{h}(T)=0$,}$ (5.6)

where the discrete elliptic operator ${\mathcal{A}}_{h}:V_{h}\to V_{h}$ is
defined by $\langle{\mathcal{A}}_{h}v_{h},q_{h}\rangle=\langle\kappa\nabla
v_{h},\nabla q_{h}\rangle$ for all $v_{h},q_{h}\in V_{h}$. We now repeat the
step in the error analysis of Section 4, with $\eta_{h},$ $\theta_{h}$,
$\psi_{h}$ and $\mathcal{A}_{h}$ playing the roles of $\eta,$ $\theta$ $\psi$
and $\mathcal{A}$, respectively, and using (5.1) and (5.4) instead of (2.4)
and (2.10), and (5.5) instead of Lemma 3.4. We notice for
$\gamma>\max\\{2/\sigma,1/\upsilon,(3-\alpha)/(2\sigma-\alpha)\\}$ and for
$\sigma,\upsilon>\alpha/2,$ that

$\|u_{h}-U_{h}\|_{L^{2}(J)}^{2}\leq C\tau^{1+\alpha}(h^{4}+\tau^{3-\alpha}).$

The proof of this theorem is completed. $\quad\Box$

##  6 Numerical results

In this section, we illustrate numerically the theoretical finding in Theorem
2.1. An $O(h^{2})$ convergence of the finite element solution was confirmed
for various choices of the given data [8, 10, 24]. In time, some numerical
convergence results (piecewise linear discontinuous Petrov-Galerkin method)
were also delivered [22]. However, we illustrate the errors and convergence
rates in the stronger $L^{\infty}((0,T),L^{2}(\Omega))$-norm on more realistic
examples. We choose $\kappa=1$, $\Omega=(0,1)$ and a uniform spatial grid
$\mathcal{T}_{h}$. In both examples, we choose $h$ so that the error from the
time discretization dominates.

### Example 1.

We choose $u_{0}(x)=x(1-x)$ and $f\equiv 0$. Thus, by separating variables,
the continuous solution has a series representation in terms of the Mittag-
Leffler function $E_{\alpha}$,

$u(x,t)=8\sum_{m=0}^{\infty}\lambda_{m}^{-3}E_{\alpha}(-\lambda_{m}^{2}t^{\alpha})\sin(\lambda_{m}x),\quad\text{where
$\lambda_{m}=(2m+1)\pi$.}$

Since $u_{0}\in\dot{H}^{2.5^{-}}(\Omega)$, the regularity estimate (2.11) is
satisfied for $\sigma=\alpha$. Thus, we expect from Theorem 2.1 that
$e_{\tau}:=\|u-U_{h}\|_{L^{2}(J)}\leq C\tau^{2}$ provided that the mesh
exponent $\gamma>\max\\{2/\alpha,(3-\alpha)/\alpha\\}=(3-\alpha)/\alpha$. The
numerical results in Table 1 indicate order $\tau^{2}$ convergence in the
stronger $L^{\infty}((0,T),L^{2}(\Omega))$-norm for $\gamma>2/\alpha$. Rates
of order $\tau^{\sigma\gamma}$ for $1\leq\gamma\leq 2/\sigma=2/\alpha$ are
observed. Thence, our imposed assumption on $\gamma$ is not sharp.

To measure the $L^{\infty}((0,T),L^{2}(\Omega))$ error $E_{\tau}:=\max_{0\leq
t\leq T}\|u-U_{h}\|$, we approximated $E_{\tau}$ by $\max_{1\leq j\leq
N}\max_{1\leq i\leq 3}\|u(t_{i,j})-U_{h}(t_{i,j})\|$ where
$t_{i,j}:=t_{j-1}+i\tau_{j}/3$. In our calculations, the $L^{2}(\Omega)$ norm,
$\|\cdot\|,$ is approximated using the two-point composite Gauss quadrature
rule.

In all tables and figures, we evaluated the series solution by truncating the
infinite series after $60$ terms. The empirical convergence rate CR is
calculated by halving $\tau$, that is,
$\text{CR}=\log_{2}(E_{\tau}/E_{\tau/2})$. Figure 1 plots the nodal errors
$\|U^{n}_{h}-u(t_{n})\|$ against $t_{n}\in[0,1]$ for different values of $N$
in the cases $\gamma=1$ and $\gamma=4$. The practical benefit of the mesh
grading is evident.

| $\gamma=1$ | $\gamma=2$ | $\gamma=3$ | $\gamma=4$
---|---|---|---|---
$N$ |  $E_{\tau}$ CR |  $E_{\tau}$ CR |  $E_{\tau}$ CR |  $E_{\tau}$ CR
8 | 1.011e-01 |  | 5.001e-02 |  | 1.861e-02 |  | 7.478e-03 |
16 | 8.337e-02 | 0.279 | 2.633e-02 | 0.926 | 6.464e-03 | 1.526 | 2.090e-03 | 1.839
32 | 6.588e-02 | 0.340 | 1.309e-02 | 1.008 | 2.259e-03 | 1.517 | 5.497e-04 | 1.927
64 | 5.001e-02 | 0.397 | 6.464e-03 | 1.018 | 7.932e-04 | 1.501 | 1.449e-04 | 1.923
128 | 3.672e-02 | 0.446 | 3.205e-03 | 1.012 | 2.817e-04 | 1.494 | 3.801e-05 | 1.931
Table 1: Errors and empirical convergence rates for Example 1 with
$\alpha=0.5$, using different choices of the time mesh-grading exponent
$\gamma$.

Figure 1: Errors for Example 1 as functions of $t$ for different choices of
$N$ when $\alpha=0.5$, taking $\gamma=1$ in the left figure and $\gamma=4$ in
the right figure.

### Example 2.

We again take $f\equiv 0$ but now choose less regular initial data, namely,
the hat function on the unit interval, $u_{0}(x)=1-2|x-\tfrac{1}{2}|$. So,

$u(x,t)=4\sum_{m=0}^{\infty}(-1)^{m}\lambda_{m}^{-2}E_{\alpha}(-\lambda_{m}^{2}t^{\alpha})\sin(\lambda_{m}x).$

Since $u_{0}\in\dot{H}^{1.5^{-}}(\Omega)$, the regularity property (2.11) is
satisfied for $\sigma=\tfrac{3}{4}\alpha$. As in Example 1, the numerical
results in Table 2 exhibit convergence of order $\tau^{\sigma\gamma}$ for
$1\leq\gamma\leq 2/\sigma$ in the stronger $\|\cdot\|_{J}$-norm. For a
graphical illustration of the impact of the graded mesh on the pointwise
error, we fixed $N=80$ in Figure 2 and plotted the error at the time nodal
points for different choices of $\gamma$.

| $\gamma=1$ | $\gamma=2$ | $\gamma=3$ | $\gamma=4$
---|---|---|---|---
$N$ |  $E_{\tau}$ CR |  $E_{\tau}$ CR |  $E_{\tau}$ CR |  $E_{\tau}$ CR
8 | 9.958e-02 |  | 2.771e-02 |  | 9.365e-03 |  | 6.670e-03 |
16 | 6.501e-02 | 0.615 | 1.345e-02 | 1.042 | 3.144e-03 | 1.575 | 1.710e-03 | 1.963
32 | 4.158e-02 | 0.644 | 6.509e-03 | 1.048 | 1.056e-03 | 1.574 | 4.300e-04 | 1.992
64 | 2.771e-02 | 0.586 | 3.144e-03 | 1.050 | 3.555e-04 | 1.571 | 1.095e-04 | 1.973
128 | 1.924e-02 | 0.526 | 1.519e-03 | 1.049 | 1.193e-04 | 1.576 |  |
Table 2: Errors and empirical convergence rates for Example 2 using different
values of $N$ and different choices of the time mesh exponent $\gamma,$ with
$\alpha=0.7.$ Figure 2: Error at $t_{n}$ for $1\leq n\leq N$ in Example 2, for
a fixed $N=80$ and different choices of the mesh exponent $\gamma$ with
$\alpha=0.7$.

##  7 Appendix: An $\alpha$-robust interpolation estimate

The purpose of this appendix is to prove Corollary 7.1, which was used in the
proof of Theorem 4.2. Suppose that $Y$ is a complex Hilbert space with inner
product $(\cdot,\cdot)_{Y}$ and norm $\|\cdot\|_{Y}$. Put

$X_{0}=L_{2}(0,T;Y)\quad\text{and}\quad X_{1}=\\{\,u\in
X_{0}:\text{$u^{\prime}\in X_{0}$ and $u(0)=0$}\,\\},$ (7.1)

equipped with the inner products

$(u,v)_{X_{0}}=\int_{0}^{T}\bigl{(}u(t),v(t)\bigr{)}_{Y}\,dt\quad\text{and}\quad(u,v)_{X_{1}}=\int_{0}^{T}\bigl{(}u^{\prime}(t),v^{\prime}(t)\bigr{)}_{Y}\,dt.$

We define an unbounded operator $D_{0}$ on $X_{0}$ with domain $X_{1}$, by
setting $D_{0}u=u^{\prime}$ on the interval $(0,T)$. Let $\Re z$ denote the
real part of a complex number $z$. Since
$\Re\bigl{(}u(t),D_{0}u(t)\bigr{)}_{Y}=\tfrac{1}{2}(d/dt)\|u(t)\|_{Y}^{2}$, we
have $\Re(u,D_{0}u)_{X_{0}}=\tfrac{1}{2}\|u(T)\|_{Y}^{2}\geq 0$ for all $u\in
X_{1}$. In addition, one easily verifies that

$(\lambda+D_{0})^{-1}u(t)=\int_{0}^{t}e^{-\lambda(t-s)}u(s)\,ds\quad\text{for
$0\leq t\leq T$ and $\lambda\in\mathbb{C}$,}$ (7.2)

showing that $D_{0}$ is m-accretive [17, p. 9], and therefore a positive
operator [17, pp. 1, 11] with
$\sup_{\lambda>0}\|\lambda(\lambda+D_{0})^{-1}\|_{X_{0}\to X_{0}}=1$, where
$\|\cdot\|_{X_{0}\to X_{0}}$ denotes the operator norm induced by the norm of
$X_{0}$.

Suppose that $0<\Re\beta<1$. The operator $D_{0}^{-\beta}$ is defined via the
integral [2, Equation (1.3)]

$D_{0}^{-\beta}u=\frac{\sin\pi\beta}{\pi}\int_{0}^{\infty}\lambda^{-\beta}(\lambda+D_{0})^{-1}u\,d\lambda\quad\text{for
$u\in X_{0}$.}$

A short calculation using (7.2) then shows [2, Theorem 2.2] that
$D_{0}^{-\beta}u=\mathcal{I}^{\beta}u$ for $0<\Re\beta<1$. One may then define
[2, Equation 1.4], [17, Equation (3.1)]

$D_{0}^{\beta}u=D_{0}^{\beta-1}D_{0}u=\frac{\sin\pi\beta}{\pi}\int_{0}^{\infty}\lambda^{\beta-1}(\lambda+D_{0})^{-1}D_{0}u\,d\lambda\quad\text{for
$u\in X_{1}$,}$

and
$D_{0}^{\beta}u=\mathcal{I}^{1-\beta}u^{\prime}=(\mathcal{I}^{1-\beta}u)^{\prime}$,
that is, $D_{0}^{\beta}$ coincides with both the Caputo and Riemann–Liouville
fractional derivatives of order $\beta$.

For $0\leq\alpha\leq 1$, let $X_{\alpha}$ denote the complex interpolation
space [4, Chap. 4] arising from $X_{0}$ and $X_{1}$. The next result is known
[17, Theorem 11.6.1], but we sketch the proof in order to verify that the
constant does not depend on $\alpha$.

###### Theorem 7.1

Let $0<\alpha<1$. For the spaces (7.1) and any complex Hilbert space $Z$, if
the linear operator $B$ satisfies $\|Bu\|_{Z}\leq M_{j}\|u\|_{X_{j}}$ for
$u\in X_{j}$ and $j\in\\{0,1\\}$, then

$\|Bu\|_{Z}\leq
e^{1+(\pi/4)^{2}}M_{0}^{1-\alpha}M_{1}^{\alpha}\|D_{0}^{\alpha}u\|_{X_{0}}\quad\text{for
$D_{0}^{\alpha}u\in X_{0}$.}$

Proof. If we set $Z_{0}=Z_{1}=Z$, then $Z_{\alpha}=Z$ with equal norms [4,
Theorem 4.2.1], so $\|Bu\|_{Z}\leq
M_{0}^{1-\alpha}M_{1}^{\alpha}\|u\|_{X_{\alpha}},$ and our task is to estimate
the interpolation norm $\|u\|_{X_{\alpha}}$ in terms of the fractional
derivative $D_{0}^{\alpha}u$. Define the closed strip
$S=\\{\,z\in\mathbb{C}:0\leq\Re z\leq 1\,\\}$ in the complex plane, and let
$\mathcal{F}$ denote the space of functions $f:S\to X_{0}+X_{1}$ that are
bounded and continuous on $S$, analytic in the interior of $S$, and such that
$f(iy)$ and $f(1+iy)$ tend to zero as $|y|\to\infty$. It can be shown [4,
Lemma 4.1.1] that $\mathcal{F}$ is a Banach space with respect to the norm

$\|f\|_{\mathcal{F}}=\max\Bigl{(}\sup_{y\in\mathbb{R}}\|f(iy)\|_{X_{0}},\sup_{y\in\mathbb{R}}\|f(1+iy)\|_{X_{1}}\Bigr{)}.$

The space $X_{\alpha}$ then consists of those $u\in X_{0}+X_{1}$ such that
$u=f(\alpha)$ for some $f\in\mathcal{F}$, with the interpolation norm defined
by

$\|u\|_{X_{\alpha}}=\inf\\{\,\|f\|_{\mathcal{F}}:\text{$u=f(\alpha)$ and
$f\in\mathcal{F}$}\,\\}.$

If we define $f(z)=e^{(z-\alpha)^{2}}D_{0}^{\alpha-z}u$, then $f(\alpha)=u$
with
$\|f(iy)\|_{X_{0}}=e^{\alpha^{2}-y^{2}}\|D_{0}^{-iy}D_{0}^{\alpha}u\|_{X_{0}}$
and

$\|f(1+iy)\|_{X_{1}}=e^{(1-\alpha)^{2}-y^{2}}\|D_{0}(D_{0}^{\alpha-1-iy}u)\|_{X_{0}}=e^{(1-\alpha)^{2}-y^{2}}\|D_{0}^{-iy}D_{0}^{\alpha}u\|_{X_{0}}.$

It is known [26, Example 2] that $\|D_{0}^{iy}u\|_{X_{0}}\leq
e^{\pi|y|/2}\|u\|_{X_{0}}$ for $u\in X_{0}$ because $D_{0}$ is m-accretive.
Thus

$\|u\|_{X_{\alpha}}\leq\|f\|_{\mathcal{F}}\leq\max\Bigl{(}\sup_{y\in\mathbb{R}}e^{\alpha^{2}-y^{2}+\pi|y|/2},\sup_{y\in\mathbb{R}}e^{(1-\alpha)^{2}-y^{2}+\pi|y|/2}\Bigr{)}\|D_{0}^{\alpha}u\|_{X_{0}},$

and since
$-y^{2}+\tfrac{1}{2}\pi|y|=(\tfrac{1}{4}\pi)^{2}-(|y|-\tfrac{1}{4}\pi)^{2}\leq(\tfrac{1}{4}\pi)^{2},$
the proof is completed. $\quad\Box$

The interpolation estimate used in the proof of Theorem 4.2 now follows by
repeating the preceding arguments with $D_{0}$ and $X_{1}$ replaced by their
time-reversed equivalents

$D_{T}u(t)=-u^{\prime}(t)\quad\text{and}\quad X_{T,1}=\\{\,u\in
X_{0}:\text{$u^{\prime}\in X_{0}$ and $u(T)=0$}\,\\}.$

In fact, $\Re(u,D_{T}u)_{X_{0}}=\tfrac{1}{2}\|u(0)\|_{Y}^{2}\geq 0$ for all
$u\in X_{T,1}$, and

$(\lambda+D_{T})^{-1}u(t)=\int_{t}^{T}e^{-\lambda(s-t)}u(s)\,ds\quad\text{for
$0\leq t\leq T$ and $\lambda\in\mathbb{C}$,}$

so $D_{T}$ is m-accretive. We find that if $0<\Re\beta<1$ then
$D_{T}^{-\beta}u=\mathcal{J}_{T}^{\beta}u$, with

$D_{T}^{\beta}u=D_{T}D_{T}^{\beta-1}u=(\mathcal{J}_{T}^{1-\beta}u)^{\prime}\quad\text{for
$u\in X_{T,1}$.}$

###### Corollary 7.1

In the statement of Theorem 7.1, we may replace $D_{0}$ and $X_{1}$ with
$D_{T}$ and $X_{T,1}$, respectively.

## References

  * [1] A. A. Alikhanov, A new difference scheme for the time fractional diffusion equation, J. Comput. Phys., 280, (2015), 424–438.
  * [2] A. Ashyralyev, A note on fractional derivatives and fractional powers of operators, J. Math. Anal. Appl. 357 (2009), 232–236.
  * [3] L. Banjai and C. G. Makridakis, A posteriori error analysis for approximations of time-fractional subdiffusion problems, Math. Comp., 91 (2022), 1711–1737.
  * [4] J. Bergh and J. Löfström, Interpolation Spaces: An Introduction, Springer-Verlag, 1976.
  * [5] Z. Cen, J. Huang, A. Le, A. Xu, An efficient numerical method for a time-fractional diffusion equation, preprint, 2018, https://arxiv.org/abs/1810.07935.
  * [6] M. Gunzburger and J. Wang, A second-order Crank-Nicolson method for time-fractional PDEs, Int. J. Numer. Anal. Model., 16, (2019), 225–239.
  * [7] D. Hou, M. T. Hasan, and C. Xu, Müntz spectral methods for the time-fractional diffusion equation, Comput. Methods Appl. Math., 18, (2018), 43–62.
  * [8] B. Jin and Z. Zhou, Correction of high-order BDF convolution quadrature for fractional evolution equations, SIAM J. Sci. Comput., 39, (2017), A3129–A3152.
  * [9] B. Jin, B. Li, and Z. Zhou, An analysis of the Crank-Nicolson method for subdiffusion, IMA J. Numer. Anal., 38 (2018), 518–541.
  * [10] B. Jin, B. Li and Z. Zhou, Subdiffusion with time-dependent coefficients: improved regularity and second-order time stepping, Numer. Math., 145, (2020), 883–913.
  * [11] B. Jin and Z. Zhou, Numerical Treatment and Analysis of Time-Fractional Evolution Equations, Springer, 2023.
  * [12] S. Karaa, Semidiscrete finite element analysis of time fractional parabolic problems: a unified approach, SIAM J. Numer. Anal., 56, (2018), 1673–1692.
  * [13] N. Kopteva, Error analysis for time-fractional semilinear parabolic equations using upper and lower solutions, SIAM J. Numer. Anal., 58, (2020), 2212–2234.
  * [14] N. Kopteva, Error analysis of an L2-type method on graded meshes for a fractional-order parabolic problem, Math. Comp., 90, (2021), 19–40.
  * [15] K. N. Le, W. McLean, and K. Mustapha, A semidiscrete finite element approximation of a time-fractional Fokker–Planck equation with non-smooth initial data, SIAM J. Sci. Comput., 40, (2018), A3831–3852.
  * [16] H.-L. Liao, D. Li, and J. Zhang, Sharp error estimate of the nonuniform L1 formula for linear reaction-subdiffusion equations, SIAM J. Numer. Anal., 56, (2018), 1112–1133.
  * [17] C. Martínez Carracedo and M. Sanz Alix, The Theory of Fractional Powers of Operators, North-Holland, 2001.
  * [18] M. Al-Maskari and S. Karaa, Numerical approximation of semilinear subdiffusion equations with nonsmooth initial data, SIAM J. Numer. Anal., 57 (2019).
  * [19] W. McLean and K. Mustapha, A second-order accurate numerical method for a fractional wave equation, Numer. Math., 105, (2007), 481–510.
  * [20] W. McLean, Regularity of solutions to a time-fractional diffusion equation, ANZIAM J., 52, (2010), 123–138.
  * [21] 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.
  * [22] K. Mustapha, B. Abdallah, and K. M. Furati, A discontinuous Petrov-Galerkin method for time-fractional diffusion equations, SIAM J. Numer. Anal., 52, (2014), 2512–2529.
  * [23] K. Mustapha and D. Schötzau, Well-posedness of $hp-$version discontinuous Galerkin methods for fractional diffusion wave equations, IMA J. Numer. Anal., 34, (2014), 1226–1246.
  * [24] K. Mustapha, FEM for time-fractional diffusion equations, novel optimal error analyses, Math. Comp., 87, (2018), 2259–2272.
  * [25] K. Mustapha, An L1 Approximation for a fractional reaction-diffusion equation, a second-order error analysis over time-graded meshes, SIAM J. Numer. Anal., 58, (2000), 1319–1338.
  * [26] J. Prüss and H. Sohr, On operators with bounded imaginary powers in Banach spaces, Math. Z. 203 (1990), 429–452.
  * [27] 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, (2011), 426–447.
  * [28] J. Shen and C.-T. Sheng, An efficient space-time method for time fractional diffusion equation, J. Sci. Comput., 81, (2019), 1088–1110.
  * [29] M. Styness, A survey of the L1 scheme in the discretisation of time-fractional problems, Numer. Math. Theor. Meth. Appl., 15, (2022), 1173–1192.
  * [30] 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, (2017), 1057–1079.
  * [31] J. Wang, J. Wang, and L. Yin, A single-step correction scheme of Crank-Nicolson convolution quadrature for the subdiffusion equation, J. Sci. Comput., 87, (2021).
  * [32] Y. Wang, Y. Yan, and Y. Yang, Two high-order time discretization schemes for subdiffusion problems with nonsmooth data, Fract. Calc. Appl. Anal., 23, (2020), 1349–1380.
  * [33] S. Wu and Z. Zhou, A parallel-in-time algorithm for high-order BDF methods for diffusion and subdiffusion equations, SIAM J. Sci. Comput., 43 (2021), A3627–A3656.
  * [34] Y. Yan, B. A. Egwu, Z. Liang, and Y. Yan, Error estimates of a continuous Galerkin time stepping method for subdiffusion problem, J. Sci. Comput., 88, (2021).
  * [35] H. Zhang, F. Zeng, X. Jiang, and Z. Zhang, Fast time-stepping discontinuous Galerkin method for the subdiffusion equation, https://arxiv.org/abs/2309.02988, 2023.
  * [36] F. Zeng, C. Li, F. Liu, and I. Turner, Numerical algorithms for time-fractional subdiffusion equation with second-order accuracy, SIAM J. Sci. Comput., 37, (2015), A55–A78.
  * [37] H. Zhu and C. Xu, A fast high order method for the time-fractional diffusion equation, SIAM J. Numer. Anal., 57, (2019), 2829–2849.


remarkRemark hypothesisHypothesis claimClaim Implicitly Restarted LanczosP.S.
Negi and C. Arratia [][nocite]ex_supplement

# Structure preserving restarts of the non-symmetric Lanczos algorithm via the
implicitly shifted LR algorithm††thanks: Submitted to the editors DATE. This
work was funded by the Nordita the Swedish Research Council Grant No.
2018-04290. Nordita is partially supported by Nordforsk.

P. S. Negi Nordita, Stockholm University, KTH Royal Institute of Technology,
Stockholm, Sweden<EMAIL_ADDRESS>C. Arratia Nordita, Stockholm
University, KTH Royal Institute of Technology, Stockholm, Sweden ().
<EMAIL_ADDRESS>
###### Abstract

The implicitly shifted QR iteration is used as a restart procedure for the
Arnoldi method for the calculation of a few dominant eigenvalues of a large
matrix. We show that the underlying idea of implicit polynomial filtering can
be utilized in much the same manner via the implicitly shifted LR iteration to
create a restart procedure for the non-symmetric Lanczos algorithm for
eigenvalue computations, which preserves the tri-diagonal structure of the
reduced matrix.

###### keywords:

LR algorithm, non-symmetric Lanczos, implicit restart

68Q25, 68R10, 68U05

##  1 Introduction

The Arnoldi iteration [1] is a popular Krylov space method for calculating a
few eigenvalues of a large matrix. The method relies on the generation of a
sequence of Krylov vectors which determine the subspace within which
approximations of the eigenvalue-eigenvector pairs are obtained. Depending on
the accuracy and number of eigenpair approximations needed, the Krylov space
size can become exceedingly large so that the quality of the results may be
limited by the available memory. Sorensen [18] introduced an elegant procedure
for restarting the Arnoldi factorization based on polynomial filters, which
are applied through the implicitly shifted QR iterations on the reduced
Hessenberg matrix obtained through the Arnoldi method. In particular, the use
of exact shifts was shown to be successful in the convergence process of the
eigenspace[18] of the specified eigenvalues. The method has subsequently found
widespread application through the ARPACK library [11]. An alternative, more
modern method for eigenvalue computations is the Krylov-Schur method
introduced by Stewart [19], which can be found implemented in the SLEPc
library [9]. Krylov-Schur restarts are known to be less sensitive than the
restarting based on the implicitly shifted QR method. However, they do not
preserve the Hessenberg structure of the reduced matrix which may sometimes be
required. The use of QR iterations ensures that the reduced matrix preserves
its Hessenberg structure through the transforms that make up the restart
process. If the underlying matrix is symmetric, the Arnoldi iteration reduces
to the Lanczos algorithm, and the Hessenberg matrix reduces to a symmetric
tridiagonal matrix. The QR iteration preserves the symmetric tridiagonal
structure as well and, as pointed out by Sorensen [18], the implicit restart
process applies equally well for the Lanczos method for symmetric matrices.

The Lanczos algorithm introduced in [10] is in fact the predecessor of the
Arnoldi method and can be used for non-symmetric matrices as well. In such
cases it is referred to as the non-symmetric or the bi-Lanczos. It has been
used in investigationss for plasma instabilities [13], linear systems solvers
[8], model reduction [12], and more recently for non-linear eigenvalue
problems [5]. For the non-symmetric Lanczos one builds two Krylov subspaces
referred to as the right and the left Krylov subspaces. The idea behind the
two methods is similar, which is to obtain a projection of the original large
matrix on to an appropriate reduced subspace such that, the eigenvalues may be
approximated via the eigenvalues of the reduced operator. The difference being
that in the Arnoldi method one obtains an orthogonal projection on to a
subspace while, in the Lanczos method one obtains an oblique projection. One
would therefore like to extend the idea of implicit restarts to the non-
symmetric Lanczos algorithm as well. However, the reduced matrix that one
obtains in such a case is a non-symmetric tridiagonal matrix, with the
tridiagonal structure being the result of the recurrence relations of the
Lanczos algorithm [17]. Since the QR iterations do not preserve the banded
structure of non-symmetric matrices, a straightforward application of the
restart procedure put forward by Sorensen will lead to a loss of this
tridiagonal structure (the Hessenberg structure will still be preserved). This
loss of structure can be circumvented if one looks to the predecessor of the
QR algorithm namely, the LR algorithm proposed by Rutishauser [14, 16, 15],
which has the attractive property of preserving the band structure of a
matrix. This property was already pointed out by Rutishauser in [14] where the
banded matrices were referred to as striped matrices. As we will show in the
next section, the shifted LR transform is the appropriate analogue of the
restart procedure in the case of non-symmetric Lanczos iteration. The process
would necessarily require refining both the right as well as the left Krylov
spaces simultaneously.

The rest of the paper is organized as follows. In section 2 we start with the
introduction of the non-symmetric Lanczos iteration and then develop the
restart procedure. In section 3 we apply the restart process to the Grcar
matrix, and make some concluding remarks in section 4.

##  2 Non-symmetric Lanczos

Lanczos first introduced his algorithm in [10] as a method for
tridiagonalizing a matrix, but also realized that the method could be used
iteratively to find eigenvalues. For an arbitry matrix $A$, the method
generates a pair of Krylov subspaces $\\{v_{1},\ldots,v_{m}\\}$ and
$\\{w_{1},\ldots,w_{m}\\}$, through repeated action of $A$ and $A^{H}$
respectively. We refer to these as the right and left Krylov spaces
respectively and they satisfy the biorthogonality relation
$w_{i}^{H}v_{j}=\delta_{ij}$. The two subspaces are generated through a three
term recurrence relation which, for a Krylov space of size $m$, can be written
in matrix form as

(1a) $\displaystyle AV_{m}-V_{m}T_{m}$ $\displaystyle=$
$\displaystyle\widetilde{v}_{m+1}e_{m}^{T},$ (1b) $\displaystyle
A^{H}W_{m}-W_{m}T^{H}_{m}$ $\displaystyle=$
$\displaystyle\widetilde{w}_{m+1}e_{m}^{T},$ (1c) $\displaystyle
W_{m}^{H}V_{m}$ $\displaystyle=$ $\displaystyle I_{m},$

where, $I_{m}$ represents the Identity matrix of size $m$, $T_{m}$ is a tri-
diagonal matrix of size $m$ with $T_{m}^{H}$ it’s Hermitian conjugate, and
$e_{m}$ is the standard unit vector. $\widetilde{v}_{m+1}$ and
$\widetilde{w}_{m+1}$ represent the residual vectors at the $m^{th}$ step. If
either $\widetilde{v}_{m+1}$ or $\widetilde{w}_{m+1}$ vanishes it represents
the convergence of the right or the left Krylov subspaces to an invariant
subspace of dimension $m$. A more serious breakdown occurs if
$\widetilde{w}_{m+1}^{H}\widetilde{v}_{m+1}=0$ with both
$\widetilde{v}_{m+1}\neq 0$ and $\widetilde{w}_{m+1}\neq 0$, in which case a
look-ahead strategy may be employed. We do not address the issues with
breakdown here since it is not specifically related to the restart procedure.
We refer the reader to [13] for the look-ahead Lanczos and to [8] for a
comprehensive overview on Lanczos type solvers and the related issues of
breakdown.

As Sorensen points out for the Arnoldi method [18], if one is interested in an
invariant subspace of dimension $m$, the starting vector of the Krylov
subspace must not contain components of the generator of a cyclic subspace of
dimension greater than $m$. This applies equally for the right and left Krylov
subspaces generated through the Lanczos recurrence relations. Hence a non-
vanishing $\widetilde{v}_{m+1}$ (respectively $\widetilde{w}_{m+1}$) implies
that $v_{1}$ (respectively $w_{1}$) contains components of an invariant
subspace of dimension greater than $m$. The idea behind restarts then is to
discard the components of the starting vector $v_{1}$ (and $w_{1}$) along the
unwanted dimensions, such that each restart process moves the Krylov space(s)
closer to being invariant. For the Arnoldi method Sorensen [18] proposed to
achieve this via polynomial filtering, i.e. replacing

(2a) $\displaystyle v_{1}\leftarrow\psi(A)v_{1},$ (2b)
$\displaystyle\psi(\lambda)=(1/\tau)\Pi_{j=1}^{p}(\lambda-\mu_{j}).$

Obviously $\psi(\lambda)$ is the filtering polynomial, $\tau$ is a
normalization constant and each $\mu_{j}$ specifies a node of the polynomial.
The polynomial acts on $v_{1}$ to filter out the part of the spectrum of $A$
that is close to each $\mu_{j}$. If a particular $\mu_{j}$ corresponds to an
exact eigenvalue of $A$, then components of the corresponding eigenvector are
completely filtered out from $v_{1}$ (at least in exact arithmetic). The node
$\mu_{j}$ is referred to as a shift since the application of the polynomial
filtering relies on the shifted QR algorithm, where $\mu_{j}$ is used as the
shift. As shown below for the case of a single shift, an analogous procedure
can be followed using a shifted LR algorithm which achieves the same effect of
applying a polynomial filter to the starting vector $v_{1}$. Starting with the
Lanczos relation for the right subspace (1a), and adding and subtracting $\mu
V_{m}$ we obtain

(3a) $\displaystyle(A-\mu I)V_{m}-V_{m}(T_{m}-\mu I)$ $\displaystyle=$
$\displaystyle\widetilde{v}_{m+1}e_{m}^{T}$ (3b) $\displaystyle(A-\mu
I)V_{m}-V_{m}(L_{1}R_{1})$ $\displaystyle=$
$\displaystyle\widetilde{v}_{m+1}e_{m}^{T}$ (3c) $\displaystyle(A-\mu
I)V_{m}L_{1}-V_{m}(L_{1}R_{1})L_{1}$ $\displaystyle=$
$\displaystyle\widetilde{v}_{m+1}e_{m}^{T}L_{1}$ (3d) $\displaystyle
A(V_{m}L_{1})-(V_{m}L_{1})(R_{1}L_{1}+\mu I)$ $\displaystyle=$
$\displaystyle\widetilde{v}_{m+1}e_{m}^{T}L_{1}$ (3e) $\displaystyle
AV^{\prime}_{m}-V^{\prime}_{m}T^{\prime}_{m}$ $\displaystyle=$
$\displaystyle\widetilde{v}_{m+1}e_{m}^{T}L_{1}.$

Here we have set $V^{\prime}_{m}=V_{m}L_{1}$ and
$T^{\prime}_{m}=(R_{1}L_{1}+\mu I)$. The matrices $L_{1},R_{1}$ are the lower
and upper triangular matrices obtained from the LU decomposition of
$(T_{m}-\mu I)$. The matrix $L_{1}$ can be required to be unit triangular (all
entries on the main diagonal are ones), in which case the LU decomposition is
unique. Furthermore, $L_{1}$ for a tridiagonal matrix only consists of one
sub-diagonal (in addition to the main diagonal). One can easily recognize that
the new reduced matrix $T^{\prime}_{m}$ is a result of one step of the shifted
LR iteration. Hence $T^{\prime}_{m}$ retains the tridiagonal structure of the
$T_{m}$ [14, 15]. The relationship between starting vectors of the two spaces
$V_{m}$ and $V^{\prime}_{m}$ can be obtained by multiplying (3b) by $e_{1}$,
i.e.

$\displaystyle(A-\mu I)V_{m}e_{1}-(V^{\prime}_{m})R_{1}e_{1}$ $\displaystyle=$
$\displaystyle\widetilde{v}_{m+1}e_{m}^{T}e_{1}$ $\displaystyle\implies(A-\mu
I)v_{1}$ $\displaystyle=$ $\displaystyle v^{\prime}_{1}\rho_{11},$

where $\rho_{11}=e^{T}_{1}R_{1}e_{1}$. This clearly shows the filtering
operation done on the original vector $v_{1}$ to generate the new vector
$v^{\prime}_{1}$.

Since the Lanczos method creates a biorthogonal basis, one must simultaneously
transform the left basis $W_{m}$ to maintain the biorthogonality property. It
is easy to see that the necessary transform to maintain biorthogonality is
$W^{\prime}_{m}=W_{m}L_{1}^{-H}$, since,

$\displaystyle W^{\prime
H}V^{\prime}_{m}=(W_{m}L_{1}^{-H})^{H}(V_{m}L_{1})=L_{1}^{-1}(W^{H}_{m}V_{m})L_{1}=I$

Substituting the relation obtained from the LU decomposition
$T_{m}^{H}=(R_{1}^{H}L_{1}^{H}+\bar{\mu}I)$ in to equation (1b) we can obtain
the modified Lanczos relation for the left Krylov space as

(4a) $\displaystyle A^{H}W_{m}-W_{m}(R^{H}_{1}L^{H}_{1}+\bar{\mu}I)$
$\displaystyle=$ $\displaystyle\widetilde{w}_{m+1}e_{m}^{T}$ (4b)
$\displaystyle
A^{H}(W_{m}L_{1}^{-H})-W_{m}(L_{1}^{-H}L_{1}^{H})(R^{H}_{1}L^{H}_{1}+\bar{\mu}I)L_{1}^{-H}$
$\displaystyle=$ $\displaystyle\widetilde{w}_{m+1}e_{m}^{T}L_{1}^{-H}$ (4c)
$\displaystyle
A^{H}(W_{m}L_{1}^{-H})-(W_{m}L_{1}^{-H})(L^{H}_{1}R^{H}_{1}+\bar{\mu}I)$
$\displaystyle=$ $\displaystyle\widetilde{w}_{m+1}e_{m}^{T}L_{1}^{-H}$ (4d)
$\displaystyle A^{H}W^{\prime}_{m}-W^{\prime}_{m}(T^{\prime}_{m})^{H}$
$\displaystyle=$ $\displaystyle\widetilde{w}_{m+1}e_{m}^{T}L_{1}^{-H}.$

Conveniently the structure of the Lanczos iteration for the left Krylov space
is also preserved. Noting that $L_{1}^{-H}$ is upper triangular, one can again
expose the relationship between the generating vectors of the two left Krylov
spaces as

$\displaystyle W^{\prime}e_{1}$ $\displaystyle=$
$\displaystyle(W_{m}L_{1}^{-H})e_{1}$ $\displaystyle\implies w^{\prime}_{1}$
$\displaystyle=$ $\displaystyle w_{1}(e_{1}^{T}L_{1}^{-H}e_{1}).$

Clearly $w^{\prime}_{1}$ is simply a scalar multiple of the old vector $w_{1}$
and no filtering of the generating vector has occurred. In order to ensure
that we filter the left Krylov space as well, we perform one step of the
shifted LR iteration with the conjugated shift $\bar{\mu}$ on the reduced
matrix $(T^{\prime}_{m})^{H}$ obtained in equation (4d).

(5a) $\displaystyle(A-\mu
I)^{H}W^{\prime}_{m}-W^{\prime}_{m}(T^{\prime}_{m}-\mu I)^{H}$
$\displaystyle=$ $\displaystyle\widetilde{w}_{m+1}e_{m}^{T}L_{1}^{-H},$ (5b)
$\displaystyle(A-\mu I)^{H}W^{\prime}_{m}-W^{\prime}_{m}(L_{2}R_{2})$
$\displaystyle=$ $\displaystyle\widetilde{w}_{m+1}e_{m}^{T}L_{1}^{-H},$ (5c)
$\displaystyle(A-\mu
I)^{H}W^{\prime}_{m}L_{2}-W^{\prime}_{m}(L_{2}R_{2})L_{2}$ $\displaystyle=$
$\displaystyle\widetilde{w}_{m+1}e_{m}^{T}L_{1}^{-H}L_{2},$ (5d)
$\displaystyle
A(W^{\prime}_{m}L_{2})-(W^{\prime}_{m}L_{2})(R_{2}L_{2}+\bar{\mu}I)$
$\displaystyle=$ $\displaystyle\widetilde{w}_{m+1}e_{m}^{T}L_{1}^{-H}L_{2},$
(5e) $\displaystyle
AW^{\prime\prime}_{m}-W^{\prime\prime}_{m}(T^{\prime\prime}_{m})^{H}$
$\displaystyle=$ $\displaystyle\widetilde{w}_{m+1}e_{m}^{T}L_{1}^{-H}L_{2}.$

One may again obtain the relation between the starting vectors of the two
spaces as

(6)
$\displaystyle(A^{H}-\bar{\mu}I)w^{\prime}_{1}=w^{\prime\prime}_{1}(e_{1}^{T}R_{2}e_{1}),$

which clearly shows the filtering operation performed on the starting vector
of the left Krylov space. Obviously the appropriate transform for the right
subspace to maintain orthogonality is
$V^{\prime\prime}_{m}=V^{\prime}_{m}L_{2}^{-H}=V_{m}L_{1}L_{2}^{-H}$. Again,
note that the upper triangular $L_{2}^{-H}$ implies that the new
$v^{\prime\prime}_{1}$ is simply the scalar multiple of $v^{\prime}_{1}$ and
no additional filtering occurs for $v_{1}$ in this step. One can write the
corresponding modified Lanczos relation as

(7) $\displaystyle
AV^{\prime\prime}_{m}-V^{\prime\prime}_{m}T^{\prime\prime}_{m}$
$\displaystyle=$ $\displaystyle\widetilde{v}_{m+1}e_{m}^{T}L_{1}L_{2}^{-H}.$

The above process can be repeated for $p$ unwanted shifts. We denote by
$L_{1}^{p}=L_{11}L_{12}\ldots L_{1p}$ as the product of the lower triangular
matrices generated due to $p$ shifted-LR steps for the right Krylov space
$V_{m}$, and by $L^{p}_{2}=L_{21}L_{22}\ldots L_{2p}$ as the product of the
lower triangular matrices due to the $p$ shifted-LR iterations for the left
Krylov space. Then for a Krylov space size of $m=k+p$ we have two modified
Lanczos relations

(8a) $\displaystyle
AV^{\prime\prime}_{k+p}-V^{\prime\prime}_{k+p}T^{\prime\prime}_{k+p}$
$\displaystyle=$
$\displaystyle\widetilde{v}_{k+p+1}e_{k+p}^{T}L_{1}^{p}(L_{2}^{p})^{-H},$ (8b)
$\displaystyle
A^{H}W^{\prime\prime}_{k+p}-W^{\prime\prime}_{k+p}(T^{\prime\prime}_{k+p})^{H}$
$\displaystyle=$
$\displaystyle\widetilde{w}_{k+p+1}e_{k+p}^{T}(L_{1}^{p})^{-H}L_{2}^{p}.$

We may take a closer look at the structure of the residual matrices on the
right hand side of equation (8a). $L_{1}^{p}$ is a product of $p$ matrices
that are lower triangular with just one subdiagonal. $L_{1}^{p}$ then is lower
triangular with $p$ non-zero subdiagonals. $(L_{2}^{p})^{-H}$ is upper
triangular and the product $L_{1}^{p}(L_{2}^{p})^{-H}$ therefore has $p$ non-
zero subdiagonals. Left multiplication by $e_{k+p}^{T}$ therefore has the form

$\displaystyle
e_{k+p}^{T}L_{1}^{p}(L_{2}^{p})^{-H}=(\underbrace{0,0\ldots,\theta_{k+p}}_{k},\underbrace{b_{1}^{T}}_{p})$

where, $\theta_{k+p}=e_{k+p}^{T}(L_{1}^{p}(L_{2}^{p})^{-H})e_{k}$. Therefore
the matrix on the right hand side of (8a) has zeros in the first $k-1$ columns
and the $k^{th}$ column is simply $\theta_{k+p}\widetilde{v}_{k+p+1}$. The
remaining columns are non-zero in general. A very similar structure is
obtained for the residual matrix in the right hand side of equation (8b) with
zeros in the first $k-1$ columns and the $k^{th}$ column being equal to
$\phi_{k+p}\widetilde{w}_{k+p+1}$, with
$\phi_{k+p}=e_{k+p}^{T}((L_{1}^{p})^{-H}L_{2}^{p})e_{k}$

Partitioning the matrices such that

$\displaystyle
V^{\prime\prime}_{k+p}=(V^{\prime\prime}_{k},V^{\prime\prime}_{p}),$
$\displaystyle\hskip
10.0ptT^{\prime\prime}_{k+p}=\begin{pmatrix}T^{\prime\prime}_{k}&\delta_{k+1}e_{k}e_{1}^{T}\\\
\beta_{k+1}e_{1}e_{k}^{T}&T^{\prime\prime}_{p}\end{pmatrix},$ $\displaystyle
W^{\prime\prime}_{k+p}=(W^{\prime\prime}_{k},W^{\prime\prime}_{p}),$
$\displaystyle\hskip
10.0pt(T^{\prime\prime}_{k+p})^{H}=\begin{pmatrix}(T^{\prime\prime}_{k})^{H}&\bar{\beta}_{k+1}e_{k}e_{1}^{T}\\\
\bar{\delta}_{k+1}e_{1}e_{k}^{T}&(T^{\prime\prime}_{p})^{H}\end{pmatrix},$

with the length of the $e_{i}$ vectors understood to be such that the
resulting matrices are consistent. We can write the modified Lanczos relations
of (8a) and (8b) as

$\displaystyle
A(V^{\prime\prime}_{k},V^{\prime\prime}_{p})=(V^{\prime\prime}_{k},V^{\prime\prime}_{p})\begin{pmatrix}T^{\prime\prime}_{k}&\delta_{k+1}e_{k}e_{1}^{T}\\\
\beta_{k+1}e_{1}e_{k}^{T}&T^{\prime\prime}_{p}\end{pmatrix}+\begin{pmatrix}\theta_{k+p}\widetilde{v}_{k+p+1}e_{k}^{T},M_{v}\end{pmatrix},$
$\displaystyle
A^{H}(W^{\prime\prime}_{k},W^{\prime\prime}_{p})=(W^{\prime\prime}_{k},W^{\prime\prime}_{p})\begin{pmatrix}(T^{\prime\prime}_{k})^{H}&\bar{\beta}_{k+1}e_{k}e_{1}^{T}\\\
\bar{\delta}_{k+1}e_{1}e_{k}^{T}&(T^{\prime\prime}_{p})^{H}\end{pmatrix}+\begin{pmatrix}\phi_{k+p}\widetilde{w}_{k+p+1}e_{k}^{T},M_{w}\end{pmatrix},$

Finally, equating the individual columns on both sides and discarding columns
$k+1,\ldots,k+p$ we are left with the new Krylov spaces of order $k$ and the
Lanczos relations

(10a) $\displaystyle
AV^{\prime\prime}_{k}-V^{\prime\prime}_{k}T^{\prime\prime}_{k}$
$\displaystyle=$ $\displaystyle\widetilde{v}^{\prime\prime}_{k+1}e_{k}^{T},$
(10b) $\displaystyle
A^{H}W^{\prime\prime}_{k}-W^{\prime\prime}_{k}(T^{\prime\prime}_{k})^{H}$
$\displaystyle=$ $\displaystyle\widetilde{w}^{\prime\prime}_{k+1}e_{k}^{T},$
(10c) $\displaystyle(W^{\prime\prime}_{k})^{H}V_{k}$ $\displaystyle=$
$\displaystyle I.$

The new residual vectors are defined as

(11a) $\displaystyle\widetilde{v}^{\prime\prime}_{k+1}$ $\displaystyle=$
$\displaystyle\beta_{k+1}V^{\prime\prime}_{p}e_{1}+\theta_{k+p}\widetilde{v}_{k+p+1},$
(11b) $\displaystyle\widetilde{w}^{\prime\prime}_{k+1}$ $\displaystyle=$
$\displaystyle\bar{\delta}_{k+1}W^{\prime\prime}_{p}e_{1}+\phi_{k+p}\widetilde{w}_{k+p+1},$

which may be normalized appropriately such that the inner product of the new
Krylov vectors is unity. The Lanczos process may now be carried out again to
generate the next $p$ vectors of the right and left Krylov spaces and the
cycle may be repeated till an adequately converged subspace has been obtained.
In the case that exact shifts (eigenvalues of $T_{m}$) are used for the
restart procedure $\beta_{k+1}$ and $\delta_{k+1}$ are both zero in exact
arithmetic.

As a final note, we mention that Della-Dora [3] introduced a class of
algorithms of the GR type of which, LR and QR are special cases. Watkins then
introduced generic bulge-chasing algorithms for the entire GR class of
methods[20]. The shifted LR algorithm used in the restart procedure outlined
above can therefore be carried out in an implicit manner through the bulge-
chase sequence of Watkins [20]. For real matrices, the operations can be
confined in the real space by using the double-shift strategy introduced by
Francis [4]. The entire restart process is put together in algorithm 1.

Algorithm 1 Restarted nonsymmetric Lanczos

Input:  $V_{k+p},\widetilde{v}_{k+p+1},W_{k+p},\widetilde{w}_{k+p+1},T_{k+p}$
Input: $\mu_{1},\mu_{2},\ldots\mu_{p}$ $\triangleright$ Unwanted shifts
Output:  $V_{k},\widetilde{v}_{k+1},W_{k},\widetilde{w}_{k+1},T_{k}$

1:$AV_{k+p}-V_{k+p}T_{k+p}=\widetilde{v}_{k+p+1}e_{k+p}^{T}$ $\triangleright$
Right Lanczos relation

2:$A^{H}W_{k+p}-W_{k+p}T^{H}_{k+p}=\widetilde{w}_{k+p+1}e_{k+p}^{T}$
$\triangleright$ Left Lanczos relation

3:$W_{k+p}^{H}V_{k+p}=I$; $W_{k+p}^{H}\widetilde{v}_{k+p+1}=0$;
$V_{k+p}^{H}\widetilde{w}_{k+p+1}=0$; $\triangleright$ Bi-orthogonality

4:procedure Refine right subspace

5: $L_{r}\leftarrow I_{k+p}$

6: for $j\leftarrow 1$ to $p$ do

7: $T_{k+p}\leftarrow L^{-1}T_{k+p}L$ $\triangleright$ Implicitly shifted LR
with shift $\mu_{j}$

8: $L_{r}\leftarrow LL_{r}$

9: end for

10:end procedure

11:procedure Refine left subspace

12: $L_{l}\leftarrow I_{k+p}$

13: for $j\leftarrow 1$ to $p$ do

14: $T^{H}_{k+p}\leftarrow L^{-1}T^{H}_{k+p}L$ $\triangleright$ Implicitly
shifted LR with shift $\bar{\mu}_{j}$

15: $L_{l}\leftarrow LL_{l}$

16: end for

17:end procedure

18:procedure Update

19: $\beta\leftarrow e^{T}_{k+1}T_{k+p}e_{k}$; $\delta\leftarrow
e^{T}_{k}T_{k+p}e_{k+1}$

20: $\theta=e_{k+p}^{T}(L_{r}(L_{l})^{-H})e_{k}$;
$\phi=e_{k+p}^{T}(L_{r}^{-H}L_{l})e_{k}$

21: $\widetilde{v}_{k+1}\leftarrow\beta
V_{k+p}e_{k+1}+\theta\widetilde{v}_{k+p+1}$ $\triangleright$ New right
residual

22:
$\widetilde{w}_{k+1}\leftarrow\bar{\delta}W_{k+p}e_{k+1}+\phi\widetilde{w}_{k+p+1}$
$\triangleright$ New left residual

23: $V_{k}\leftarrow V_{k+p}(L_{r}L_{l}^{-H})\begin{pmatrix}I_{k}\\\
0_{p}\end{pmatrix}$ $\triangleright$ New right Krylov space

24: $W_{k}\leftarrow W_{k+p}(L_{r}^{-H}L_{l})\begin{pmatrix}I_{k}\\\
0_{p}\end{pmatrix}$ $\triangleright$ New left Krylov space

25:
$T_{k}\leftarrow\begin{pmatrix}I_{k}&0_{p}\end{pmatrix}T_{k+p}\begin{pmatrix}I_{k}\\\
0_{p}\end{pmatrix}$ $\triangleright$ New tridiagonal matrix

26:end procedure

27:return $V_{k},\widetilde{v}_{k+1},W_{k},\widetilde{w}_{k+1},T_{k}$

##  3 Computational Results

We present the results of computational tests performed in the Grcar matrix
which is highly non-normal and has presented problems with convergence in
previous studies. In his implicit restart work Sorensen [18] indeed points out
that the restarted Arnoldi has trouble converging to the left-most part of the
spectrum and even though the algorithm claimed convergence, what was obtained
was in fact part of the pseudospectrum.

We present the results of the restarted Lanczos method applied to the
$50\times 50$ Grcar matrix for $k=10$ and $p=10$. Figure 2 shows the spectrum
obtained after performing eleven restarts of the Lanczos algorithm. The
restarts were performed such that the $k$ eigenvalues with the largest
imaginary part were retained for the reduced operator and the remaining were
used as shifts to be discarded from the subspace. The same part of the
spectrum was sought by Sorensen in [18]. In this particular case we obtained
$\widetilde{w}_{m+1}^{T}\widetilde{v}_{m+1}\sim O(10^{-11})$ even though the
individual residuals were both of order $10^{-3}-10^{-4}$. At this particular
moment, one would need to employ a look-ahead step of the Lanczos however, we
have shown the spectrum of the reduced matrix at this point. One may think of
the spectrum to have converged to $O(10^{-4})$. As seen from the figure, the
spectrum of the reduced operator matches quite well with the original spectrum
of the Grcar matrix. In figure 2, we show the convergence history of the
individual eigenvalues with each restart step as the Lanczos algorithm
progresses. Clearly the restart process is working well to shift the Krylov
spaces towards the wanted region of the spectrum. The error in the eigenvalues
is of $O(10^{-7})$ even though the perturbation to the reduced matrix is of
order $10^{-4}$. This is in contrast to the results obtained in [18] which
reported large perturbations to the eigenvalues even when the Arnoldi method
reported convergence. We suspect this is due to the fact that the non-
symmetric Lanczos approximates both the right and left eigenspaces
simultaneously and leads to a lower error in the truncated matrix. We expect
this to be particularly useful in hydrodynamic problems where highly non-
normal matrices are a routine occurrence and where problems with convergence
of the spectrum have often been reported (see [2] for example).

Figure 1: Spectrum obtained from the restarted Lanczos after $11$ restart
steps.

Figure 2: Error in the eigenvalues after each restart $i$.

At this point we make a note of caution that for large Grcar matrices our
attempts have been somewhat less successful and the convergence sometimes
fails. We expect this is due to the naive implementations of the implicit LR
method that we have done in Julia where, the standard checks for small sub-
diagonal elements have not been performed. In our investigations of such cases
we indeed do find small sub diagonal elements which lead to loss of precision.
Very small sub diagonal elements also lead to rapid loss of biorthogonality of
the two subspaces. Similarly we have not paid attention to the issue arising
out of the Lanczos algorithm itself, except for employing a double (two-sided)
Gram-Schmidt to ensure biorthogonality. We also note that for large number of
restarts, the biorthogonality property of the two subspaces is progressively
lost, probably due to accruing floating point errors. Hence we expect some
method of reorthogonalization would be required to ensure stability for very
long calculations, which has not been done in the current work. The (implicit)
LR decomposition is not unique. Uniqueness is ensured for unit main diagonal
of the lower triangular matrix however, this does not pay any heed to
conditioning of the transforming matrices. A better strategy of building the
matrices could be pursued which has better conditioning while at the same time
preserves the tridiagonal structure. These issues would require careful
implementation of all the individual components and we do not address those in
the current work.

Finally we note that very similar work has been reported in [7] where
hyperbolic transforms (HR) are used for restarts of the bi-Lanczos method, and
in [6] where the LR transformations are used, albeit in the context of model
reduction.

##  4 Conclusion

We present an algorithm to restart the non-symmetric Lanczos method which is
based on the idea of polynomial filtering via the implicit QR method proposed
by Sorensen [18] for restarting the Arnoldi iteration. The (implicitly)
shifted LR method is shown to be the appropriate analogue for restarting the
non-symmetric Lanczos algorithm for structure preserving restarts. It is shown
that the polynomial filtering process needs to be carried out for both the
right and the left Krylov spaces and the appropriate transforms for
maintaining biorthogonality of the two spaces are highlighted. Computational
results are shown for a Grcar matrix which is known to be highly non-normal
and the spectrum is found to converge adequately even when the residual error
is relatively large.

## Acknowledgments

The authors would like to thank Professor Elias Jarlebring for his helpful
comments on the manuscript. The authors acknowledge support of Nordita and the
Swedish Research Council Grant No. 2018-04290. Nordita is partially supported
by Nordforsk.

## References

  * [1] W. E. Arnoldi, The principle of minimized iterations in the solution of the matrix eigenvalue problem, Quarterly of applied mathematics, 9 (1951), pp. 17–29.
  * [2] M. Brynjell-Rahkola, N. Shahriari, P. Schlatter, A. Hanifi, and D. S. Henningson, Stability and sensitivity of a cross-flow-dominated falkner–skan–cooke boundary layer with discrete surface roughness, Journal of Fluid Mechanics, 826 (2017), p. 830–850.
  * [3] J. Della-Dora, Numerical linear algorithms and group theory, Linear Algebra and its Applications, 10 (1975), pp. 267–283.
  * [4] J. G. F. Francis, The QR Transformation—Part 2, The Computer Journal, 4 (1962), pp. 332–345.
  * [5] S. W. Gaaf and E. Jarlebring, The infinite bi-lanczos method for nonlinear eigenvalue problems, SIAM Journal on Scientific Computing, 39 (2017), pp. S898–S919.
  * [6] E. Grimme, D. Sorensen, and P. Van Dooren, Stable partial realizations via an implicitly restarted lanczos method, in Proceedings of 1994 American Control Conference - ACC ’94, vol. 3, 1994, pp. 2814–2818 vol.3.
  * [7] E. J. Grimme, D. C. Sorensen, and P. Van Dooren, Model reduction of state space systems via an implicitly restarted lanczos method, Numerical algorithms, 12 (1996), pp. 1–31.
  * [8] M. H. Gutknecht, Lanczos-type solvers for nonsymmetric linear systems of equations, Acta numerica, 6 (1997), pp. 271–397.
  * [9] V. Hernandez, J. E. Roman, and V. Vidal, SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems, ACM Trans. Math. Software, 31 (2005), pp. 351–362.
  * [10] C. Lanczos, An iteration method for the solution of the eigenvalue problem of linear differential and integral operators, (1950).
  * [11] R. B. Lehoucq, D. C. Sorensen, and C. Yang, ARPACK users’ guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods, SIAM, 1998.
  * [12] V. Papakos and I. Jaimoukha, A deflated implicitly restarted lanczos algorithm for model reduction, in 42nd IEEE International Conference on Decision and Control (IEEE Cat. No.03CH37475), vol. 3, 2003, pp. 2902–2907 Vol.3.
  * [13] B. N. Parlett, D. R. Taylor, and Z. A. Liu, A look-ahead lanczos algorithm for unsymmetric matrices, Mathematics of Computation, 44 (1985), pp. 105–124.
  * [14] H. Rutishauser, Solution of eigenvalue problems with the lr-transformation, National Bureau of Standards, Applied Mathematics Series, 49 (1958), pp. 47–81.
  * [15] H. Rutishauser, Lectures on Numerical Mathematics, Birkhauser, 1990\.
  * [16] H. Rutishauser and H. Schwarz, The lr transformation method for symmetric matrices, Numerische Mathematik, 5 (1963), pp. 273–289.
  * [17] Y. Saad, The lanczos biorthogonalization algorithm and other oblique projection methods for solving large unsymmetric systems, SIAM Journal on Numerical Analysis, 19 (1982), pp. 485–506.
  * [18] D. C. Sorensen, Implicit application of polynomial filters in a k-step arnoldi method, Siam journal on matrix analysis and applications, 13 (1992), pp. 357–385.
  * [19] G. W. Stewart, A krylov–schur algorithm for large eigenproblems, SIAM Journal on Matrix Analysis and Applications, 23 (2002), pp. 601–614.
  * [20] D. S. Watkins and L. Elsner, Chasing algorithms for the eigenvalue problem, SIAM Journal on Matrix Analysis and Applications, 12 (1991), pp. 374–384.