$\epsilon^{\mathrm{P}}$ can be estimated using single and two-qubit RB
information.

QubitError Rate (a) Average readout error. The readout error is the
probability the state of a given qubit is incorrectly labelled.

QubitError Rate (b) Average error per ${\sf U_{2}}$ gate. The error per gate
is a measure of how accurately the ${\sf U_{2}}$ gate is applied.

Figure 15: Error per single qubit operations on the devices used in this work.
Bars indicate the average error rates; error bars are one standard deviation.
Data aggregated based on calibration data collected over the course of our
experiments. Devices shown here are: ibmqx2 [ ], ibmq_ourense [ ],
ibmq_singapore [ ], ibmq_16_melbourne [ ]. A logarithmic scale is used. Error
RateQubits Figure 16: Average error per CX operation on the devices used in
this work. The error per CX gate is a measure of how accurately the CX gate is
applied. Bars indicate the average error rates; error bars are one standard
deviation. Data aggregated based on calibration data collected over the course
of our experiments. Devices shown here are: ibmqx2 [ ], ibmq_ourense [ ],
ibmq_singapore [ ], ibmq_16_melbourne [ ]. A logarithmic scale is used.

Several important noise channels, most notably cross-talk, are not included in
the device calibration data. As shown in Section 5, the effects of this noise
can be inferred through the application-motivated benchmarks we introduce in
this work, by showing the trade-off between connectivity of the device and
cross-talk [75].

##  Appendix D Empirical Relationship Between Heavy Output Generation
Probability and $\ell_{1}$-Norm Distance

As discussed in Section 2.4, the theoretical foundations for believing that
implementing shallow circuits to within a fixed $\ell_{1}$-norm distance
constitutes a demonstration of quantum computational supremacy are stronger
than for implementations with high heavy output generation probability. That
said, Figure 3 and Figure 11 contain similar features. For example,
ibmq_16_melbourne consistently performs the worst, with ibmq_singapore and
ibmq_ourense performing the best by both figures of merit. How then do these
two figures of merit generally relate to one another.

If the $\ell_{1}$-norm distance was $0$, the experimental outcome frequencies
would equal the ideal outcome probabilities. Consequently, the heavy output
probabilities would be the same between the device and an ideal quantum
computer. Because the heavy output probability depends on the circuit in
question, it is useful to normalise the latter by the heavy output probability
of an ideally-implemented circuit. We define _the normalised heavy output
generation probability_ as the ratio of the heavy output probability of the
device and the heavy output probability from an ideal quantum computer. Hence
if the $\ell_{1}$-norm distance was 0, the normalised heavy output generation
probability would be 1.

As the $\ell_{1}$-norm distance increases, the experimental frequencies
increasingly differ from the ideal outcome probabilities. Two things then may
happen: heavy outputs are produced more regularly, in which case the
normalised heavy output generation probability will grow above $1$; or less
regularly, in which case the normalised heavy output generation probability
will fall below $1$. In practice, we expect the distribution produced by the
device to converge to the uniform one as the noise increases, so we expect the
normalised heavy output generation probability to fall with increasing
$\ell_{1}$-norm distance.

Normalised Heavy Output Probability$\ell_{1}$-norm distanceCircuit
Classshallow circuitssquare circuitsdeep
circuitsDeviceibmq_ourenseibmq_16_melbourneibmqx2ibmq_singapore Figure 17:
Scatter plot and linear regression line comparing the normalised heavy output
generation probability and $\ell_{1}$-norm distance. Each point corresponds to
one circuit of the class and width as labelled. Colours correspond to numbers
of qubits in the following way: 2 [ ], 3 [ ], 4 [ ], 5 [ ], 6 [ ], 7 [ ].

The empirical relationship between the normalised heavy output generation
probability and $\ell_{1}$-norm distance is shown in Figure 17. For each
circuit, Figure 17 plots the $\ell_{1}$-norm distance of the distribution
produced by a real device against the normalised heavy output generation
probability. As expected, a negative correlation exists between these two
figures of merit. For the deepest circuits, and in particular the widest
circuits from the deep circuits class, the cluster of points indicates that
the the normalised heavy output generation probability falls more slowly as
the $\ell_{1}$-norm distance becomes larger. This is because the minimum value
of heavy output generation probability is being reached, which is to say that
the output distribution from the real device has converged to the uniform one,
while more detail can be extracted by considering the $\ell_{1}$-norm
distance.

This correlation is encouraging as, in the regime where it becomes impossible
to calculate the $\ell_{1}$-norm distance, we can be justified in believing
the correlation between the features present in Section 5 persist. This is to
say that we may be able to make lose predictions about the performance of a
quantum computing stack as measured by $\ell_{1}$-norm distance based on it’s
performance as measured by Heavy Output Generation Benchmarking. This
reasoning is similar to that used when Cross-Entropy Benchmarking is used to
predict demonstrations of quantum computational supremacy in the regime when
it becomes impossible to calculate [3].


# Sharp pointwise-in-time error estimate of L1 scheme for nonlinear
subdiffusion equations

Dongfang Li School of Mathematics and Statistics, Huazhong University of
Science and Technology, Wuhan, 430074, China. (dfli@mail.hust.edu.cn). The
research of this author was supported by the National Natural Science
Foundation of China under grants No. 11771162\.  Hongyu Qin School of
Mathematics and Statistics, Wuhan University, Wuhan, 430072, China.
(qinhongyuhust@sina.com).  Jiwei Zhang School of Mathematics and Statistics,
Wuhan University, Wuhan, 430072, China. (jiweizhang@whu.edu.cn). The research
of this author was supported by the National Natural Science Foundation of
China under grants No. 11771035\.

###### Abstract

An essential feature of the subdiffusion equations with the $\alpha$-order
time fractional derivative is the weak singularity at the initial time. The
weak regularity of the solution is usually characterized by a regularity
parameter $\sigma\in(0,1)\cup(1,2)$. Under this general regularity assumption,
we here obtain the pointwise-in-time error estimate of the widely used L1
scheme for nonlinear subdiffusion equations. To the end, we present a refined
discrete fractional-type Grönwall inequality and a rigorous analysis for the
truncation errors. Numerical experiments are provided to demonstrate the
effectiveness of our theoretical analysis.

###### keywords:

Sharp pointwise-in-time error estimate, L1 scheme, nonlinear subdiffusion
equations, non-smooth solutions

###### AMS:

65M06,65M12,65M15,35R11

##  1 Introduction

In this paper, we consider sharp pointwise-in-time error estimate of L1 scheme
in time for solving the following nonlinear subdiffusion equations

$\displaystyle~{}\partial_{t}^{\alpha}u-\Delta u=f(u),\quad
x\in\Omega\times(0,T]$ (1.1)

with the initial and boundary conditions

$\displaystyle\begin{array}[]{ll}u(x,0)=u_{0}(x),&x\in\Omega,\\\
u(x,t)=0,&x\in\partial\Omega\times[0,T],\end{array}$ (1.4)

where $\Omega=(0,L)^{d}\subset\mathbb{R}^{d}$ ($d\geq 1$). The time fractional
Caputo derivative is defined as

$~{}\partial_{t}^{\alpha}u(x,t)=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}\frac{\partial
u(x,s)}{\partial s}\frac{1}{(t-s)^{\alpha}}ds,\quad 0<\alpha<1.$ (1.5)

Here $\Gamma(\cdot)$ denotes the Gamma function. The equations provide a
powerful tool to describe anomalous diffusion in different physical
situations, see e.g., [1, 2, 3]. Hence, the theoretical and numerical analysis
of the models have attracted the interest of plenty of researchers.

In developing numerical methods for solving the subdiffusion problem (1.1), an
important consideration is that the solution $u$ is typically less regular
than in the case of a classical parabolic PDE (as the limiting case $\alpha\to
1$). For instance, Jin et al. [4] show that if the initial condition $u_{0}\in
H_{0}^{1}(\Omega)\cap H^{2}(\Omega)$, the solution to problem (1.1) satisfies
$\|\partial_{t}u(t)\|_{L^{2}(\Omega)}\leq Ct^{\alpha-1}$. Maskari and Karaa
[5] obtain that if $u_{0}\in\dot{H}^{\nu}(\Omega)$ with $\nu\in(0,2]$, the
solution of problem (1.1) satisfies $\|\partial_{t}u(t)\|_{L^{2}(\Omega)}\leq
Ct^{\nu\alpha/2-1}$, which implies that there exists a parameter
$\sigma\in(0,\alpha]$ and $u_{t}\rightarrow\infty$ as $t\rightarrow 0^{+}$.
One can refer to more works [6, 7, 8, 9, 10, 11, 12] on the discussion of the
regularity of solutions. With loss of generality, we assume the solution
regularity satisfies

$\|{\partial^{m}_{t}}u\|_{L^{2}(\Omega)}\leq Ct^{\sigma-m},\quad\text{ for
}m=1,2,\text{ and }\sigma\in(0,1)\cup(1,2].$ (1.6)

Under the regularity assumption $\sigma=\alpha$ in (1.6), many works indicate
that the convergence order with the maximum norm in time is
$\mathcal{O}(\tau^{\alpha})$, where $\tau$ is the temporal stepsize. Here we
refer readers to [13, 14, 15, 16] for the L1 and L2-type schemes on the
uniform meshes, [4, 17] for the the convolution quadrature (CQ) Euler method
and [4, 18] for the CQ BDF methods. In addition, numerical simulations show an
interesting phenomenon that the convergence order of the L1 scheme is
$\mathcal{O}(\tau^{\alpha})$ as $t$ tends to $0$, and $\mathcal{O}(\tau)$ at
the final time $t=T$. It motivates much works to consider the pointwise error
estimate. For linear subdiffusion equations (i.e., $f(u)=0$), Gracia et al.
[19] proved the temporal error of L1 scheme is of $\tau t_{n}^{\alpha-1}$. Yan
et el. [20] considered time-stepping error estimates of the modified L1
scheme. Jin et al. [21] showed if the initial condition $u_{0}(x)\in
L_{2}(\Omega)$, the temporal error of L1 scheme is of $\tau t_{n}^{-1}$. After
that, they further obtained [22] time-stepping error estimates of some high-
order BDF convolution quadrature methods. Mustapha and McLean [23, 24]
investigated time-stepping error bounds of discontinuous Galerkin methods for
fractional diffusion problems. For the nonlinear subdiffusion equations(1.1),
Maskari and Karaa [5] study the optimal pointwise-in-time error estimates
based on the CQ Euler method. As far as we know, it still remains open to
obtain the pointwise error estimate of L1 scheme for nonlinear subdiffsuion
equations.

The main goal of this paper is to present the sharp pointwise-in-time error
estimate of the widely used L1 scheme for the nonlinear subdiffusion equations
under the regularity assumption (1.6). Generally, our goals are theoretically
challenging mainly due to two reasons. On the one hand, the numerical Caputo
formula always has a form of discrete convolutional summation. The temporal
truncation error varies at each time level and has a close relation with the
time level $n$, which leads to the consistency analysis becomes too cumbersome
to implement in practice. On the other hand, due to the nonlinearity, the
evolutions of the solutions and errors are non-monotone decreasing, which
requires the rigorous and refined analysis step by step.

In this paper, we consider the L1 scheme on the uniform meshes to approximate
the time fractional derivative and central finite difference scheme with the
stepsize $h$ to discretize the diffusion term. We overcome the mentioned
difficulties and paint a full picture for sharp pointwise error estimate of
the L1 scheme for the nonlinear time fractional parabolic problems under
assumption (1.6). The given results indicate that

  * •

Suppose $\sigma\in(0,\alpha)$, one has

$\displaystyle\|u^{n}-U^{n}\|_{\infty}\lesssim\tau^{\sigma+1-\alpha}t_{n}^{\alpha-1}+t_{n}^{\alpha}h^{2},n=1,2,\cdots,N,$

where $u^{n}$ and $U^{n}$ are theoretical and numerical solutions. The result
has never been found before, even for the linear time fractional problems.

  * •

Suppose that $\sigma=\alpha$, one has

$\displaystyle\|u^{n}-U^{n}\|_{\infty}\lesssim\tau
t_{n}^{\alpha-1}+t_{n}^{\alpha}h^{2},n=1,2,\cdots,N.$

The pointwise error estimates are firstly given for the nonlinear problems.
Especially, when $f(u)=0$, problem (1.1) is reduced to the linear subdiffusion
models. Then, the given result agrees with the one in [19].

  * •

Suppose the solution is smoother (i.e., $\sigma>\alpha$), one has

$\displaystyle\|u^{n}-U^{n}\|_{\infty}\lesssim\begin{cases}\tau^{\sigma+1-\alpha}t_{n}^{\alpha-1}+t_{n}^{\alpha}h^{2},~{}~{}~{}~{}~{}~{}&\alpha<\sigma<1,\\\
\tau^{2-\alpha}t_{n}^{\alpha+\sigma-2}+t_{n}^{\alpha}h^{2},~{}~{}~{}~{}~{}~{}&1<\sigma\leq
2,\end{cases}n=1,2,\cdots,N.$

From the results, one can see that when $t_{n}\rightarrow 0$, the error
estimate is of $\tau^{\sigma}$. When $t_{n}$ is far away from $0$, different
convergence results can be found. Especially, the maximum error in the whole
domain is of $\tau^{\min(\sigma,2-\alpha)}$. The results on the maximum errors
agree with the error estimates in [13].

The rest of the paper is organized as follows. In Section 2, we present the
fully-discrete scheme and the time-stepping error estimates. In Section 3, we
present the discrete fractional Grönwall inequalities and a detailed proof of
the main results. In Section 4, we give several numerical examples to confirm
the theoretical results. Finally, we have the conclusions and discussions in
Section 5.

Throughout the paper, notation $A\lesssim B$ implies that there exists a
constant $c$ such that $A\leq cB$.

##  2 Fully discrete scheme and main convergence results

In this section, we present the fully-discrete scheme and the main convergence
results.

###  2.1 Fully discrete scheme

In this section, we will present the fully-discrete scheme and the main
convergence results. In this study, we always assume that $u(\cdot,t)\in
C^{4}(\Omega)$ for every $t$ so that the central finite difference method is
available to approximate the diffusion term.

Let $\mathcal{T}_{\tau}=\\{t_{n}|t_{n}=n\tau;\;0\leq n\leq N\\}$ be a uniform
partition of $[0,T]$ with the time step $\tau=T/N$. Let $h=L/M$ be the spatial
stepsize with $M$ a give positive integer. Denote
$u^{n}_{i_{1},i_{2},\cdots,i_{d}}=u(i_{1}h,i_{2}h,\cdots,i_{d}h,t_{n})$, where
$i_{j}=1,\cdots,M-1,\;j=1,\cdots,d$ and $n=0,1,\cdots,N$. At the grid point
$(i_{1}h,i_{2}h,\cdots,i_{d}h)$ , the central finite difference method is
given by

$\displaystyle\sum_{j=1}^{d}u_{x_{j}x_{j}}^{n}$ $\displaystyle=$
$\displaystyle\frac{1}{h^{2}}\sum_{j=1}^{d}(u^{n}_{i_{1},\cdots,x_{j}-1,\cdots,i_{d}}-2u^{n}_{i_{1},i_{2},\cdots,i_{d}}+u^{n}_{i_{1},\cdots,x_{j}+1,\cdots,i_{d}})+R^{n}_{i_{1},i_{2},\cdots,i_{d}}$
(2.1) $\displaystyle:=$
$\displaystyle\sum_{j=1}^{d}\delta_{x_{j}}^{2}u^{n}_{i_{1},i_{2},\cdots,i_{d}}+R^{n}_{i_{1},i_{2},\cdots,i_{d}},$

where $R^{n}_{i_{1},i_{2},\cdots,i_{d}}=\mathcal{O}(h^{2})$ is the spatial
truncation error.

The standard $L1$-approximation to the time fractional derivative is given by

$\displaystyle\partial_{t_{n}}^{\alpha}u$ $\displaystyle=$
$\displaystyle\frac{1}{\Gamma(1-\alpha)}\sum_{j=1}^{n}\frac{u^{j}_{i_{1},i_{2},\cdots,i_{d}}-u^{j-1}_{i_{1},i_{2},\cdots,i_{d}}}{\tau}\int_{t_{j-1}}^{t_{j}}\frac{1}{(t_{n}-s)^{\alpha}}ds+r^{n}_{i_{1},i_{2},\cdots,i_{d}}$
(2.2) $\displaystyle=$
$\displaystyle\frac{\tau^{-\alpha}}{\Gamma(2-\alpha)}\sum_{j=1}^{n}a_{n-j}(u^{j}_{i_{1},i_{2},\cdots,i_{d}}-u^{j-1}_{i_{1},i_{2},\cdots,i_{d}})+r^{n}_{i_{1},i_{2},\cdots,i_{d}}$
$\displaystyle:=$ $\displaystyle
D_{\tau}^{\alpha}u^{n}_{i_{1},i_{2},\cdots,i_{d}}+r^{n}_{i_{1},i_{2},\cdots,i_{d}},$

where $r^{n}$ represents the truncation error and
$a_{i}=(i+1)^{1-\alpha}-i^{1-\alpha},\;i\geq 0.$

Therefore, it holds that

$\displaystyle D_{\tau}^{\alpha}u^{n}_{i_{1},i_{2},\cdots,i_{d}}=$
$\displaystyle\sum_{j=1}^{d}\delta_{x_{j}}^{2}u^{n}_{i_{1},i_{2},\cdots,i_{d}}+f(u^{n}_{i_{1},i_{2},\cdots,i_{d}})+r^{n}_{i_{1},i_{2},\cdots,i_{d}}+R^{n}_{i_{1},i_{2},\cdots,i_{d}}.$
(2.3)

Let $U^{n}_{i_{1},i_{2},\cdots,i_{d}}$ be numerical approximation to
$u^{n}_{i_{1},i_{2},\cdots,i_{d}}$. Omitting the truncation errors in (2.3)
and replacing $u^{n}_{i_{1},i_{2},\cdots,i_{d}}$ by
$U^{n}_{i_{1},i_{2},\cdots,i_{d}}$, we get the fully-discrete scheme

$\displaystyle D_{\tau}^{\alpha}U^{n}_{i_{1},i_{2},\cdots,i_{d}}=$
$\displaystyle\sum_{j=1}^{d}\delta_{x_{j}}^{2}U^{n}_{i_{1},i_{2},\cdots,i_{d}}+f(U^{n}_{i_{1},i_{2},\cdots,i_{d}}).$
(2.4)

###  2.2 Convergence

The convergence analysis depends heavily on the consistency error and the
discrete fractional-type Grönwall inequality. For the purpose of the
readability, we present the lemmas here and leave the proof in the next
section.

We now present the truncated error of L1 scheme for Caputo derivative.

######  Lemma 2.1.

Suppose that $u$ satisfies (1.6). Then, it holds that

$\displaystyle|r_{n}|:=\left|D_{\tau}^{\alpha}u^{n}-\partial_{t_{n}}^{\alpha}u\right|\lesssim\begin{cases}\tau^{\sigma-\alpha}n^{-\min(1+\alpha,2-\sigma)},~{}~{}~{}~{}~{}&0<\sigma<1,\\\
\tau^{\sigma-\alpha}n^{-2+\sigma},&1<\sigma<2.\end{cases}$ (2.5)

The bound in (2.5) is sharp. It improves and generalizes the truncation error
bound proved in [14, Lemma 5.1]. Especially, when $\sigma=\alpha$, the bound
is consistent to the estimate in [19, Lemma 1].

We now consider a refined discrete fractional-type Grönwall inequality, which
will be used for convergence analysis (pointwise error esimate) of the L1
scheme studied in this paper.

######  Lemma 2.2 (A refined discrete fractional Grönwall inequality).

Suppose $0<\alpha<1$ and $\tau>0$. Let $y_{i}$, $0\leq i\leq N$, be a sequence
of non-negative real numbers satisfying

$\displaystyle D_{\tau}^{\alpha}y_{n}\leq\lambda
y_{n}+\mu_{1}n^{-\sigma_{1}}+\mu_{2}n^{-\sigma_{2}}+\eta,\quad\quad\
\text{for}~{}n=1,\dots,N,$ (2.6)

where $\sigma_{1}>1$,$\sigma_{2}<1$,$\lambda>0$, and $\mu_{1},\mu_{2}\geq 0$.
Then there exists a constant
$\tau_{*}=\sqrt[\alpha]{\frac{1}{2\Gamma(2-\alpha)\lambda}}$ such that
$\tau<\tau^{*}$, it holds

$\displaystyle y_{n}\lesssim
y_{0}+\mu_{1}\tau^{\alpha}n^{\alpha-1}+\mu_{2}\tau^{\alpha}n^{\alpha-\sigma_{2}}+\eta\tau^{\alpha}n^{\alpha}.$
(2.7)

######  Remark 2.1.

The case of $\mu_{1}=0$ and $\mu_{2}=0$ is investigated in [4, Theorem 2.6]
and [25, Lemma 3.1] and [13], respectively. The above Grönwall inequality is a
refined version, which is suitable to the pointwise error estimate for a class
of scheme.

In this paper, we mainly focus on the pointwise error estimate in time.
Without loss of generality, we assume the nonlinear term
$f:\mathbb{R}\rightarrow\mathbb{R}$ is Lipschitz continuous as

$\displaystyle|f(\xi_{1})-f(\xi_{2})|\leq
L|\xi_{1}-\xi_{2}|,~{}~{}\quad\textrm{for }\quad\xi_{1},\xi_{2}\in\mathbb{R},$
(2.8)

where $L$ denotes the Lipschitz coefficient. In this assumption, we have
following sharp pointwise error estimate of scheme (2.4).

######  Theorem 2.1.

Suppose that the system (1.1)-(1.4) has a unique solution satisfying (1.6).
Then, there exists a positive constant $\tau_{0}$, such that when
$\tau\leq\tau_{0}$, the finite difference system (2.4) admits a unique
solution $U^{n}$, satisfying, for $n=1,2,\cdots,N$, it holds

$\displaystyle\|u^{n}-U^{n}\|_{\infty}\lesssim\begin{cases}\tau^{\sigma+1-\alpha}t_{n}^{\alpha-1}+t_{n}^{\alpha}h^{2},~{}~{}~{}~{}~{}~{}&0<\sigma<1,\\\
\tau^{2-\alpha}t_{n}^{\alpha+\sigma-2}+t_{n}^{\alpha}h^{2},~{}~{}~{}~{}~{}~{}&1<\sigma\leq
2,\end{cases}$ (2.9)

where
$u^{n}=[u_{1,1,\cdots,1},u_{2,1,\cdots,1},\cdots,u_{M-1,1,\cdots,1},u_{1,2,\cdots,1},u_{2,2,\cdots,1},\cdots,u_{M-1,2,\cdots,1},\cdots,\\\
u_{1,1,\cdots,M-1},u_{2,1,\cdots,M-1},\cdots,u_{M-1,1,\cdots,M-1}]^{T}$ and
$U^{n}=[U_{1,1,\cdots,1},U_{2,1,\cdots,1},\cdots,\\\
U_{M-1,1,\cdots,1},U_{1,2,\cdots,1},U_{2,2,\cdots,1},\cdots,U_{M-1,2,\cdots,1},\cdots,U_{1,1,\cdots,M-1},U_{2,1,\cdots,M-1},\cdots,\\\
U_{M-1,1,\cdots,M-1}]^{T}.$

###### Proof.

We firstly take the proof of one-dimensional problem for an example. The exact
solution of problem (1.1) satisfies

$\displaystyle
D_{\tau}^{\alpha}u^{n}_{i}=\delta_{x}^{2}u_{i}^{n}+f(u_{i}^{n})+r^{n}_{i}+R_{i}^{n},$
(2.10)

where the temporal truncation error
$r^{n}_{i}=\partial_{t_{n}}^{\alpha}u-D_{\tau}^{\alpha}u^{n}_{i}$ and spatial
truncation error
$R_{i}^{n}=u_{xx}(x_{i},t^{n})-\delta_{x}^{2}u_{i}^{n}=\mathcal{O}(h^{2}).$

Suppose that $\|e^{n}\|_{\infty}=e^{n}_{i_{0}}$ and let
$e^{n}_{i}=u^{n}_{i}-U_{n}^{i}$. The error equation at the grid point
$(x_{i_{0}},t_{n})$ satisfies

$\displaystyle
D_{\tau}^{\alpha}e^{n}_{i_{0}}=\delta_{x}^{2}e_{i_{0}}^{n}+f(u^{n}_{i_{0}})-f(U^{n}_{i_{0}})+r^{n}_{i_{0}}+R_{i_{0}}^{n},$

which can be rewritten as

$\displaystyle\Big{(}\frac{\tau^{\\!-\\!\alpha}}{\Gamma(2\\!-\\!\alpha)}a_{0}\\!+\\!\frac{2}{h^{2}}\Big{)}e_{i_{0}}^{n}\\!$
$\displaystyle=$
$\displaystyle\\!\frac{e_{i_{0}\\!-\\!1}^{n}+e_{i_{0}+1}^{n}}{h^{2}}+\frac{\tau^{\\!-\\!\alpha}}{\Gamma(2\\!-\\!\alpha)}\sum_{j=1}^{n-1}(a_{n-j-1}-a_{n-j})e^{j}_{i_{0}}$
$\displaystyle+f(u^{n}_{i_{0}})-f(U^{n}_{i_{0}})+r^{n}_{i_{0}}+R_{i_{0}}^{n}.$

Therefore, we get

$\displaystyle\Big{(}\frac{\tau^{-\alpha}}{\Gamma(2-\alpha)}a_{0}+\frac{2}{h^{2}}\Big{)}|e_{i_{0}}^{n}|$
(2.11) $\displaystyle=$
$\displaystyle\Big{|}\frac{e_{i_{0}-1}^{n}+e_{i_{0}+1}^{n}}{h^{2}}+\frac{\tau^{-\alpha}}{\Gamma(2-\alpha)}\sum_{j=1}^{n-1}(a_{n-j-1}-a_{n-j})e^{j}_{i_{0}}+f(u^{n}_{i_{0}})-f(U^{n}_{i_{0}})+r^{n}_{i_{0}}+R_{i_{0}}^{n}\Big{|}$
$\displaystyle=$
$\displaystyle\Big{|}\frac{e_{i_{0}-1}^{n}+e_{i_{0}+1}^{n}}{h^{2}}+\frac{\tau^{-\alpha}}{\Gamma(2-\alpha)}\sum_{j=1}^{n-1}(a_{n-j-1}-a_{n-j})e^{j}_{i_{0}}+f(u^{n}_{i_{0}})-f(U^{n}_{i_{0}})+r^{n}_{i_{0}}+R_{i_{0}}^{n}\Big{|}$
$\displaystyle\leq$
$\displaystyle\Big{|}\frac{2e_{i_{0}}^{n}}{h^{2}}\Big{|}+\frac{\tau^{-\alpha}}{\Gamma(2-\alpha)}\sum_{j=1}^{n-1}(a_{n-j-1}-a_{n-j})|e^{j}_{i_{0}}|+L|e_{i_{0}}^{n}|+|r^{n}_{i_{0}}|+|R_{i_{0}}^{n}|,$

where we use the facts $a_{i}>a_{i+1}$.

Equation (2.11) further implies that

$\displaystyle\Big{(}\frac{\tau^{-\alpha}}{\Gamma(2-\alpha)}a_{0}\Big{)}|e_{i_{0}}^{n}|$
$\displaystyle\leq$
$\displaystyle\frac{\tau^{-\alpha}}{\Gamma(2-\alpha)}\sum_{j=1}^{n-1}(a_{n-j-1}-a_{n-j})|e^{j}_{i_{0}}|+L|e_{i_{0}}^{n}|+|r^{n}_{i_{0}}|+|R_{i_{0}}^{n}|,$

which can be rewritten as

$\displaystyle D_{\tau}^{\alpha}|e_{i_{0}}^{n}|\leq
L|e_{i_{0}}^{n}|+|r^{n}_{i_{0}}|+|R_{i_{0}}^{n}|.$ (2.12)

For $0\leq\sigma<1$, it holds that
$|r^{n}_{i_{0}}|\lesssim\tau^{\sigma-\alpha}n^{-\min(1+\alpha,2-\alpha)}$ with
$\min(1+\alpha,2-\alpha)>1$ (by lemma 3.2). Together with (2.12) and Lemma
2.2, we have

$\|e^{n}\|_{\infty}=|e_{i_{0}}^{n}|\lesssim\tau^{\sigma}n^{\alpha-1}+h^{2}=\tau^{\sigma+1-\alpha}t_{n}^{\alpha-1}+t_{n}^{\alpha}h^{2}.$

For $1<\sigma\leq 2$, it holds that
$|r^{n}_{i_{0}}|\lesssim\tau^{\sigma-\alpha}n^{\sigma-2}$ with $2-\sigma<1$
(by lemma 3.2). Together with (2.12) and Theorem 2.2, we have

$\|e^{n}\|_{\infty}=|e_{i_{0}}^{n}|\lesssim\tau^{\sigma}n^{\alpha+\sigma-2}+h^{2}=\tau^{2-\alpha}t_{n}^{\alpha+\sigma-2}+t_{n}^{\alpha}h^{2}.$

The extension to the multi-dimensional problems can be obtained similarly.
This completes the proof of the main results. ∎

######  Remark 2.2.

The convergence results imply that, when $t_{n}\rightarrow 0$, it holds that

$\displaystyle\max_{1\leq n\leq
N}\|u^{n}-U^{n}\|_{\infty}\lesssim\tau^{\sigma}+\tau^{\alpha}h^{2}.$ (2.13)

When $t$ is far away form $0$, it holds that

$\displaystyle\|u^{n}-U^{n}\|_{\infty}\lesssim\begin{cases}\tau^{\sigma+1-\alpha}+h^{2},~{}~{}~{}~{}~{}~{}&\sigma\in(0,1),\\\
\tau^{2-\alpha}+h^{2},~{}~{}~{}~{}~{}~{}&\sigma\in(1,2].\end{cases}$ (2.14)

Meanwhile, if we test the maximum error in the whole domain
$\Omega\times[0,T]$, we have

$\displaystyle\max_{1\leq n\leq
N}\|u^{n}-U^{n}\|_{\infty}\lesssim\begin{cases}\tau^{\sigma}+h^{2},~{}~{}~{}~{}~{}~{}&\sigma\in(0,1)\cup(1,2-\alpha),\\\
\tau^{2-\alpha}+h^{2},~{}~{}~{}~{}~{}~{}&\sigma\in[2-\alpha,2].\end{cases}$
(2.15)

######  Remark 2.3.

Scheme (2.4) is fully implicit. If the nonlinear term is approximated by using
the Newton linearized method, we have the following fully-discrete scheme,
i.e.,

$\displaystyle D_{\tau}^{\alpha}U^{n}_{i_{1},i_{2},\cdots,i_{d}}$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{d}\delta_{x_{j}}^{2}U^{n}_{i_{1},i_{2},\cdots,i_{d}}+f_{1}(U^{n-1}_{i_{1},i_{2},\cdots,i_{d}})$
(2.16)
$\displaystyle+f(U^{n-1}_{i_{1},i_{2},\cdots,i_{d}})(U^{n}_{i_{1},i_{2},\cdots,i_{d}}-U^{n-1}_{i_{1},i_{2},\cdots,i_{d}}),$

where $f_{1}(U^{n-1}_{i_{1},i_{2},\cdots,i_{d}})=\frac{\partial}{\partial
u}f|_{u=U^{n-1}_{i_{1},i_{2},\cdots,i_{d}}}$. We can have the same error
estimates as (2.9). The proof is similar to that done above. Meanwhile, the
proof of the convergence results is based on Lipschitz condition. If $f\in
C^{1}(\mathbb{R})$, the main results still hold. This is because

$\displaystyle\|U^{n}\|_{\infty}\leq\|u^{n}\|_{\infty}+\|u^{n}-U^{n}\|_{\infty}\lesssim\|u^{n}\|_{\infty}+1,$
(2.17)

whenever the stepsizes are sufficiently small. Due to the boundedness of
$\|U^{n}\|_{\infty}$, we have

$|f(u^{n}_{i_{1},i_{2},\cdots,i_{d}})-f(U^{n}_{i_{1},i_{2},\cdots,i_{d}})|\lesssim|u^{n}_{i_{1},i_{2},\cdots,i_{d}}-U^{n}_{i_{1},i_{2},\cdots,i_{d}}|.$

Then, the results can be proved by using similar method.

######  Remark 2.4.

The assumption $\sigma\in(0,\alpha]$ is reasonable and widely accepted. Here
we refer readers to [5, 18] for the detailed theoretical results. Suppose that
the solution is smoother, i.e., $\sigma>\alpha$, some additional hypothesis
should be added (see [14]). However, it is still quite restrictive due to the
strong hypotheisis. In this work, we assume that $\sigma\in(0,1)\cup(1,2)$ in
order to make the current analysis extendable.

######  Remark 2.5.

In the previous results, most time-stepping error estimates of L1 scheme focus
on the linear subdiffusion equations with reaction terms in the type of
$f(u)=-u$ and under the assumption $\sigma=\alpha$ [14, 19, 18] . For the
reaction terms in the type of $f(u)=u$ or nonlinear problems, the proof is
essentially different, which requires a refined discrete fractional-type
Grönwall inequality presented in Lemma 2.2 above.

######  Remark 2.6.

There are some error estimates of L1 scheme for nonlinear problems under the
assumption $\sigma=2$, e.g., [9, 25]. In such cases, the truncation error is
independent of $n$. Then, the convergence order is $2-\alpha$ and unchanged.
Such results can be concluded from the present results, but their proof is
totally different.

######  Remark 2.7.

In [13], another version of the discrete fractional-type Grönwall inequality
is developed. Based on the inequality and the regularity assumption (1.6), one
has the optimal error estimate $\mathcal{O}(\tau^{\sigma})$ for uniform time
mesh, and $\mathcal{O}(\tau^{\min\\{\gamma\sigma,2-\alpha\\}})$ for graded
mesh, where $\gamma$ represents the parameter reflecting the graded mesh
$t_{k}=T(k/N)^{\gamma}$. Note that in [13] there is a maximum in the discrete
fractional-type Grönwall inequality and the kernels are hard to explicitly be
expressed, it fails to obtain the pointwise error estimate.

##  3 Proof of the lemmas

In this section, we present a detailed proof of the lemmas in the previous
section.

###  3.1  Proof of Lemma 2.1

###### Proof.

We now estimate the truncation errors (2.5) in two cases.

Case A: n=1. It holds that

$\displaystyle\Big{|}D_{\tau}^{\alpha}u^{1}-\partial_{t_{1}}^{\alpha}u\Big{|}$
(3.1) $\displaystyle=$
$\displaystyle\Big{|}\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t_{1}}(t_{1}-s)^{-\alpha}[u(x,\xi_{1})-\frac{\partial}{\partial
s}u(x,s)]ds\Big{|}~{}~{}~{}\quad(\xi_{1}\in(0,t_{1}))$ $\displaystyle\leq$
$\displaystyle\left|\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t_{1}}(t_{1}-s)^{-\alpha}\frac{\partial}{\partial\xi_{1}}u(x,\xi_{1})ds\right|+\left|\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t_{1}}(t_{1}-s)^{-\alpha}\frac{\partial}{\partial
s}u(x,s)ds\right|$ $\displaystyle\lesssim$
$\displaystyle\left|\int_{0}^{t_{1}}(t_{1}-s)^{-\alpha}s^{\sigma-1}ds\right|$
$\displaystyle=$
$\displaystyle\left|t_{1}^{\sigma-\alpha}\beta(1-\alpha,\sigma)\right|$
$\displaystyle\lesssim$ $\displaystyle\tau^{\sigma-\alpha}.$

Case B: $n>1$. It follows from the definitions of the discrete and continuous
fractional operator that

$\displaystyle D_{\tau}^{\alpha}u^{n}-\partial_{t_{n}}^{\alpha}u$
$\displaystyle=\sum\limits_{k=0}^{n\\!-\\!1}\frac{1}{\Gamma(1\\!-\\!\alpha)}\int_{t_{k}}^{t_{k\\!+\\!1}}(t_{n}\\!-\\!s)^{\\!-\\!\alpha}[\frac{u^{k\\!+\\!1}\\!-\\!u^{k}}{\tau}-\frac{\partial}{\partial
s}u(x,s)]ds$
$\displaystyle=\sum\limits_{k=0}^{n\\!-\\!1}\frac{\\!-\\!\alpha}{\Gamma(1\\!-\\!\alpha)}\int_{t_{k}}^{t_{k+1}}(t_{n}\\!-\\!s)^{\\!-\\!\alpha\\!-\\!1}[\frac{u^{k+1}-u^{k}}{\tau}(s-t_{k})-(u(x,s)-u^{k})]ds$
$\displaystyle:$ $\displaystyle=\sum\limits_{k=0}^{n-1}R_{nk}.$ (3.2)

For $1\leq k\leq n-1$, by standard interpolation theory, it holds that

$\displaystyle\left|R_{nk}\right|$
$\displaystyle\lesssim\tau^{2}(\max\limits_{s\in[t_{k},t_{k+1}]}\left|u_{tt}(x,s)\right|\int_{t_{k}}^{t_{k+1}}(t_{n}-s)^{-\alpha-1}ds)$
$\displaystyle\lesssim\tau^{3}t_{k}^{\sigma-2}(t_{n}-t_{k+1})^{-\alpha-1}$
$\displaystyle=\tau^{\sigma-\alpha}k^{\sigma-2}(n-k+1)^{-\alpha-1}.$

As a result, we get

$\displaystyle\sum\limits_{k=1}^{\lceil n/2\rceil-1}\left|R_{nk}\right|$
$\displaystyle\lesssim\tau^{\sigma-\alpha}n^{-(1+\alpha)}\sum\limits_{k=1}^{\lceil\frac{n}{2}\rceil-1}k^{\sigma-2}\lesssim\begin{cases}\tau^{\sigma-\alpha}n^{-(1+\alpha)},~{}~{}~{}~{}~{}~{}~{}&0<\sigma<1,\\\
\tau^{\sigma-\alpha}n^{-2-\alpha+\sigma},~{}~{}~{}~{}&1<\sigma\leq
2.\end{cases}$ (3.3)

For $\lceil n/2\rceil\leq k\leq n-1$, it holds that

$\displaystyle\left|R_{nk}\right|$
$\displaystyle\lesssim\tau^{2}t_{k}^{\sigma-2}\int_{t_{k}}^{t_{k+1}}(t_{n}-s)^{-\alpha-1}ds$
$\displaystyle=\tau^{\sigma}k^{\sigma-2}\int_{t_{k}}^{t_{k+1}}(t_{n}-s)^{-\alpha-1}ds$
$\displaystyle\lesssim\tau^{\sigma}n^{\sigma-2}\int_{t_{k}}^{t_{k+1}}(t_{n}-s)^{-\alpha-1}ds.$
(3.4)

Therefore,

$\displaystyle\sum\limits_{k=\lceil
n/2\rceil}^{n-2}\left|R_{nk}\right|\lesssim\tau^{\sigma}n^{\sigma-2}\int_{t_{\ulcorner
n/2\rceil}}^{t_{n-2}}(t_{n}-s)^{-\alpha-1}ds\lesssim\tau^{\sigma-\alpha}n^{\sigma-2}.$
(3.5)

For $k=0$, noting that

$\displaystyle\left|\frac{u^{1}-u^{0}}{\tau}(s-0)-(u(x,s)-u^{0})\right|\leq
2\int_{0}^{t_{1}}\left|u_{t}(x,s)\right|ds\leq 2\tau^{\sigma},$

we get

$\displaystyle\left|R_{n0}\right|$
$\displaystyle=\left|\frac{-\alpha}{\Gamma(1-\alpha)}\int_{0}^{t_{1}}(t_{n}-s)^{-\alpha-1}[\frac{u^{1}-u^{0}}{\tau}(s-0)-(u(x,s)-u^{0})]ds\right|$
$\displaystyle\lesssim\tau^{\sigma}\int_{0}^{t_{1}}(t_{n}-s)^{-\alpha-1}ds$
$\displaystyle=\tau^{\sigma}\tau(t_{n}-t_{1})^{-\alpha-1}$
$\displaystyle=\tau^{\sigma-\alpha}n^{-\alpha-1}.$ (3.6)

For $k=n-1$, it holds that

$\displaystyle\left|T_{n,n-1}\right|$
$\displaystyle=\frac{1}{\Gamma(1-\alpha)}\int_{t_{n-1}}^{t_{n}}(t_{n}-s)^{-\alpha}[\frac{u^{n}-u^{n-1}}{\tau}-\frac{\partial}{\partial
s}u(x,s)]ds$
$\displaystyle\leq\frac{1}{\Gamma(1-\alpha)}\Big{|}\frac{\partial}{\partial
s}u(x,\xi_{n-1}^{*})-\frac{\partial}{\partial
s}u(x,\xi_{n-1}^{**})\Big{|}\int_{t_{n-1}}^{t_{n}}(t_{n}-s)^{-\alpha}ds$
$\displaystyle\lesssim\tau\Big{|}\frac{\partial^{2}}{\partial
s^{2}}u(x,\xi_{n-1}^{***})\Big{|}\int_{t_{n-1}}^{t_{n}}(t_{n}-s)^{-\alpha}ds$
$\displaystyle\leq\tau^{2-\alpha}t_{n-1}^{\sigma-2}$
$\displaystyle=\tau^{\sigma-\alpha}n^{\sigma-2},$ (3.7)

where $\xi_{n-1}^{*},\xi_{n-1}^{**}$ and $\xi_{n-1}^{***}$ are constants
between $t_{n-1}$ and $t_{n}$.

Now, together with the estimates (3.1)-(3.1), we obtain the finial results. ∎

###  3.2 Discrete complementary convolutions

In order to prove the refined discrete fractional-type Grönwall inequality, we
first introduce the useful properties of the discrete complementary
convolutions. They are given in [13] can be reduced to the following lemma by
removing the factor with respect to $\tau$.

######  Lemma 3.1 ([19]).

Let $\\{p_{n}\\}$ be a sequence defined by

$\displaystyle p_{0}=1,\quad p_{n}=\sum_{j=1}^{n}(a_{j-1}-a_{j})p_{n-j},\quad
n\geq 1.$ (3.8)

Then it holds that

  * (i)

$p_{n}<(n+1)^{\alpha-1}$, $1\leq k\leq n$,

  * (ii)

$\sum_{j=k}^{n}p_{n-j}a_{j-k}=1$, $1\leq k\leq n$,

  * (iii)

$\Gamma(2-\alpha)\sum_{j=1}^{n}p_{n-j}\leq\frac{n^{\alpha}}{\Gamma(1+\alpha)}$.

######  Lemma 3.2.

Let $\\{p_{n}\\}$ be a sequence defined in (3.8). Then, it holds that

$\displaystyle\frac{1}{\Gamma(1-\gamma)}\sum_{j=1}^{n}p_{n-j}j^{-\gamma}\lesssim\begin{cases}\frac{n^{\alpha-1}}{\Gamma(1-\gamma)},~{}~{}~{}~{}~{}&\gamma>1,\\\
\frac{n^{\alpha-\gamma}}{\Gamma(1-\gamma+\alpha)},&\gamma<1.\end{cases}$ (3.9)

###### Proof.

It follows from ((i)) in Lemma 3.1 and Lemma 2.1 that

$\displaystyle\sum_{j=1}^{n}p_{n-j}j^{-\gamma}$ $\displaystyle\lesssim$
$\displaystyle\sum_{j=1}^{n-1}(n-j+1)^{\alpha-1}j^{-\gamma}+n^{-\gamma}$
$\displaystyle\lesssim$
$\displaystyle(\frac{n}{2})^{\alpha-1}\sum_{j=1}^{\lceil
n/2\rceil}j^{-\gamma}+(\frac{n}{2})^{\alpha-\gamma}\sum_{j={\lceil
n/2\rceil}+1}^{n-1}(n-j+1)^{\alpha-1}j^{-\gamma}+\tau^{\sigma-\alpha}n^{-\gamma}$
$\displaystyle\lesssim$
$\displaystyle\Big{(}n^{\alpha-1}+n^{\alpha-\gamma}\int_{0}^{n}s^{-\alpha}(n-s)^{\alpha-1}ds+n^{-\gamma}\Big{)}$
$\displaystyle\lesssim$ $\displaystyle n^{\alpha-1},$

where we have noted that

$\sum_{j=1}^{n}j^{-\gamma}<\infty,\quad\gamma>1.$

For $\gamma<1$, it holds that

$\displaystyle\frac{1}{\Gamma(1-\gamma)}\sum_{j=1}^{n}p_{n-j}j^{-\gamma}$
$\displaystyle\lesssim$
$\displaystyle\frac{1}{\Gamma(1-\gamma)}\sum_{j=1}^{n-1}(n-j+1)^{\alpha-1}j^{-\gamma}+\frac{1}{\Gamma(1-\gamma)}n^{-\gamma}$
(3.10) $\displaystyle\lesssim$
$\displaystyle\frac{1}{\Gamma(1-\gamma)}\sum_{j=1}^{n-1}\int_{j-1}^{j}s^{-\gamma}(n+1-s)^{\alpha-1}+\frac{1}{\Gamma(1-\gamma)}n^{-\gamma}$
$\displaystyle\lesssim$
$\displaystyle\frac{1}{\Gamma(1-\gamma)}\int_{0}^{n}s^{-\gamma}(n-s)^{\alpha-1}+\frac{1}{\Gamma(1-\gamma)}n^{-\gamma}$
$\displaystyle\lesssim$
$\displaystyle\frac{n^{\alpha-\gamma}}{\Gamma(1-\gamma)}\beta{(1-\gamma,\alpha)}+\frac{1}{\Gamma(1-\gamma)}n^{-\gamma}$
$\displaystyle\lesssim$
$\displaystyle\frac{n^{\alpha-\gamma}}{\Gamma(1-\gamma+\alpha)}.$

This completes the proof. ∎

######  Lemma 3.3.

Let
$Z_{1}=(n^{-\sigma_{1}},(n-1)^{-\sigma_{1}},\cdots,1^{-\sigma_{1}})^{T},~{}Z_{2}=(n^{-\sigma_{2}},(n-1)^{-\sigma_{2}},\cdots,1^{-\sigma_{2}})^{T},Z_{3}=(1,1,\cdots,1)^{T}\in
R^{n}$ with $\sigma_{1}>1$ and $\sigma_{2}<1$, and

$J=2\Gamma(2-\alpha)\lambda\tau^{\alpha}\left[\begin{matrix}0&~{}p_{1}&~{}\cdots&p_{n-2}&p_{n-1}\\\
0&~{}0&~{}\cdots&~{}p_{n-3}&~{}p_{n-2}\\\
\vdots&~{}\vdots&~{}\ddots&~{}\vdots&~{}\vdots\\\
0&~{}0&~{}\cdots&~{}0&~{}p_{1}\\\ 0&~{}0&~{}\cdots&~{}0&~{}0\\\
\end{matrix}\right]_{n\times n}.$ (3.11)

Then, it holds that

  * (i)

$J^{i}=0,~{}~{}i\geq n$;

  * (ii)

$J^{m}Z_{2}\leq\frac{\Gamma(1-\sigma_{2})(2\Gamma(2-\alpha)\lambda\tau^{\alpha})^{m}}{\Gamma(1-\sigma_{2}+m\alpha)}\Big{(}n^{m\alpha-\sigma_{2}},(n-1)^{m\alpha-\sigma_{2}},\cdots,1^{m\alpha-\sigma_{2}}\Big{)}^{T}$,
$m=0,1,2,\cdots$;

  * (iii)

$\sum\limits_{j=1}^{i}J^{j}Z_{2}\lesssim\tau^{\alpha}\Big{(}n^{\alpha-\sigma_{2}}E_{\alpha,1-\sigma_{2}}(2\Gamma(2-\alpha)\lambda
t_{n}^{\alpha}),(n-1)^{\alpha-\sigma_{2}}E_{\alpha,1-\sigma_{2}}(2\Gamma(2-\alpha)\lambda
t_{n-1}^{\alpha}),\cdots,\\\
1^{\alpha-\sigma_{2}}E_{\alpha,1-\sigma_{2}}(2\Gamma(2-\alpha)\lambda
t_{1}^{\alpha}))\Big{)}^{T}$, $i\geq n$,

  * (iv)

$JZ_{1}\lesssim
2\Gamma(2-\alpha)\lambda\tau^{\alpha}\Big{(}n^{\alpha-1},(n-1)^{\alpha-1},\cdots,1^{\alpha-1}\Big{)}^{T}$;

  * (v)

$\sum\limits_{j=1}^{\infty}J^{j}Z_{1}\lesssim
2\Gamma(2-\alpha)\lambda\tau^{\alpha}\Big{(}n^{\alpha-1}E_{\alpha,\alpha}(2\Gamma(2-\alpha)\lambda
t_{n}^{\alpha}),(n-1)^{\alpha-1}E_{\alpha,\alpha}(2\Gamma(2-\alpha)\lambda
t_{n-1}^{\alpha}),\cdots,1^{\alpha-1}E_{\alpha,\alpha}(2\Gamma(2-\alpha)\lambda
t_{1}^{\alpha}))\Big{)}^{T}$ , $i\geq n$;

  * (vi)

$\sum\limits_{j=1}^{\infty}J^{j}Z_{3}\lesssim\tau^{\alpha}\Big{(}n^{\alpha}E_{\alpha,1-\sigma_{2}}(2\Gamma(2-\alpha)\lambda
t_{n}^{\alpha}),(n-1)^{\alpha}E_{\alpha,1-\sigma_{2}}(2\Gamma(2-\alpha)\lambda
t_{n-1}^{\alpha}),\cdots,\\\
1^{\alpha}E_{\alpha,1-\sigma_{2}}(2\Gamma(2-\alpha)\lambda
t_{1}^{\alpha}))\Big{)}^{T}$.

###### Proof.

Since that $J$ is an upper triangular matrix, one can check (i) holds .

Now, we prove (ii) by using the mathematical induction. Firstly, (ii) holds
for $m=0$. Suppose that (ii) holds for $m=k$.

$\displaystyle J^{k+1}Z_{2}$ $\displaystyle=$ $\displaystyle
J(J^{k}Z_{2})\leq\frac{\Gamma(1-\sigma_{2})(2\Gamma(2-\alpha)\lambda\tau^{\alpha})^{k}}{\Gamma(1-\sigma_{2}+k\alpha)}J\Big{(}n^{k\alpha-\sigma_{2}},(n-1)^{k\alpha-\sigma_{2}},\cdots,1^{k\alpha-\sigma_{2}}\Big{)}^{T}$
$\displaystyle=$
$\displaystyle\frac{\Gamma(1-\sigma_{2})(2\Gamma(2-\alpha)\lambda\tau^{\alpha})^{k+1}}{\Gamma(1-\sigma_{2}+k\alpha)}\Big{(}\sum_{i=1}^{n-1}p_{n-i}i^{1-\sigma_{2}+k\alpha},\sum_{i=1}^{n-2}p_{n-1-i}i^{1-\sigma_{2}+k\alpha},\cdots,0\Big{)}^{T}$
$\displaystyle\leq$
$\displaystyle\frac{\Gamma(1-\sigma_{2})(2\Gamma(2-\alpha)\lambda\tau^{\alpha})^{k+1}}{\Gamma(1-\sigma_{2}+k\alpha+\alpha)}\Big{(}n^{k\alpha-\sigma_{2}+\alpha},(n-1)^{k\alpha-\sigma_{2}+\alpha},\cdots,1^{k\alpha-\sigma_{2}+\alpha}\Big{)}^{T},$

where in the last inequality, we have noted Lemma 3.2. Therefore, (ii) holds
for $m=k+1$. This completes the mathematical induction and the proof.

Next, we show (iii) holds. Note that (i) implies that
$\sum_{j=1}^{i}J^{j}Z_{2}=\sum_{j=1}^{n-1}J^{j}Z_{2}$ for $i\geq n$, and by
(iii), we get

$\displaystyle\sum_{j=1}^{n-1}J^{j}Z_{2}$ $\displaystyle\lesssim$
$\displaystyle\tau^{\alpha}\sum_{j=1}^{n-1}\Big{(}\frac{(2\Gamma(2-\alpha)\lambda\tau^{\alpha}n^{\alpha})^{j}}{\Gamma(1-\sigma_{2}+j\alpha)}n^{\alpha-\sigma_{2}},\frac{(2\Gamma(2-\alpha)\lambda\tau^{\alpha}(n-1)^{\alpha})^{j}}{\Gamma(1-\sigma_{2}+j\alpha)}(n-1)^{\alpha-\sigma_{2}},$
$\displaystyle\cdots,\frac{(2\Gamma(2-\alpha)\lambda\tau^{\alpha}1^{\alpha})^{j}}{\Gamma(1-\sigma_{2}+j\alpha)}1^{\alpha-\sigma_{2}}\Big{)}^{T}$
$\displaystyle\lesssim$
$\displaystyle\tau^{\alpha}\Big{(}n^{\alpha-\sigma_{2}}E_{\alpha,1-\sigma_{2}}(2\Gamma(2-\alpha)\lambda
t_{n}^{\alpha}),(n-1)^{\alpha-\sigma_{2}}E_{\alpha,1-\sigma_{2}}(2\Gamma(2-\alpha)\lambda
t_{n-1}^{\alpha}),\cdots,$ $\displaystyle
1^{\alpha-\sigma_{2}}E_{\alpha,1-\sigma_{2}}(2\Gamma(2-\alpha)\lambda
t_{1}^{\alpha}))\Big{)}^{T}$

where

$E_{\alpha,\beta}(z)=\sum_{k=0}^{\infty}\frac{z^{k}}{\Gamma(k\alpha+\beta)}.$

For inequality (iv), one can check it holds by using Lemma 3.2.

For inequality (v), we have

$\sum\limits_{j=1}^{\infty}J^{j}Z_{1}=\sum\limits_{j=0}^{\infty}J^{j}(JZ_{1})\lesssim
2\Gamma(2-\alpha)\lambda\tau^{\alpha}\sum\limits_{j=0}^{\infty}J^{j}\Big{(}n^{\alpha-1},(n-1)^{\alpha-1},\cdots,1^{\alpha-1}\Big{)}^{T}.$

Then, the result (v) holds by using (iii) and the fact $1-\alpha<1$.

Finally, let $\sigma_{2}=0$ in (iii), one can obtain (vi). This completes the
proof. ∎

###  3.3 Proof of the Grönwall inequality

Now, we are ready to prove the refined Grönwall inequality in Lemma 2.2.

###### Proof.

It follows from the definition of $L1$-approximation (2.2) that

$\sum_{k=1}^{j}a_{j-k}{\delta_{t}y_{k}}\leq\Gamma(2-\alpha)\tau^{\alpha}\lambda_{1}y_{j}+\Gamma(2-\alpha)\tau^{\alpha}(\mu_{1}j^{-\sigma_{1}}+\mu_{1}j^{-\sigma_{2}}+\eta).$
(3.13)

Multiplying (3.13) by $p_{n-j}$ and then summing over for $j$ from $1$ to $n$,
we get

$\displaystyle\sum_{j=1}^{n}p_{n\\!-\\!j}\sum_{k=1}^{j}a_{j\\!-\\!k}{\delta_{t}y_{k}}\\!\leq\\!\Gamma(2\\!-\\!\alpha)\tau^{\alpha}\sum_{j=1}^{n}p_{n\\!-\\!j}\lambda_{1}y_{j}+\Gamma(2\\!-\\!\alpha)\tau^{\alpha}\sum_{j=1}^{n}p_{n\\!-\\!j}(j^{-\sigma_{1}}+j^{-\sigma_{2}}+\eta).$
(3.14)

Applying ((ii)) in Lemma 3.1 gives

$\displaystyle\sum_{j=1}^{n}p_{n-j}\sum_{k=1}^{j}a_{j-k}{\delta_{t}y_{k}}=\sum_{k=1}^{n}\delta_{t}y_{k}\sum_{j=k}^{n}p_{n-j}a_{j-k}=\sum_{k=1}^{n}\delta_{t}y_{k}=y_{n}-y_{0},\quad
n\geq 1.$ (3.15)

Substituting (3.15) into (3.14) gives

$\displaystyle y_{n}$ $\displaystyle\leq
y_{0}+\Gamma(2\\!-\\!\alpha)\tau^{\alpha}\sum_{j=1}^{n}p_{n\\!-\\!j}\lambda_{1}y_{j}+\Gamma(2\\!-\\!\alpha)\tau^{\alpha}\sum_{j=1}^{n}p_{n\\!-\\!j}(j^{-\sigma_{1}}+j^{-\sigma_{2}}+\eta)$
$\displaystyle=y_{0}+\Gamma(2-\alpha)\tau^{\alpha}y_{n}+\Gamma(2-\alpha)\tau^{\alpha}\lambda_{1}\sum_{j=1}^{n-1}p_{n-j}y_{j}$
$\displaystyle~{}~{}~{}+\Gamma(2\\!-\\!\alpha)\tau^{\alpha}\sum_{j=1}^{n}p_{n\\!-\\!j}(j^{-\sigma_{1}}+j^{-\sigma_{2}}+\eta)$
$\displaystyle\leq
y_{0}+\frac{1}{2}y_{n}+\Gamma(2-\alpha)\tau^{\alpha}\lambda_{1}\sum_{j=1}^{n-1}p_{n-j}y_{j}+\Gamma(2\\!-\\!\alpha)\tau^{\alpha}\sum_{j=1}^{n}p_{n\\!-\\!j}(j^{-\sigma_{1}}+j^{-\sigma_{2}}+\eta),$
(3.16)

whenever $\tau\leq\sqrt[\alpha]{\frac{1}{2\Gamma(2-\alpha)\lambda_{1}}}$.

Therefore, we get

$\displaystyle
y_{n}\\!\leq\\!2y_{0}\\!+\\!2\Gamma(2\\!-\\!\alpha)\tau^{\alpha}\lambda_{1}\sum_{j=1}^{n-1}p_{n-j}y_{j}\\!+\\!\Gamma(2\\!-\\!\alpha)\tau^{\alpha}\sum_{j=1}^{n}p_{n\\!-\\!j}(j^{-\sigma_{1}}\\!+\\!j^{-\sigma_{2}}\\!+\\!\eta).$
(3.17)

Let $Y=(y_{n},y_{n-1},\cdots,y_{1})^{T}$ and
$Z=C\Gamma(2\\!-\\!\alpha)\tau^{\alpha}(\mu_{1}Z_{1}+\mu_{2}Z_{2}+\eta
Z_{3})$. Then, (3.17) can be rewritten in a matrix form as follows

$Y\leq JY+JZ.$ (3.18)

As a result, we have

$\displaystyle Y$ $\displaystyle\leq$ $\displaystyle JV+JZ\leq
J(JY+JZ)+JZ=J^{2}Y+\sum_{j=1}^{2}J^{j}Z$ (3.19) $\displaystyle\leq$
$\displaystyle\cdots\leq
J^{n}Y+\sum_{j=1}^{n}J^{j}Z=\mu_{1}\sum_{j=1}^{n-1}J^{j}Z_{1}+\mu_{2}\sum_{j=1}^{n-1}J^{j}Z_{2}+\eta\sum_{j=1}^{n-1}J^{j}Z_{3}\;.$

Now, by (iv),(v) and (vi) in Lemma 3.3, we obtain (2.7) and complete the
proof. ∎

##  4 Numerical examples

In this section, we present several numerical results to illustrate the sharp
pointwise error estimates.
Example 1. Consider the one-dimensional nonlinear subdiffusion problems

$\displaystyle~{}\partial_{t}^{\alpha}u=u_{xx}+\sqrt{1+u^{2}}+g(x,t),\quad(x,t)\in(0,\pi)\times(0,t_{N}],$
(4.1)

where the initial condition and $g(x,t)$ are specially chosen such that the
problem admits an exact solution in the form of

$u(x,t)=t^{\sigma}\sin(x).$

Table 1: Maximum errors at $t=1$ and convergence orders in temporal direction with $M=1000$ for Example 1 $\alpha$ |  $\sigma$ $\backslash N$ | $10$ | $20$ | $40$ | $80$ | $160$ | $Rate$ | expected order
---|---|---|---|---|---|---|---|---
0.4 | $0.1$ | 3.50e-2 | 2.05e-2 | 1.22e-2 | 7.37e-3 | 4.46e-3 | 0.73 | $\sigma+1-\alpha$
| $0.4$ | 1.12e-2 | 5.21e-3 | 2.46e-3 | 1.18e-3 | 5.71e-4 | 1.07 | $\sigma+1-\alpha$
| $0.6$ | 5.13e-3 | 2.10e-3 | 8.70e-4 | 3.63e-4 | 1.53e-4 | 1.25 | $\sigma+1-\alpha$
| $1.2$ | 3.98e-3 | 1.57e-3 | 6.17e-4 | 2.41e-4 | 9.40e-5 | 1.36 | $2-\alpha$
| $1.8$ | 1.30e-2 | 5.03e-3 | 1.93e-3 | 7.38e-4 | 2.81e-4 | 1.39 | $2-\alpha$
0.6 | $0.4$ | 3.27e-2 | 1.78e-2 | 9.89e-3 | 5.54e-3 | 3.13e-3 | 0.84 | $\sigma+1-\alpha$
| $0.6$ | 1.51e-2 | 7.28e-3 | 3.54e-3 | 1.73e-3 | 8.51e-4 | 1.03 | $\sigma+1-\alpha$
| $0.8$ | 5.68e-3 | 2.49e-3 | 1.09e-3 | 4.76e-4 | 2.08e-4 | 1.19 | $\sigma+1-\alpha$
| $1.2$ | 3.98e-3 | 1.57e-3 | 6.17e-4 | 2.41e-4 | 9.40e-5 | 1.36 | $2-\alpha$
| $1.8$ | 1.30e-2 | 5.03e-3 | 1.93e-3 | 7.38e-4 | 2.81e-4 | 1.39 | $2-\alpha$
Table 2: Maximum errors as $t\rightarrow 0$ and convergence orders in temporal direction with $M=1000$ and $N=10$ for Example 1 $\alpha$ |  $\sigma$ $\backslash t_{N}$ | $1e-3$ | $1e-4$ | $1e-5$ | $1e-6$ | $1e-7$ | $Rate$ | expected order
---|---|---|---|---|---|---|---|---
0.4 | $0.1$ | 2.90e-2 | 2.39e-2 | 1.92e-2 | 1.54e-2 | 1.23e-2 | 0.09 | $\sigma$
| $0.4$ | 1.20e-3 | 5.02e-4 | 2.05e-4 | 8.22e-5 | 3.28e-5 | 0.40 | $\sigma$
| $0.6$ | 1.38e-4 | 3.67e-5 | 9.43e-6 | 2.39e-6 | 6.03e-7 | 0.59 | $\sigma$
| $1.2$ | 5.61e-7 | 3.70e-8 | 2.38e-9 | 1.51e-10 | 9.57e-12 | 1.19 | $\sigma$
| $1.8$ | 2.87e-8 | 4.71e-10 | 7.58e-12 | 1.21e-13 | 1.92e-15 | 1.80 | $\sigma$
0.6 | $0.4$ | 3.67e-3 | 1.49e-3 | 5.94e-4 | 2.37e-4 | 9.42e-5 | 0.40 | $\sigma$
| $0.6$ | 4.29e-4 | 1.09e-4 | 2.76e-5 | 6.95e-6 | 1.75e-6 | 0.60 | $\sigma$
| $0.8$ | 3.99e-5 | 6.43e-6 | 1.02e-6 | 1.62e-7 | 2.57e-8 | 0.80 | $\sigma$
| $1.2$ | 1.64e-6 | 1.05e-7 | 6.64e-9 | 4.49e-10 | 2.55e-11 | 1.20 | $\sigma$
| $1.8$ | 7.58e-7 | 1.21e-9 | 1.92e-11 | 3.05e-13 | 4.85e-14 | 1.80 | $\sigma$
Table 3: Maximum errors at $t=1$ and convergence orders in spatial direction with $N=1000$ for Example 1 $\alpha$ |  $\sigma$ $\backslash M$ | $8$ | $16$ | $24$ | $32$ | $40$ | $Rate$ | expected order
---|---|---|---|---|---|---|---|---
0.4 | $0.4$ | 9.25E-3 | 2.30e-3 | 1.02e-3 | 5.67e-4 | 3.60e-4 | 2.02 | $2$
| $1.2$ | 7.68e-3 | 1.92E-3 | 8.51e-4 | 4.79e-4 | 3.06e-4 | 2.00 | $2$
0.6 | $0.6$ | 8.20e-3 | 2.03e-3 | 8.96e-4 | 4.98e-4 | 3.14e-4 | 2.03 | $2$
| $0.2$ | 6.79e-3 | 1.70e-3 | 7.54e-4 | 4.24e-4 | 2.71e-4 | 2.00 | $2$

We investigate the convergence orders in temporal direction by setting
$M=1000$ and $t_{N}=1$. We show the maximum errors at $t_{N}=1$ in Table
LABEL:1d11. It can be seen from the results that when $t$ is far away from
$0$, the convergence order tends to $\sigma+1-\alpha$ when $0<\sigma<1$ and to
$2-\alpha$ when $1<\sigma<2$. We also show the maximum errors as $t\rightarrow
0$ in Table LABEL:1d12. The convergence results indicate that the convergence
order tends to $\sigma$. Then, we investigate the spatial convergence orders
by setting $N$=1000 and $t_{N}=1$. We show the numerical results in Table
LABEL:1d13. Clearly, the convergence order tends to $2$. All the numerical
results agree with the theoretical findings well.

Example 2. Consider the nonlinear subdiffusion problems (1.1) with the
following conditions

$\displaystyle(a)~{}d=1,f(u)=\sqrt{1+u^{2}}~{}\textrm{and}~{}u_{0}(x)=x(1-x),\Omega=(0,1),$
$\displaystyle(b)~{}d=1,f(u)=u-u^{3}~{}\textrm{and}~{}u_{0}(x)=\sin(\pi
x),\Omega=(0,1),$
$\displaystyle(c)~{}d=2,f(u)=\sqrt{1+u^{2}}~{}\textrm{and}~{}u_{0}(x)=\sin(\pi
x)\sin(\pi y),\Omega=(0,1)^{2},$
$\displaystyle(d)~{}d=2,f(u)=u-u^{3}~{}\textrm{and}~{}u_{0}(x)=x(1-x)y(1-y),\Omega=(0,\pi)^{2}.$

One can check that the initial condition $u_{0}\in H_{0}^{1}(\Omega)\cap
H^{2}(\Omega)$. Therefore, we have $\sigma=\alpha$. In the numerical
experiments, we take $M=1000$ for one-dimensional problem and $M=10$ for two
dimensional problem. The reference solutions are obtained by using small
temporal stepsizes. We list the maximum errors at $t=1$ and convergence orders
in Table LABEL:1d21. One can see that the convergence order tends to $1$. We
also show the maximum errors as $t\rightarrow 0$ and convergence orders in
Table LABEL:1d22. Clearly, the convergence order tends to $\alpha$. The
results further confirm the theoretical results of the present paper.

Table 4: Maximum errors at $t=1$ and convergence orders in temporal direction for Example 2 $\alpha$ | $case\backslash N$ | $10$ | $20$ | $40$ | $80$ | $160$ | $Rate$ | expected order
---|---|---|---|---|---|---|---|---
$0.4$ | $(a)$ | 1.74e-4 | 8.36e-5 | 4.04e-5 | 1.94e-5 | 8.99e-6 | 1.06 | 1
| $(b)$ | 1.49e-3 | 7.14e-4 | 3.46e-4 | 1.66e-4 | 7.69e-5 | 1.06 | 1
| $(c)$ | 9.63e-5 | 4.77e-5 | 2.36e-5 | 1.16e-5 | 5.61e-6 | 1.03 | 1
| $(d)$ | 7.20e-6 | 3.57e-6 | 1.76e-6 | 8.67e-7 | 4.19e-7 | 1.02 | 1
$0.6$ | $(a)$ | 2.15e-4 | 1.01e-4 | 4.84e-5 | 2.30e-5 | 1.06e-5 | 1.08 | 1
| $(b)$ | 1.88e-3 | 8.84e-4 | 4.23e-4 | 2.01e-4 | 9.27e-5 | 1.08 | 1
| $(c)$ | 1.06e-4 | 5.21e-5 | 2.57e-5 | 1.26e-5 | 6.09e-6 | 1.03 | 1
| $(d)$ | 7.94e-6 | 3.92e-6 | 1.93e-6 | 9.47e-7 | 4.57e-7 | 1.03 | 1
$0.8$ | $(a)$ | 2.13e-4 | 9.53e-5 | 4.44e-5 | 2.08e-5 | 9.48e-6 | 1.10 | 1
| $(b)$ | 1.96e-3 | 8.74e-4 | 4.07e-5 | 1.90e-5 | 8.66e-6 | 1.12 | 1
| $(c)$ | 7.91e-5 | 3.88e-5 | 1.91e-5 | 9.31e-6 | 4.47e-6 | 1.04 | 1
| $(d)$ | 6.02e-6 | 2.95e-6 | 1.45e-6 | 7.06e-7 | 3.40e-7 | 1.03 | 1
Table 5: Maximum errors as $t\rightarrow 0$ and convergence orders in temporal direction for Example 2 $\alpha$ | $case\backslash t_{N}$ | $1e-4$ | $1e-5$ | $1e-6$ | $1e-7$ | $1e-8$ | $Rate$ | expected order
---|---|---|---|---|---|---|---|---
$0.4$ | $(a)$ | 4.70e-4 | 2.21e-5 | 9.00e-5 | 3.61e-5 | 1.45e-5 | 0.38 | $\alpha$
| $(b)$ | 3.65e-3 | 1.87e-3 | 8.36e-4 | 3.50e-4 | 1.43e-4 | 0.35 | $\alpha$
| $(c)$ | 4.96e-3 | 3.05e-3 | 1.47e-3 | 6.34e-4 | 2.61e-4 | 0.32 | $\alpha$
| $(d)$ | 3.31e-4 | 1.87e-4 | 8.26e-5 | 3.39e-5 | 1.37e-5 | 0.35 | $\alpha$
$0.6$ | $(a)$ | 1.20e-4 | 3.05e-5 | 7.70e-6 | 1.94e-6 | 4.87e-7 | 0.60 | $\alpha$
| $(b)$ | 1.14e-3 | 2.99e-4 | 7.59e-5 | 1.91e-5 | 4.81e-6 | 0.59 | $\alpha$
| $(c)$ | 2.03e-3 | 5.46e-4 | 1.39e-4 | 3.52e-5 | 8.84e-6 | 0.59 | $\alpha$
| $(d)$ | 1.11e-4 | 2.86e-5 | 7.23e-6 | 1.82e-6 | 4.57e-7 | 0.60 | $\alpha$
$0.8$ | $(a)$ | 1.73e-5 | 2.75e-6 | 4.35e-7 | 6.86e-8 | 1.09e-8 | 0.80 | $\alpha$
| $(b)$ | 1.63e-4 | 2.59e-5 | 4.10e-6 | 6.51e-7 | 1.03e-8 | 0.80 | $\alpha$
| $(c)$ | 3.06e-4 | 4.77e-5 | 7.56e-6 | 1.20e-6 | 1.90e-7 | 0.80 | $\alpha$
| $(d)$ | 1.55e-5 | 2.46e-6 | 3.90e-7 | 6.18e-8 | 9.80e-9 | 0.80 | $\alpha$

##  5 Conclusions

In this paper, we consider the numerical solutions of the fully discrete
scheme for solving the nonlinear subdiffusion problems. The scheme is
constructed by using the L1-scheme in temporal directions and central finite
difference method for the diffusion term. The pointwise-in-time error
estimates of the L1 sheme are obtained under the regularity parameter
$\sigma\in(0,1)\cup(1,2)$ by carefully using the refined discrete fractional-
type Grönwall inequality. The present error estimates are sharp and
illustrated by several numerical experiments.

## References

  * [1] R. Gorenflo, F. Mainardi, D. Moretti, and P. Paradisi, _Time fractional diffusion: A discrete random walk approach_ , Nonlinear Dynam., 29 (2002), pp. 129–143.
  * [2] B. I. Henry and S. L. Wearne, _Fractional reaction-diffusion_ , Phys. A, 276 (2000), pp. 448–455.
  * [3] H. Liao, T. Tang, T. Zhou, _A second-order and nonuniform time-stepping maximum-principle preserving scheme for time-fractional Allen-Cahn equations_ , J. Comput. Phys., 414 (2020) 109473.
  * [4] B. Jin, B. Li, Z. Zhou, _Numerical analysis of nonlinear subdiffusion equations_ , SIAM J. Numer. Anal. 56(1) (2018) 1-23.
  * [5] M. Maskari, S. Karaa, _Numerical approximation of semilinear subdiffusion equations with nonsmooth intial data_ , SIAM J. Numer. Anal. 57 (2019) 1524–1544.
  * [6] N. Kopteva,_Error analysis of the L1 method on graded and uniform meshes for a fractional-derivative problem in two and three dimensions_ , Math. Comp., 88(2019), 2135–2155.
  * [7] W. Cao, F. Zeng, Z. Zhang, G.E.Karniadakis,_Implicit $-$explicit difference schemes for nonlinear fractional differential equations with nonsmooth solutions_, SIAM J. Sci. Comput., 38 (2016), A3070–A3093.
  * [8] B. Jin, R. Lazarov, Z. Zhou, _Two schemes for fractional diffusion and diffusion-wave equations_ , SIAM. J. Sci. Comput., 38 (2016), 146–170.
  * [9] D. Li, J. Wang, J. Zhang,_Unconditionally convergent $L1$-Galerkin FEMs for nonlinear time-fractional Schrödinger equations,_ SIAM. J. Sci. Comput., 39 (2017), A3067–A3088.
  * [10] D. Li, C. Wu, Z. Zhang. _Linearized Galerkin FEMs for nonlinear time fractional parabolic problems with nonsmooth solutions in time direction_ , J. Sci. Comput., 80 (2019) 403–419.
  * [11] W. McLean, K. Mustapha, R. Ali, O. M. Knio, _Well-posedness of time-fractional advection-diffusion-reaction equations._ Fract. Calc. Appl. Anal. 22 (2019), no. 4, 918–944.
  * [12] W. McLean, K. Mustapha, R. Ali, O. M. Knio, _Regularity theory for time-fractional advection-diffusion-reaction equations._ Comput. Math. Appl. 79 (2020), no. 4, 947–961.
  * [13] H. Liao, D. Li, J. Zhang, _Sharp error estimate of the nonuniform L1 formula for linear reaction-subdiffusion equations,_ SIAM J. Numer. Anal. 56 (2), (2018) 1112–1133.
  * [14] M. Stynes, E. O’Riordan, 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.
  * [15] H. Liao, W. Mclean, J. Zhang, _A discrete Grönwall inequality with application to numerical schemes for subdiffusion problems_ , SIAM J. Numer. Anal. 57 (2019) 218–237.
  * [16] N. Kopteva, _Error analysis of an L2–type method on graded meshes for a fractional–order parabolic problem_. Math. Comp. 90 (2021) 19–40.
  * [17] C. Lubich, _Convolution quadrature and discretized operational calculus I_ , Numer. Math. 52 (1988), pp. 129–145.
  * [18] B. Jin, R. Lazarov, Z. Zhou, _Numerical methods for time-fractional evolution equations with nonsmooth data: a concise overview_ , Comput. Methods Appl. Mech. Engrg. 346 (2019) 332–358.
  * [19] J. Gracia, E. ÓRiordan, M. Stynes, _Convergence in positive time for a finite difference method applied to a fractional convection-diffusion problem_ , Comput. Methods Appl. Math. 18 (2018) 33–42.
  * [20] Y. Yan, M. Khan, N.J. Ford, _An analysis of the modified L1 scheme for time-fractional partial differential equations with nonsmooth data_ , SIAM J. Numer. Anal. 56 (1) (2018) 210–227.
  * [21] B. Jin, R. Lazarov, Z. Zhou, _An analysis of the L1 scheme for the subdiffusion equation with nonsmooth data_ , IMA J. Numer. Anal. 36 (1) (2016) 197–221.
  * [22] B. Jin, B. Li, Z. Zhou, _Correction of high-order BDF convolution quadrature for fractional evolution equations_ , SIAM J. Sci. Comput. 39 (6) (2017) A3129–A3152.
  * [23] K. Mustapha, _Time-stepping discontinuous Galerkin methods for fractional diffusion problems_ , Numer. Math. 130 (3) (2015) 497–516.
  * [24] W. McLean, K. Mustapha, _Time-stepping error bounds for fractional diffusion problems with non-smooth initial data_ , J. Comput. Phys. 293 (2015) 201–217.
  * [25] D. Li, H. Liao, J. Wang, W. Sun, J. Zhang, _Analysis of L1-Galerkin FEMs for time-fractional nonlinear parabolic problems_. Commun. Comput. Phys., 24 (1) (2018), 86–103.