remarkRemark hypothesisHypothesis claimClaim Nonuniform approximations of
semilinear diffusion-wave equationsP. Lyu and S. Vong ex_supplement

# A symmetric fractional-order reduction method for direct nonuniform
approximations of semilinear diffusion-wave equations††thanks: This work was
partially supported by the Fundamental Research Funds for the Central
Universities (JBK2102010), the National Natural Science Foundation of China
(12071373), The Science and Technology Development Fund, Macau SAR (File no.
0005/2019/A) and the grant MYRG2018-00047-FST from University of Macau.

Pin Lyu School of Economic Mathematics, Southwestern University of Finance and
Economics, Chengdu, China<EMAIL_ADDRESS>Seakweng Vong Corresponding
author. Department of Mathematics, University of Macau, Macao, China. ().
<EMAIL_ADDRESS>
###### Abstract

We introduce a symmetric fractional-order reduction (SFOR) method to construct
numerical algorithms on general nonuniform temporal meshes for semilinear
fractional diffusion-wave equations. By using the novel order reduction
method, the governing problem is transformed to an equivalent coupled system,
where the explicit orders of time-fractional derivatives involved are all
$\alpha/2$ $(1<\alpha<2)$. The linearized L1 scheme and Alikhanov scheme are
then proposed on general time meshes. Under some reasonable regularity
assumptions and weak restrictions on meshes, the optimal convergence is
derived for the two kinds of difference schemes by $H^{2}$ energy method. An
adaptive time stepping strategy which based on the (fast linearized) L1 and
Alikhanov algorithms is designed for the semilinear diffusion-wave equations.
Numerical examples are provided to confirm the accuracy and efficiency of
proposed algorithms.

###### keywords:

diffusion-wave equation, weak singularity, nonuniform mesh, adaptive mesh

65M06, 65M12, 35B65, 35R11

##  1 Introduction

In this paper, we consider numerical methods of the semilinear diffusion-wave
equation:

(1) $\displaystyle{\cal D}_{t}^{\alpha}u=\nu^{2}\Delta u+f(u,{\bf
x},t),\quad{\bf x}\in\Omega,~{}t\in(0,T],$

subject to the initial conditions $u({\bf x},0)=\varphi({\bf x})$ and
$u_{t}({\bf x},0)=\tilde{\varphi}({\bf x})$ for ${\bf x}\in\Omega$, and the
homogeneous boundary condition $u({\bf x},t)=0$ for ${\bf
x}\in\partial\Omega$; where $\Omega=(x_{l},x_{r})\times(y_{l},y_{r})$,
$1<\alpha<2$, $\nu$ is a constant, and ${\cal D}_{t}^{\delta}$ denotes the
Caputo derivative of order $\delta$:

${\cal D}_{t}^{\delta}u(t):=({\cal
I}^{n-\delta}u^{(n)})(t)\quad\mbox{for}~{}t>0~{}\mbox{and}~{}n-1<\delta<n,$

in which ${\cal I}^{\beta}$ represents the Riemann-Liouville fractional
integral of order $\beta$:

${\cal
I}^{\beta}u(t):=\int_{0}^{t}\omega_{\beta}(t-s)u(s)\,\mathrm{d}s\quad\mbox{with}\quad\omega_{\beta}(t)=\frac{t^{\beta-1}}{\Gamma(\beta)}.$

The diffusion-wave equation, which is also called the time-fractional wave
equation, can be applied to describe evolution processes intermediate between
diffusion and wave propagation. For example, it governs the propagation of
mechanical waves in viscoelastic media [24, 25]. The practical applications of
equation (1) span diversely many disciplines, such as the image processing
[39, 3], the universal electromagnetic, acoustic and mechanical response [31].

It is well known that the solutions of the sub-diffusion equations (also
called the time-fractional diffusion equations) typically exhibit weak initial
singularities [9, 35, 34], and it causes that the traditional time-stepping
methods fail to preserve their desired convergence rate [9]. The same
phenomenon occurs for the diffusion-wave equations. For example, Jin, Lazarov
and Zhou [10, Theorem A.4] show that the solution of the linear diffusion-wave
equation ($f=f({\bf x},t)$) satisfies that
$\|\partial_{t}^{m}u\|_{L^{2}(\Omega)}\leq
C_{T}t^{\alpha-m}\|f\|_{W^{m-1,\infty}(0,T;L^{2}(\Omega))}$, $m=1,2$, if $f\in
W^{1,\infty}(0,T;L^{2}(\Omega))$ and $\varphi={\tilde{\varphi}}=0$. Other
studies on regularities can be found in [10, 26, 34]. Recently some excellent
works have been done on the numerical approximation of linear diffusion-wave
equations taking the weak initial singularities into account. The convolution
quadrature methods generated by backward difference formulas are rigorously
discussed in [10], where the first- and second-order temporal convergence
rates are obtained under proper assumptions of the given data, and their
discrete maximal regularities are further studied by Jin, Li and Zhou [11].
Lately, for the problem with nonsmooth data, a Petrov-Galerkin method and a
time-stepping discontinuous Galerkin method are proposed in [22] (Luo, Li and
Xie) and [14] (Li, Wang and Xie), where the temporal convergence rate is
$(3-\alpha)/2$-order and about first-order respectively. Numerical schemes
with classical L1 approximation in time and the standard P1-element in space
are also implemented in [13] to have the temporal accuracies of ${\cal
O}(\tau^{3-\alpha})$ and ${\cal O}(\tau^{2})$ provided the ratio
$\tau^{\alpha}/h^{2}_{\min}$ is uniformly bounded. We note that the numerical
methods in the above works [10, 11, 22, 14, 13] are implemented on uniform
temporal steps. On the other hand, Mustapha & McLean [29] and Mustapha &
Schötzau [30] considered the time-stepping discontinuous Galerkin methods on
nonuniform temporal meshes to solve the following kind of fractional wave
equation:

(2) $\displaystyle u_{t}+{\cal
I}^{\beta}Au(t)=f(t),\quad\mbox{for}\quad\beta\in(0,1)\quad\mbox{and}\quad
t\in(0,T],$

where $A$ is a self-adjoint linear elliptic spatial operator. It can be
observed that the above integro-differential problem is (mathematically)
equivalent to the linear case of (1) under suitable assumptions on $f$ and
initial data. Their methods are illuminating and efficient with good temporal
accuracies. Laplace transform methods and convolution quadrature methods on
uniform temporal steps are also discussed respectively by McLean &
Thom$\acute{\mbox{e}}$e [27, 28] and Cuesta et al. [4, 5, 6] for the above
integro-differential problem, where the reference [4] is for the semilinear
case $f(t)=f(u,\nabla u,x,t)$. However, the above numerical methods for
solving (2) may not be easily extended to the semilinear problem (1) due to
the nonlinearity ${\cal D}_{t}^{\beta}f(u,t)$.

To the best of our knowledge, there are still challenges for numerical methods
of the diffusion-wave equation. In this paper, we will address the following
issues: (i) establishing and analyzing difference schemes by the classical L1
[32] and Alikhanov [1] approximations on nonuniform temporal meshes
(especially on more general meshes) for the semilinear diffusion-wave equation
with typical weak singular solutions; (ii) studying efficient numerical
algorithms, such as the adaptive time-stepping algorithm, for the semilinear
diffusion-wave equation in order to deal with the highly oscillatory
variations in time since the problem (1) leads to a mixed behavior of
diffusion and wave propagation.

Before introducing our main approach, we review two classical and popular
algorithms. The first one is the L1 algorithm [32], which was generated by
Lagrange linear interpolation formula, it is a direct and convenient
approximation formula in constructing numerical methods for sub-diffusion
problems, e.g., [36, 21, 38] where it was employed on uniform temporal grids.
Recently, the L1 method on graded temporal meshes, with monotonically
increasing step sizes, was analyzed in Stynes, O’Riordan & Gracia [35] and
Kopteva [12] to resolve the sub-diffusion equations with weakly singular
solutions. The other one is the Alikhanov algorithm, which was firstly
proposed by Alikhanov [1] by combining linear and quadratic interpolations
skillfully at an off-set time point on uniform mesh for the sub-diffusion
problem with sufficiently smooth solution. Implementation of this algorithm on
graded mesh was discussed by Chen and Stynes [2] and the second-order
convergence concerning with the weak initial singularities was established.
Particularly, Liao, Li and Zhang [15] presented a novel and technical
framework to derive the optimal convergence result of the nonuniform L1
scheme. The techniques were then generalized in Liao, McLean and Zhang [17],
which were further extended to a linearized scheme for the semilinear sub-
diffusion equations [20] and the Alikhanov scheme on more general nonuniform
meshes [16]. We remark that the above methods [35, 12, 2, 15, 17, 20, 16] on
nonuniform meshes are all for sub-diffusion problems.

In view of the high efficiency and broad potential applications of the L1 and
the Alikhanov algorithms, it is of high scientific value to consider their
nonuniform versions for resolving (at least) the weak initial singularities of
the diffusion-wave problem. In [36], Sun and Wu utilized a fixed-order
reduction method, i.e., by taking an auxiliary function $v=u_{t}$, to rewrite
a linear case of equation (1) to the following coupled equations:

(3) $\displaystyle{\cal D}_{t}^{\alpha-1}v=\nu^{2}\Delta u+f({\bf x},t),$ (4)
$\displaystyle v=u_{t},$

for ${\bf x}\in\Omega,~{}t\in(0,T]$. We note that the time-fractional
derivative on the auxiliary function $v$ in equation (3) is of order
$\alpha-1$ which belongs to $(0,1)$, so the system (3)–(4) is not structure
consistency in the time derivative order point of view. The diffusion-wave
equation can be solved following the standard framework of the L1 method on
uniform temporal meshes. Although it is easy to extend the above order
reduction method to the corresponding nonuniform L2 scheme (see [36] for its
uniform version), we find that it may be difficult to establish its stability
and convergence on more general time meshes. Therefore, for the first time, we
present a new order reduction method by introducing a novel auxiliary function

${\bf v}={\cal D}_{t}^{\frac{\alpha}{2}}{\bf u},$

where ${\bf u}=u-t{\tilde{\varphi}}$, which is a non-fixed-order reduction
technique, and we call the symmetric fractional-order reduction (SFOR) method.
The semilinear diffusion-wave equation (1) is then skillfully rewritten as
coupled equations with nice structure or having the feature of structure
consistency, i.e. (8)–(9), see Section 2 for more details. Basing on this
equivalent formulation, we can construct the implicit and linearized L1 and
Alikhanov algorithms on possible nonuniform time partitions
$0=t_{0}<t_{1}<\cdots<t_{N}=T$ for a given positive integer $N$, and discuss
their unconditional convergence by utilizing the framework of [15, 17, 16,
20]. Throughout this paper we assume that the solution satisfies the following
regularity:

(5) $\displaystyle\|\partial_{t}^{(k)}u\|_{H^{4}(\Omega)}\leq
C_{u}(1+t^{\sigma_{1}-k})\quad\mbox{and}\quad\|\partial_{t}^{(k)}{\bf
v}\|_{H^{4}(\Omega)}\leq C_{u}(1+t^{\sigma_{2}-k}),\quad k=0,1,2,3,$

for $t\in(0,T]$, where $\sigma_{1}\in(1,2)\cup(2,3)$ and
$\sigma_{2}\in(\alpha/2,1)\cup(1,2)$ are two regularity parameters. Our
analysis are under the weak mesh assumption:

  * MA.

There is a constant $C_{\gamma}>0$ such that $\tau_{k}\leq
C_{\gamma}\tau\min\\{1,t_{k}^{1-1/\gamma}\\}$ for $1\leq k\leq N$, with
$t_{k}\leq C_{\gamma}t_{k-1}$ and $\tau_{k}/t_{k}\leq
C_{\gamma}\tau_{k-1}/t_{k-1}$ for $2\leq k\leq N$,

where $\gamma\geq 1$ is the mesh parameter, $\tau_{k}:=t_{k}-t_{k-1}$ denotes
the $k$-th time step size for $1\leq k\leq N$ and $\tau:=\max_{1\leq k\leq
N}\\{\tau_{k}\\}$.

We prove that our method can achieve the desired optimal temporal convergence
orders (see Theorem 3.6), that is ${\cal
O}(\tau^{\min\\{2-\frac{\alpha}{2},\gamma\sigma_{1},\gamma\sigma_{2}\\}})$ for
the nonuniform L1 algorithm and ${\cal
O}(\tau^{\min\\{2,\gamma\sigma_{1},\gamma\sigma_{2}\\}})$ for the nonuniform
Alikhanov algorithm. We note that the sum-of-exponentials approximations [8,
19] can also be directly adopted to the proposed nonuniform L1 and Alikhanov
algorithms to reduce the memory storage and computational costs. We further
design an adaptive time-stepping strategy according to the two kinds of
algorithms, which is robust and accurate for dealing with not only the weak
initial singularities but also the rapid temporal oscillations of the
semilinear diffusion-wave problem.

The main contributions of this paper are summarized below:

  * •

We propose a novel order reduction method (SFOR) which enables the nonuniform
L1 and Alikhanov algorithms for the semilinear diffusion-wave equation can be
constructed and analyzed.

  * •

Based on some reasonable regularity assumptions and weak mesh restrictions, we
obtain the optimal convergence orders: the temporal convergence rate is up to
$(2-\alpha/2)$-order for the L1 algorithm and second-order for the Alikhanov
algorithm.

  * •

An adaptive time-stepping strategy is designed for the semilinear diffusion-
wave equation to efficiently resolve possible oscillations of the solution.

The rest of the paper is organized as follows. In section 2, we present the
novel SFOR method and equivalently rewrite the semilinear diffusion-wave
equation into coupled equations. In Section 3, we construct and analyze the
linearized nonuniform L1 and the nonuniform Alikhanov algorithms, and obtain
their optimal convergences unconditionally by $H^{2}$ energy method.
Furthermore, we design an adaptive time-stepping method by combining the
proposed nonuniform (fast linearized) L1 and Alikhanov algorithms. Numerical
examples are provided in Section 4 to demonstrate the accuracy and efficiency.
A brief conclusion is followed in Section 5, and the analysis of truncation
errors is given in Section 7.

Throughout the paper, we use $C$ to denote a generic constant which may
depends on the data of the governing problem but is independent of time and
space step sizes (or nodes).

##  2 The SFOR method

In this section, we propose a symmetric fractional-order reduction (SFOR)
method such that the technical framework proposed in [15, 17, 16, 20] can be
adopted to analyze implicit numerical schemes for solving the diffusion-wave
equation (1) on temporal nonuniform mesh.

The basic idea of the SFOR method is presented in the following lemma.

######  Lemma 2.1.

For $\alpha\in(1,2)$ and $u(t)\in{\cal C}^{1}([0,T])\cap{\cal C}^{2}((0,T])$,
it holds that

(6) $\displaystyle{\cal D}_{t}^{\alpha}u(t)={\cal
D}_{t}^{\frac{\alpha}{2}}\left({\cal
D}_{t}^{\frac{\alpha}{2}}u(t)\right)-u^{\prime}(0)\omega_{2-\alpha}(t).$

Moreover, if we take ${\bf u}(t):=u(t)-tu^{\prime}(0)$, then

(7) $\displaystyle{\cal D}_{t}^{\alpha}{u}(t)={\cal D}_{t}^{\alpha}{\bf
u}(t)={\cal D}_{t}^{\frac{\alpha}{2}}\left({\cal D}_{t}^{\frac{\alpha}{2}}{\bf
u}(t)\right).$

######  Proof 2.2.

Taking $v(t):={\cal D}_{t}^{\frac{\alpha}{2}}u(t)$, one has

$\displaystyle v={\cal D}_{t}^{\frac{\alpha}{2}}u(t)=({\cal
I}^{\frac{\alpha}{2}}u^{\prime})(t)=$
$\displaystyle\int_{0}^{t}\omega_{1-\frac{\alpha}{2}}(t-s)u^{\prime}(s)\,\mathrm{d}s$
$\displaystyle=$
$\displaystyle-u^{\prime}(s)\omega_{2-\frac{\alpha}{2}}(t-s)|_{0}^{t}+\int_{0}^{t}\omega_{2-\frac{\alpha}{2}}(t-s)u^{\prime\prime}(s)\,\mathrm{d}s$
$\displaystyle=$ $\displaystyle
u^{\prime}(0)\omega_{2-\frac{\alpha}{2}}(t)+\int_{0}^{t}\omega_{2-\frac{\alpha}{2}}(t-s)u^{\prime\prime}(s)\,\mathrm{d}s,$

where the integration by parts has been utilized. Then

$\displaystyle v_{t}=\frac{\,\mathrm{d}}{\,\mathrm{d}t}{\cal
D}_{t}^{\frac{\alpha}{2}}u(t)=u^{\prime}(0)\omega_{1-\frac{\alpha}{2}}(t)+\int_{0}^{t}\omega_{1-\frac{\alpha}{2}}(t-s)u^{\prime\prime}(s)\,\mathrm{d}s=u^{\prime}(0)\omega_{1-\frac{\alpha}{2}}(t)+({\cal
I}^{1-\frac{\alpha}{2}}u^{\prime\prime})(t).$

Hence, using the composition property ${\cal I}^{p}{\cal I}^{q}g(t)={\cal
I}^{p+q}g(t)~{}(p,q>0)$ [33, pp. 59], we get

$\displaystyle{\cal D}_{t}^{\frac{\alpha}{2}}\left({\cal
D}_{t}^{\frac{\alpha}{2}}u(t)\right)={\cal
D}_{t}^{\frac{\alpha}{2}}v(t)=({\cal I}^{1-\frac{\alpha}{2}}v^{\prime})(t)=$
$\displaystyle u^{\prime}(0)({\cal
I}^{1-\frac{\alpha}{2}}\omega_{1-\frac{\alpha}{2}})(t)+{\cal
I}^{1-\frac{\alpha}{2}}{\cal I}^{1-\frac{\alpha}{2}}u^{\prime\prime}(t)$
$\displaystyle=$ $\displaystyle u^{\prime}(0)\omega_{2-\alpha}(t)+{\cal
I}^{2-\alpha}u^{\prime\prime}(t)$ $\displaystyle=$ $\displaystyle
u^{\prime}(0)\omega_{2-\alpha}(t)+{\cal D}_{t}^{\alpha}u(t),$

implying (6) is true.

The equality (7) can be obtained directly by taking ${\bf v}:={\cal
D}_{t}^{\frac{\alpha}{2}}{\bf u}(t)$ in the above derivations.

Now we take

${\bf u}({\bf x},t):=u({\bf x},t)-t{\tilde{\varphi}}({\bf
x})\quad\mbox{and}\quad{\bf v}({\bf x},t):={\cal D}_{t}^{\frac{\alpha}{2}}{\bf
u}({\bf x},t).$

From (5) and using the Sobolev embedding theorem, we have $\|{\bf u}_{t}({\bf
x},t)\|_{\infty}\leq C(1+t^{\sigma_{1}-1})$ for $t\in(0,T]$. Then utilizing
the Comparison theorem for integrals (see pp. 400–401 in [40]), one has

$\displaystyle|{\bf v}({\bf x},0)|=$ $\displaystyle\left|\lim_{t\rightarrow
0}{\cal D}_{t}^{\frac{\alpha}{2}}{\bf u}({\bf
x},t)\right|\leq\frac{1}{\Gamma(1-\frac{\alpha}{2})}\lim_{t\rightarrow
0}\int_{0}^{t}(t-s)^{-\frac{\alpha}{2}}\left|{\bf u}_{t}({\bf
x},s)\right|\,\mathrm{d}s$ $\displaystyle\leq$ $\displaystyle
C\lim_{t\rightarrow
0}\int_{0}^{t}(t-s)^{-\frac{\alpha}{2}}(1+s^{\sigma_{1}-1})\,\mathrm{d}s\leq
C\lim_{t\rightarrow
0}(t^{1-\frac{\alpha}{2}}+t^{\sigma_{1}-\frac{\alpha}{2}})=0,$

which gives ${\bf v}({\bf x},0)=0.$

Thus, by Lemma 2.1, the equation (1) can be equivalently solved by the
following coupled equations:

(8) $\displaystyle{\cal D}_{t}^{\frac{\alpha}{2}}{\bf v}=\nu^{2}\Delta{\bf
u}+{f}(u,{\bf x},t)+t\Delta{\tilde{\varphi}},\quad{\bf
x}\in\Omega,~{}t\in(0,T],$ (9) $\displaystyle{\bf v}={\cal
D}_{t}^{\frac{\alpha}{2}}{\bf u},\quad{\bf x}\in\Omega,~{}t\in(0,T],$

provided $u={\bf u}+t{\tilde{\varphi}}$, the initial conditions ${\bf u}({\bf
x},0)=\varphi({\bf x}),~{}{\bf v}({\bf x},0)=0$ for ${\bf x}\in\Omega$, and
boundary conditions ${\bf u}({\bf x},t)={\bf v}({\bf x},t)=0$ for ${\bf
x}\in\partial\Omega$.

One can observe that, by utilizing the proposed SFOR method, the explicit
orders of the time-fractional derivatives in the resulting coupled equations
(8) and (9) are all $\alpha/2$. Therefore, they can be discretized by the same
strategy (e.g., the L1 or Alikhanov approximations).

######  Remark 2.3.

We observe from our numerical experiments that, by extracting the singular
term $u^{\prime}(0)\omega_{2-\alpha}(t)$ in (6), the proposed algorithms will
have more regular accuracy due to the regularity of the remaining part. This
is the reason why we define the auxiliary function ${\bf v}={\cal
D}_{t}^{\frac{\alpha}{2}}{\bf u}$ with ${\bf u}=u-t{\tilde{\varphi}}$, instead
of ${v}={\cal D}_{t}^{\frac{\alpha}{2}}{u}$.

##  3 Numerical algorithms

###  3.1 Preliminary

Our main concern is the time approximation of (1). Here and hereafter, $g^{k}$
and $g_{h}^{k}$ denotes the numerical approximations of $g(t_{k})$ and $g({\bf
x}_{h},t_{k})$, respectively. Define the off-set time points and grid
functions

$t_{n-\theta}:=\theta t_{n-1}+(1-\theta)t_{n}\quad\mbox{and}\quad
g^{n-\theta}:=\theta g^{n-1}+(1-\theta)g^{n},\quad 1\leq n\leq N.$

Denote $\beta:=\alpha/2$. The Caputo derivative ${\cal
D}_{t}^{\beta}g(t_{n-\theta})$ can be formally approximated by the following
discrete Caputo derivative with convolution structure:

(10) $\displaystyle({\cal
D}_{\tau}^{\beta}g)^{n-\theta}:=\sum_{k=1}^{n}A_{n-k}^{(n)}\nabla_{\tau}g^{k},\quad\mbox{where}~{}\nabla_{\tau}g^{k}=g^{k}-g^{k-1}.$

The general discretization (10) includes two practical ones. It leads to the
L1 formula while $\theta=0$ (see also (11)) and yields the Alikhanov formula
while $\theta=\beta/2$ (see also (13)). To efficiently solve the semilinear
diffusion-wave equation with possible weak singular or more complicated
solutions, we next give more explicit formulations of these two classical
approximations on possible nonuniform meshes, which have also been rigorously
studied in [15, 17, 16].

Nonuniform L1 formula. The L1 formula on general mesh for the approximation of
the Caputo derivative ${\cal D}_{t}^{\beta}g(t_{n})$ is given as:

(11) $\displaystyle({\cal
D}_{\tau}^{\beta}g)^{n}:=\sum_{k=1}^{n}\int_{t_{k-1}}^{t_{k}}\omega_{1-\beta}(t_{n}-s)(\Pi_{1,k}g(s))^{\prime}\,\mathrm{d}s=\sum_{k=1}^{n}A_{n-k}^{(n)}\nabla_{\tau}g^{k},$

where $\Pi_{1,k}$ represents the linear interpolation operator, and

(12) $\displaystyle
A_{n-k}^{(n)}:=\int_{t_{k-1}}^{t_{k}}\frac{\omega_{1-\beta}(t_{n}-s)}{\tau_{k}}\,\mathrm{d}s.$

Nonuniform Alikhanov formula. Denote $\theta:=\beta/2=\alpha/4$, and define
the discrete coefficients

$\displaystyle
a_{n-k}^{(n)}:=\frac{1}{\tau_{k}}\int_{t_{k-1}}^{\min\\{t_{k},t_{n-\theta}\\}}\omega_{1-\beta}(t_{n-\theta}-s)\,\mathrm{d}s,~{}1\leq
k\leq n;$ $\displaystyle
b_{n-k}^{(n)}:=\frac{2}{\tau_{k}(\tau_{k}+\tau_{k+1})}\int_{t_{k-1}}^{t_{k}}\omega_{1-\beta}(t_{n-\theta}-s)(s-t_{k-\frac{1}{2}})\,\mathrm{d}s,~{}1\leq
k\leq n-1.$

Referring to [16], the Alikhanov formula on general mesh for the approximation
of the Caputo derivative ${\cal D}_{t}^{\beta}g(t_{n-\theta})$ is

$\displaystyle({\cal D}_{\tau}^{\beta}g)^{n-\theta}:=$
$\displaystyle\sum_{k=1}^{n-1}\int_{t_{k-1}}^{t_{k}}\omega_{1-\beta}(t_{n-\theta}-s)(\Pi_{2,k}g(s))^{\prime}\,\mathrm{d}s+\int_{t_{n-1}}^{t_{n-\theta}}\omega_{1-\beta}(t_{n-\theta}-s)(\Pi_{1,n}g(s))^{\prime}\,\mathrm{d}s$
(13) $\displaystyle=$
$\displaystyle\sum_{k=1}^{n}A_{n-k}^{(n)}\nabla_{\tau}g^{k},$

where $\Pi_{2,k}$ denotes the quadratic interpolation operator, and the
discrete convolution kernels $A_{n-k}^{(n)}$ here are given as follows:
$A_{0}^{(1)}:=a_{0}^{(1)}$ for $n=1$, and

(17) $\displaystyle
A_{n-k}^{(n)}:=\left\\{\begin{array}[]{ll}a_{0}^{(n)}+\rho_{n-1}b_{1}^{(n)},&k=n,\\\
a_{n-k}^{(n)}+\rho_{k-1}b_{n-k+1}^{(n)}-b_{n-k}^{(n)},&2\leq k\leq n-1,\\\
a_{n-1}^{(n)}-b_{n-1}^{(n)},&k=1,\end{array}\right.\quad\mbox{for}~{}n\geq 2,$

with $\rho_{k}:=\tau_{k}/\tau_{k+1}$ and $\rho:=\max_{k}\\{\rho_{k}\\}$ being
the local time step-size ratios and the maximum ratio, respectively.

######  Remark 3.1.

In the rest of this paper, we will use the general form (10) to represent the
nonuniform L1 formula and Alikhanov formula. The discrete coefficients
$A_{n-k}^{(n)}$ and the related properties studied later correspondingly refer
to those of the nonuniform L1 formula and the Alikhanov formula while
$\theta=0$ and $\theta=\beta/2$, respectively.

The following two basic properties have been verified in [17, 16] for the
discrete coefficients of the nonuniform L1 formula (with $\pi_{A}=1$) and the
nonuniform Alikhanov formula (with $\pi_{A}=11/4$ and $\rho=7/4$), which are
required in the numerical analysis of corresponding algorithms:

  * A1.

The discrete kernels are positive and monotone: $A_{0}^{(n)}\geq
A_{1}^{(n)}\geq\cdots\geq A_{n-1}^{(n)}>0$;

  * A2.

There is a constant $\pi_{A}>0$ such that
$A_{n-k}^{(n)}\geq\frac{1}{\pi_{A}}\int_{t_{k-1}}^{t_{k}}\frac{\omega_{1-\beta}(t_{n}-s)}{\tau_{k}}\,\mathrm{d}s$
for $1\leq k\leq n\leq N$.

With A1–A2, a natural and important property is valid for the nonuniform L1
formula [15, proof of Theorem 2.1] and the nonuniform Alikhanov formula [16,
Corollary 2.3]:

(18) $\displaystyle\left\langle({\cal
D}_{\tau}^{\beta}g)^{n-\theta},g^{n-\theta}\right\rangle\geq\frac{1}{2}\sum_{k=1}^{n}A_{n-k}^{(n)}\nabla_{\tau}(\|g^{k}\|^{2})\quad\mbox{for}~{}1\leq
n\leq N.$

A discrete fractional Grönwall inequality proposed in [17, Theorem 3.1] is a
crucial tool in the numerical analysis of fractional problems. As required in
the analysis later, we present a slightly modified version in the following.
It is easy to trace the proof of [17, Theorem 3.1] to justify the
modification, here we skip its trivial derivations.

######  Lemma 3.2.

Let $(g^{n})_{n=1}^{N}$ and $(\lambda_{l})_{l=0}^{N-1}$ be given nonnegative
sequences. Assume that there exists a constant $\Lambda$ (independent of the
step sizes) such that $\Lambda\geq\sum_{l=0}^{N-1}\lambda_{l}$, and that the
maximum step size satisfies

$\max_{1\leq n\leq
N}\tau_{n}\leq\frac{1}{{{}^{\beta}\sqrt{4\pi_{A}\Gamma(2-\beta)\Lambda}}}.$

Then, for any nonnegative sequence $(v^{k})_{k=0}^{N}$ and $(w^{k})_{k=0}^{N}$
satisfying

$\displaystyle\sum_{k=1}^{n}A_{n-k}^{(n)}\nabla_{\tau}\left[(v^{k})^{2}+(w^{k})^{2}\right]\leq\sum_{k=1}^{n}\lambda_{n-k}\left(v^{k-\theta}+w^{k-\theta}\right)^{2}+(v^{n-\theta}+w^{n-\theta})g^{n},~{}1\leq
n\leq N,$

it holds that

(19) $\displaystyle v^{n}+w^{n}\leq 4E_{\beta}(4\max(1,\rho)\pi_{A}\Lambda
t_{n}^{\beta})\left(v^{0}+w^{0}+\max_{1\leq k\leq
n}\sum_{j=1}^{k}P_{k-j}^{(k)}g^{j}\right)\quad\mbox{for}~{}1\leq n\leq N,$

where $E_{\beta}(z)=\sum_{k=0}^{\infty}\frac{z^{k}}{\Gamma(1+k\beta)}$ is the
Mittag-Leffler function.

The coefficients $P_{n-j}^{(n)}$ in (19) are called the complementary discrete
kernels ([17]) which are defined on the convolution coefficients
$A_{n-k}^{(n)}$ :

(20) $\displaystyle P_{0}^{(n)}:=\frac{1}{A_{0}^{(n)}},\quad
P_{n-j}^{(n)}:=\frac{1}{A_{0}^{(j)}}\sum_{k=j+1}^{n}\left(A_{k-j-1}^{(k)}-A_{k-j}^{(k)}\right)P_{n-k}^{(n)},\quad
1\leq j\leq n-1.$

It has been shown in [17, Lemmma 2.1] that the kernels satisfy

(21) $\displaystyle 0\leq
P_{n-j}^{(n)}\leq\pi_{A}\Gamma(2-\beta)\tau_{j}^{\beta},\quad\sum_{j=1}^{n}P_{n-j}^{(n)}\omega_{1-\beta}(t_{j})\leq\pi_{A},\quad
1\leq j\leq n\leq N.$

###  3.2 Nonuniform L1 and Alikhanov algorithms

We now implement linearized algorithms on temporal nonuniform meshes to solve
the coupled equations (8)–(9) based on the nonuniform L1 and Alikhanov
formulas.

Some basic notations in the spatial direction are needed. The uniform spatial
step sizes are denoted by $h_{x}:=(x_{r}-x_{l})/M_{x}$ and
$h_{y}:=(y_{r}-y_{l})/M_{y}$ respectively, where $M_{x},M_{y}$ are positive
integers. The mesh space is given by $\bar{\Omega}_{h}:=\\{{\bf
x}_{h}=(x_{l}+ih_{x},y_{l}+jh_{y})|0\leq i\leq M_{x},0\leq j\leq M_{y}\\}$.
For any grid functions
$u_{h}:=\\{u_{i,j}=u(x_{i},t_{j})|(x_{i},t_{j})\in\bar{\Omega}_{h}\\}$, we
employ standard five-point finite difference operator
$\Delta_{h}:=\delta_{x}^{2}+\delta_{y}^{2}$ on $\bar{\Omega}_{h}$ to
discretize the Laplacian operator $\Delta$, where
$\delta_{x}^{2}u_{i,j}:=(u_{i+1,j}-2u_{i,j}+u_{i-1,j})/{h_{x}^{2}}$ and
$\delta_{y}^{2}u_{i,j}$ is defined similarly.

Denote

$F(u_{h}^{n-\theta}):=f(u_{h}^{n-1},{\bf
x}_{h},t_{n-\theta})+(1-\theta)\partial_{u}f(u_{h}^{n-1},{\bf
x}_{h},t_{n-\theta})(u_{h}^{n}-u_{h}^{n-1}),\quad 1\leq n\leq N.$

The linearized and implicit difference schemes which based on the L1 and the
Alikhanov approximations on general nonuniform temporal meshes for the problem
(8)–(9) or the problem (1) are constructed as follows:

(22) $\displaystyle({\cal D}_{\tau}^{\beta}{\bf
v}_{h})^{n-\theta}=\nu^{2}\Delta_{h}{\bf
u}_{h}^{n-\theta}+F({u}_{h}^{n-\theta})+t_{n-\theta}\Delta{\tilde{\varphi}}_{h},\quad{\bf
x}\in\Omega_{h},~{}1\leq n\leq N;$ (23) $\displaystyle{\bf
v}_{h}^{n-\theta}=({\cal D}_{\tau}^{\beta}{\bf u}_{h})^{n-\theta},\quad{\bf
x}\in\Omega_{h},~{}1\leq n\leq N;$ (24) $\displaystyle u_{h}^{n}={\bf
u}_{h}^{n}+t_{n}{\tilde{\varphi}}_{h},\quad{\bf x}\in\Omega_{h},~{}0\leq n\leq
N,$

equipped with the initial conditions ${\bf u}_{h}^{0}=\varphi({\bf x}_{h})$
and ${\bf v}_{h}^{0}=0$ for ${\bf x}\in\Omega_{h}$, and the boundary
conditions ${\bf u}_{h}^{n}={\bf v}_{h}^{n}=0$ for ${\bf
x}\in\partial\Omega_{h},~{}1\leq n\leq N$.

######  Remark 3.3.

The equations (22)–(24) represent two different numerical algorithms for
solving the semilinear diffusion-wave equation. It is the nonuniform L1
algorithm while $\theta=0$ and is the nonuniform Alikhanov algorithm while
$\theta=\beta/2=\alpha/4$.

In order to analyze the two proposed algorithms, we consider an equivalent
form of (22)–(24). Firstly, denote $w:={\cal D}_{t}^{\beta}{\bf v}-f(u,{\bf
x},t)+t\Delta{\tilde{\varphi}}$ with the initial condition $w({\bf
x},0):=\nu^{2}\Delta\varphi$ and the boundary $w({\bf x},t):=-f(0,{\bf x},t)$.
Then (8)–(9) can be rewritten as

$\displaystyle w={\cal D}_{t}^{\beta}{\bf v}-f(u,{\bf
x},t)+t\Delta{\tilde{\varphi}},\quad{\bf x}\in\Omega,~{}t\in(0,T];$
$\displaystyle w=\nu^{2}\Delta{\bf u},\quad{\bf x}\in\Omega,~{}t\in(0,T];$
$\displaystyle{\bf v}={\cal D}_{t}^{\beta}{\bf u},\quad{\bf
x}\in\Omega,~{}t\in(0,T].$

Utilizing the nonuniform L1 and Alikhanov formulas and the linearized
technique to approximate the above equations, we obtain an auxiliary system of
(22)–(24):

(25) $\displaystyle w_{h}^{n-\theta}=({\cal D}_{\tau}^{\beta}{\bf
v}_{h})^{n-\theta}-F(u_{h}^{n-\theta})+t_{n-\theta}\Delta{\tilde{\varphi}}_{h},\quad{\bf
x}_{h}\in\Omega_{h},~{}1\leq n\leq N;$ (26) $\displaystyle
w_{h}^{n}=\nu^{2}\Delta_{h}{\bf u}_{h}^{n},\quad{\bf
x}_{h}\in\Omega_{h},~{}0\leq n\leq N;$ (27) $\displaystyle{\bf
v}_{h}^{n-\theta}=({\cal D}_{\tau}^{\beta}{\bf u}_{h})^{n-\theta},\quad{\bf
x}_{h}\in\Omega_{h},~{}1\leq n\leq N;$ (28) $\displaystyle u_{h}^{n}={\bf
u}_{h}^{n}+t_{n}{\tilde{\varphi}}_{h},\quad{\bf x}\in\Omega_{h},~{}0\leq n\leq
N.$

As $w_{h}^{n-\theta}=(1-\theta)w_{h}^{n}+\theta w_{h}^{n-1}$, one can see that
equations (22)–(24) are equivalent to (25)–(28) by eliminating the functions
$w_{h}^{n-\theta}$ and $w_{h}^{n}$.

###  3.3 Unconditional convergence

In this subsection, we show the unconditional convergence of the proposed
nonuniform L1 and Alikhanov algorithms (22)–(24) according to their auxiliary
system (25)–(28).

Take $\Omega_{h}=\bar{\Omega}_{h}\cap\Omega$ and
$\partial\Omega_{h}=\bar{\Omega}_{h}\cap\partial\Omega$. For $u_{h},v_{h}$
belonging to the space of grid functions which vanish on $\partial\Omega_{h}$,
we introduce the discrete inner product $\langle
u,v\rangle:=h_{x}h_{y}\sum_{{\bf x}_{h}\in\Omega_{h}}u_{h}v_{h}$, the discrete
$L_{2}$-norm $\|u\|:=\sqrt{\langle u,u\rangle}$, the discrete
$L_{\infty}$-norm $\|u\|_{\infty}:=\max\\{|u_{h}|\\}$, the discrete $H^{1}$
seminorms $\|\delta_{x}u\|$ and $\|\delta_{y}u\|$, and
$\|\nabla_{h}u\|:=\sqrt{\|\delta_{x}u\|^{2}+\|\delta_{y}u\|^{2}}$, where
$\delta_{x}u_{i-\frac{1}{2},j}:=(u_{i,j}-u_{i-1,j})/h_{x}$ and similar
definition works for $\delta_{y}u_{i,j-\frac{1}{2}}$. One can easily check
that $\langle\Delta_{h}u,u\rangle=-\|\nabla_{h}u\|^{2}$, and, for some
positive constants ${\tilde{C}}_{\Omega},{\hat{C}}_{\Omega}$, the embedding
inequalities are valid ([18]): $\|u\|\leq{\tilde{C}}_{\Omega}\|\nabla_{h}u\|$
and
$\max\\{\|\nabla_{h}u\|,\|u\|_{\infty}\\}\leq{\hat{C}}_{\Omega}\|\Delta_{h}u\|$.
For simplicity of presentation, we take $h:=\max\\{h_{x},h_{y}\\}$ and
$C_{\Omega}=\max\\{{\tilde{C}}_{\Omega},{\hat{C}}_{\Omega}\\}$.

For $\vartheta\in(0,1]$, let $U_{h}^{n-\vartheta}:=u({\bf
x}_{h},t_{n-\vartheta})$ and denote the solution errors

${\tilde{u}}_{h}^{n}:=U_{h}^{n}-u_{h}^{n}={\bf u}({\bf x}_{h},t_{n})-{\bf
u}_{h}^{n},\quad{\tilde{v}}_{h}^{n}:={\bf v}({\bf x}_{h},t_{n})-{\bf
v}_{h}^{n},\quad\mbox{and}\quad{\tilde{w}}_{h}^{n}:=w({\bf
x}_{h},t_{n})-w_{h}^{n}.$

One can obtain the error system of (25)–(28):

(29) $\displaystyle{\tilde{w}}_{h}^{n-\theta}=({\cal
D}_{\tau}^{\beta}{\tilde{v}}_{h})^{n-\theta}-{\cal N}_{h}^{n-\theta}-({\cal
T}_{f})_{h}^{n-\theta}+({\cal T}_{v1})_{h}^{n-\theta}-({\cal
T}_{w})_{h}^{n-\theta},\quad{\bf x}_{h}\in\Omega_{h},~{}1\leq n\leq N;$ (30)
$\displaystyle{\tilde{w}}_{h}^{n}=\nu^{2}\Delta_{h}{\tilde{u}}_{h}^{n}+\nu^{2}{\cal
S}_{h}^{n},\quad{\bf x}_{h}\in\Omega_{h},~{}1\leq n\leq N;$ (31)
$\displaystyle{\tilde{v}}_{h}^{n-\theta}=({\cal
D}_{\tau}^{\beta}{\tilde{u}}_{h})^{n-\theta}+({\cal
T}_{u})_{h}^{n-\theta}-({\cal T}_{v2})_{h}^{n-\theta},\quad{\bf
x}_{h}\in\Omega_{h},~{}1\leq n\leq N;$
$\displaystyle{\tilde{u}}_{h}^{0}={\tilde{v}}_{h}^{0}={\tilde{w}}_{h}^{0}=0,\quad{\bf
x}_{h}\in\bar{\Omega}_{h};\qquad{\tilde{u}}_{h}^{n}={\tilde{v}}_{h}^{n}=0,\quad{\bf
x}_{h}\in\partial\Omega_{h},~{}1\leq n\leq N.$

where $({\cal T}_{f})_{h}^{n-\theta}$, $({\cal T}_{v1})_{h}^{n-\theta}$,
$({\cal T}_{w})_{h}^{n-\theta}$, $({\cal T}_{u})_{h}^{n-\theta}$, $({\cal
T}_{v2})_{h}^{n-\theta}$ and ${\cal S}_{h}^{n}$ are the temporal and spatial
truncation errors, see more details in Appendix (Section 7); and

$\displaystyle{\cal N}_{h}^{n-\theta}:=$ $\displaystyle
F(U_{h}^{n-\theta})-F(u_{h}^{n-\theta})$ $\displaystyle=$
$\displaystyle(1-\theta)\left[\partial_{u}f(u_{h}^{n-1},{\bf
x}_{h},t_{n-\theta})\nabla_{\tau}{\tilde{u}}_{h}^{n}+{\tilde{u}}_{h}^{n-1}\nabla_{\tau}U_{h}^{n}\int_{0}^{1}\partial^{2}_{u}f(sU_{h}^{n-1}+(1-s)u_{h}^{n-1},{\bf
x}_{h},t_{n-\theta})\,\mathrm{d}s\right]$
$\displaystyle+{\tilde{u}}_{h}^{n-1}\int_{0}^{1}\partial_{u}f(sU_{h}^{n-1}+(1-s)u_{h}^{n-1},{\bf
x}_{h},t_{n-\theta})\,\mathrm{d}s.$

It can be deduced from (30) that
$\sum_{k=1}^{n}A_{n-k}^{(n)}\nabla_{\tau}{\tilde{w}}_{h}^{k}=\nu^{2}\sum_{k=1}^{n}A_{n-k}^{(n)}\nabla_{\tau}\left(\Delta_{h}{\tilde{u}}_{h}^{k}+{\cal
S}_{h}^{k}\right)$. Then, performing the operator $\Delta_{h}$ on (29) and
(31), it follows

$\displaystyle\Delta_{h}{\tilde{w}}_{h}^{n-\theta}=({\cal
D}_{\tau}^{\beta}\Delta_{h}{\tilde{v}}_{h})^{n-\theta}+\Delta_{h}{\cal
N}_{h}^{n-\theta}-\Delta_{h}({\cal T}_{f})_{h}^{n-\theta}+\Delta_{h}({\cal
T}_{v1})_{h}^{n-\theta}-\Delta_{h}({\cal T}_{w})_{h}^{n-\theta},$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}{\bf
x}_{h}\in\Omega_{h},~{}1\leq n\leq N;$ (32) $\displaystyle({\cal
D}_{\tau}^{\beta}{\tilde{w}}_{h})^{n-\theta}=\nu^{2}({\cal
D}_{\tau}^{\beta}\Delta_{h}{\tilde{u}}_{h})^{n-\theta}+\nu^{2}({\cal
D}_{\tau}^{\beta}{\cal S}_{h})^{n-\theta},\quad{\bf
x}_{h}\in\Omega_{h},~{}1\leq n\leq N;$ (33)
$\displaystyle\Delta_{h}{\tilde{v}}_{h}^{n-\theta}=({\cal
D}_{\tau}^{\beta}\Delta_{h}{\tilde{u}}_{h})^{n-\theta}+\Delta_{h}({\cal
T}_{u})_{h}^{n-\theta}-\Delta_{h}({\cal T}_{v2})_{h}^{n-\theta},\quad{\bf
x}_{h}\in\Omega_{h},~{}1\leq n\leq N;$
$\displaystyle{\tilde{u}}_{h}^{0}={\tilde{v}}_{h}^{0}={\tilde{w}}_{h}^{0}=0,\quad{\bf
x}_{h}\in\bar{\Omega}_{h};\qquad{\tilde{u}}_{h}^{n}={\tilde{v}}_{h}^{n}=0,\quad{\bf
x}_{h}\in\partial\Omega_{h},~{}1\leq n\leq N.$

By eliminating the term $({\cal
D}_{\tau}^{\beta}\Delta_{h}{\tilde{u}}_{h})^{n-\theta}$ in (32) and (33), we
get

$\displaystyle\Delta_{h}{\tilde{w}}_{h}^{n-\theta}=({\cal
D}_{\tau}^{\beta}\Delta_{h}{\tilde{v}}_{h})^{n-\theta}+\Delta_{h}{\cal
N}_{h}^{n-\theta}-\Delta_{h}({\cal T}_{f})_{h}^{n-\theta}+\Delta_{h}({\cal
T}_{v1})_{h}^{n-\theta}-\Delta_{h}({\cal T}_{w})_{h}^{n-\theta},$ (34)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}{\bf
x}_{h}\in\Omega_{h},~{}1\leq n\leq N;$ (35)
$\displaystyle\frac{1}{\nu^{2}}({\cal
D}_{\tau}^{\beta}{\tilde{w}}_{h})^{n-\theta}=\Delta_{h}{\tilde{v}}_{h}^{n-\theta}-\Delta_{h}({\cal
T}_{u})_{h}^{n-\theta}+\Delta_{h}({\cal T}_{v2})_{h}^{n-\theta}+({\cal
D}_{\tau}^{\beta}{\cal S}_{h})^{n-\theta},\quad{\bf
x}_{h}\in\Omega_{h},~{}1\leq n\leq N;$
$\displaystyle{\tilde{u}}_{h}^{0}={\tilde{v}}_{h}^{0}={\tilde{w}}_{h}^{0}=0,\quad{\bf
x}_{h}\in\bar{\Omega}_{h};\qquad{\tilde{u}}_{h}^{n}={\tilde{v}}_{h}^{n}=0,\quad{\bf
x}_{h}\in\partial\Omega_{h},~{}1\leq n\leq N.$

######  Lemma 3.4.

Let ${\cal F}(\psi({\bf x}),{\bf x})\in C^{2}(\mathbb{R}\times\Omega)$, and
$\\{\psi_{h}\\}$ be a grid function which satisfy
$\max\\{\|\psi\|_{\infty},\|\nabla_{h}\psi\|,\|\Delta_{h}\psi\|\\}\leq
C_{\psi}$. Then there is a constant $C_{F}>0$ dependent on $C_{\psi}$ and
$C_{\Omega}$ such that

$\|\Delta_{h}[{\cal F}(\psi,{\bf x})v]\|\leq C_{F}\|\Delta_{h}v\|.$

######  Proof 3.5.

The proof can be worked out following that of [20, Lemma 4.1] just by routine
computations on $\delta_{x}{\cal
F}(\psi_{i-\frac{1}{2},j},(x_{i-\frac{1}{2}},y_{j}))$, $\delta_{y}{\cal
F}(\psi_{i,j-\frac{1}{2}},(x_{i},y_{j-\frac{1}{2}}))$, $\delta_{x}^{2}{\cal
F}(\psi_{i,j},(x_{i},y_{j}))$ and $\delta_{y}^{2}{\cal
F}(\psi_{i,j},(x_{i},y_{j}))$ using the Taylor formula with integral
remainder.

We next show the unconditional convergence of the proposed linearized scheme
(22)–(24) based on the $H^{2}$ energy method ([20, 18]).

######  Theorem 3.6.

Let $f\in C^{(4,2,0)}({\mathbb{R}}\times\Omega\times[0,T])$. If the
assumptions in (5) and the mesh assumption MA hold, the linearized schemes
(22)–(24) are unconditional convergent with

(36)
$\|\Delta_{h}{\tilde{u}}^{n}\|+\|\nabla_{h}{\tilde{v}}^{n}\|\leq\left\\{\begin{array}[]{ll}\vskip
3.0pt plus 1.0pt minus
1.0ptC(\tau^{\min\\{2-\beta,\gamma\sigma_{1},\gamma\sigma_{2}\\}}+h^{2}),\quad\mbox{if}\quad\theta=0;\\\
C(\tau^{\min\\{2,\gamma\sigma_{1},\gamma\sigma_{2}\\}}+h^{2}),\quad\mbox{if}\quad\theta=\frac{\beta}{2};\end{array}\quad\mbox{for}\quad
1\leq n\leq N.\right.$

######  Proof 3.7.

Taking inner product of equations (34) and (35) with
${\tilde{v}}_{h}^{n-\theta}$ and ${\tilde{w}}_{h}^{n-\theta}$ respectively, we
have

$\displaystyle\left\langle\Delta_{h}{\tilde{w}}^{n-\theta},{\tilde{v}}^{n-\theta}\right\rangle=$
$\displaystyle\left\langle({\cal
D}_{\tau}^{\beta}\Delta_{h}{\tilde{v}})^{n-\theta},{\tilde{v}}^{n-\theta}\right\rangle$
(37) $\displaystyle+\left\langle\Delta_{h}{\cal N}^{n-\theta}-\Delta_{h}({\cal
T}_{f})^{n-\theta}+\Delta_{h}({\cal T}_{v1})^{n-\theta}-\Delta_{h}({\cal
T}_{w})^{n-\theta},{\tilde{v}}^{n-\theta}\right\rangle$

and

(38) $\displaystyle\left\langle\frac{1}{\nu^{2}}({\cal
D}_{\tau}^{\beta}{\tilde{w}})^{n-\theta},{\tilde{w}}^{n-\theta}\right\rangle=\left\langle\Delta_{h}{\tilde{v}}^{n-\theta},{\tilde{w}}^{n-\theta}\right\rangle+\left\langle-\Delta_{h}({\cal
T}_{u})^{n-\theta}+\Delta_{h}({\cal T}_{v2})^{n-\theta}+({\cal
D}_{\tau}^{\beta}{\cal S})^{n-\theta},{\tilde{w}}^{n-\theta}\right\rangle.$

With the identity
$\left\langle\Delta_{h}{\tilde{w}}^{n-\theta},{\tilde{v}}^{n-\theta}\right\rangle=\left\langle\Delta_{h}{\tilde{v}}^{n-\theta},{\tilde{w}}^{n-\theta}\right\rangle$
and the zero boundary conditions of ${\tilde{v}}_{h}^{n-\theta}$, it follows
form (37)–(38) that

$\displaystyle\left\langle\frac{1}{\nu^{2}}({\cal
D}_{\tau}^{\beta}{\tilde{w}})^{n-\theta},{\tilde{w}}^{n-\theta}\right\rangle+\left\langle({\cal
D}_{\tau}^{\beta}\nabla_{h}{\tilde{v}})^{n-\theta},\nabla_{h}{\tilde{v}}^{n-\theta}\right\rangle$
$\displaystyle=$ $\displaystyle\left\langle\Delta_{h}{\cal
N}^{n-\theta}-\Delta_{h}({\cal T}_{f})^{n-\theta}+\Delta_{h}({\cal
T}_{v1})^{n-\theta}-\Delta_{h}({\cal
T}_{w})^{n-\theta},{\tilde{v}}^{n-\theta}\right\rangle$
$\displaystyle+\left\langle-\Delta_{h}({\cal
T}_{u})^{n-\theta}+\Delta_{h}({\cal T}_{v2})^{n-\theta}+({\cal
D}_{\tau}^{\beta}{\cal S})^{n-\theta},{\tilde{w}}^{n-\theta}\right\rangle.$

Utilizing (18) and the Cauchy-Schwarz inequality, the above equation leads to

$\displaystyle\sum_{k=1}^{n}A_{n-k}^{(n)}\nabla_{\tau}(\|{\tilde{w}}^{n}\|^{2}+\nu^{2}\|\nabla_{h}{\tilde{v}}^{n}\|^{2})$
$\displaystyle\leq$ $\displaystyle 2\nu^{2}\left(\|\Delta_{h}{\cal
N}^{n-\theta}\|+\|\Delta_{h}({\cal T}_{f})^{n-\theta}\|+\|\Delta_{h}({\cal
T}_{v1})^{n-\theta}\|+\|\Delta_{h}({\cal
T}_{w})^{n-\theta}\|\right)\|{\tilde{v}}^{n-\theta}\|$
$\displaystyle+2\nu^{2}\left(\|\Delta_{h}({\cal
T}_{u})^{n-\theta}\|+\|\Delta_{h}({\cal T}_{v2})^{n-\theta}\|+\|({\cal
D}_{\tau}^{\beta}{\cal S})^{n-\theta}\|\right)\|{\tilde{w}}^{n-\theta}\|$ (39)
$\displaystyle\leq$ $\displaystyle
2\nu(1+\nu)(1+C_{\Omega})\left(\|\Delta_{h}{\cal N}^{n-\theta}\|+{\cal
T}^{n-\theta}\right)\left(\|{\tilde{w}}^{n-\theta}\|+\nu\|\nabla_{h}{\tilde{v}}^{n-\theta}\|\right),$

where

${\cal T}^{n-\theta}:=\|\Delta_{h}({\cal
T}_{f})^{n-\theta}\|+\|\Delta_{h}({\cal
T}_{v1})^{n-\theta}\|+\|\Delta_{h}({\cal
T}_{w})^{n-\theta}\|+\|\Delta_{h}({\cal
T}_{u})^{n-\theta}\|+\|\Delta_{h}({\cal T}_{v2})^{n-\theta}\|+\|({\cal
D}_{\tau}^{\beta}{\cal S})^{n-\theta}\|.$

From (45) and (47)–(52), there exist positive constant $C_{r}$ such that

(40) $\sum_{j=1}^{n}P_{n-j}^{(n)}{\cal
T}^{j-\theta}\leq\left\\{\begin{array}[]{ll}\vskip 3.0pt plus 1.0pt minus
1.0ptC_{r}(\tau^{\min\\{2-\beta,\gamma\sigma_{1},\gamma\sigma_{2}\\}}+h^{2}),\quad\mbox{for}\quad\theta=0;\\\
C_{r}(\tau^{\min\\{2,\gamma\sigma_{1},\gamma\sigma_{2}\\}}+h^{2}),\quad\mbox{for}\quad\theta=\frac{\beta}{2}.\end{array}\right.$

From the regularity assumptions in (5), we introduce the following constant

(41) $\displaystyle C_{0}=\max_{0\leq n\leq
N}\\{\|U^{n}\|_{\infty},\|\nabla_{h}U^{n}\|,\|\Delta_{h}U^{n}\|\\}.$

The mathematical induction method will be applied to show that

(42)
$\displaystyle\|{\tilde{w}}^{n}\|+\nu\|\nabla_{h}{\tilde{v}}^{n}\|\leq{\cal
E}_{n}{\tilde{\cal T}}^{n-\theta},\quad 1\leq n\leq N,$

where ${\cal E}_{n}:=4E_{\beta}(2\max(1,\rho)\pi_{A}\Lambda t_{n}^{\beta})$
with $\Lambda=2\nu(1+\nu)C_{f}(1+C_{\Omega})$ and

${\tilde{\cal
T}}^{n-\theta}:=2\nu(1+\nu)(1+C_{\Omega})\times\left\\{\begin{array}[]{ll}\vskip
3.0pt plus 1.0pt minus
1.0ptC_{r}(\tau^{\min\\{2-\beta,\gamma\sigma_{1},\gamma\sigma_{2}\\}}+h^{2})+\pi_{A}\Gamma(1-\beta)C_{f}C_{u}t_{n}^{\beta}h^{2},\quad\mbox{for}\quad\theta=0;\\\
C_{r}(\tau^{\min\\{2,\gamma\sigma_{1},\gamma\sigma_{2}\\}}+h^{2})+\pi_{A}\Gamma(1-\beta)C_{f}C_{u}t_{n}^{\beta}h^{2},\quad\mbox{for}\quad\theta=\frac{\beta}{2},\end{array}\right.$

in which
$C_{f}=\max\\{(1-\theta)C_{1},C_{2}[1+((1-\theta)(C_{2}+1)+1)/\theta]\\}$ with
$C_{1}$ and $C_{2}$ being two proper positive constants which depend on
$C_{0}$ and $C_{\Omega}$.

While $n=1$, it holds that ${\tilde{u}}_{h}^{0}=0$ and
$u_{h}^{0}=U_{h}^{0}\leq C_{0}$. Suppose $f\in
C^{(3,2,0)}({\mathbb{R}}\times\Omega\times[0,T])$, by Lemma 3.4 and (46),
there exists a positive constant $C_{1}$ such that

(43) $\displaystyle\|\Delta_{h}{\cal
N}^{1-\theta}\|=(1-\theta)\|\Delta_{h}f_{u}^{\prime}(u_{h}^{0},{\bf
x},t_{1-\theta}){\tilde{u}}_{h}^{1}\|\leq(1-\theta)C_{1}\|\Delta_{h}{\tilde{u}}_{h}^{1}\|\leq(1-\theta)C_{1}(\|{\tilde{w}}^{1}\|+C_{u}h^{2}).$

For simplicity, denote
$\|{\tilde{w}}^{(n-\theta)}\|:=(1-\theta)\|{\tilde{w}}^{n}\|+\theta\|{\tilde{w}}^{n-1}\|$.
Similarly, we define $\|\nabla_{h}{\tilde{v}}^{(n-\theta)}\|$. The triangle
inequality gives $\|{\tilde{w}}^{n-\theta}\|\leq\|{\tilde{w}}^{(n-\theta)}\|$
and
$\|\nabla_{h}{\tilde{v}}^{n-\theta}\|\leq\|\nabla_{h}{\tilde{v}}^{(n-\theta)}\|$.

Then, it follows from (39) and (43) that

$\displaystyle
A_{0}^{(1)}\nabla_{\tau}(\|{\tilde{w}}^{1}\|^{2}+\nu^{2}\|\nabla_{h}{\tilde{v}}^{1}\|^{2})\leq
2\nu(1+\nu)C_{1}(1+C_{\Omega})\left(\|{\tilde{w}}^{(1-\theta)}\|+\nu\|\nabla_{h}{\tilde{v}}^{(1-\theta)}\|\right)^{2}$
$\displaystyle+2\nu(1+\nu)(1+C_{\Omega})\left({\cal
T}^{1-\theta}+(1-\theta)C_{1}C_{u}h^{2}\right)\left(\|{\tilde{w}}^{(1-\theta)}\|+\nu\|\nabla_{h}{\tilde{v}}^{(1-\theta)}\|\right).$

Thus, applying Lemma 3.2 on the above inequality, and utilizing (40), we get

$\displaystyle\|{\tilde{w}}^{1}\|+\nu\|\nabla_{h}{\tilde{v}}^{1}\|\leq{\cal
E}_{1}\left[2\nu(1+\nu)(1+C_{\Omega})P_{0}^{(1)}\left({\cal
T}^{1-\theta}+(1-\theta)C_{1}C_{u}h^{2}\right)\right]\leq{\cal
E}_{1}{\tilde{\cal T}}^{1-\theta},$

which means that (42) holds for $n=1$.

Assume that (42) is valid for $1\leq k\leq n-1~{}(n\geq 2)$. The eq. (30) and
discrete embedding inequalities imply that

$\displaystyle\max\\{\|{\tilde{u}}^{k}\|_{\infty},\|\nabla_{h}{\tilde{u}}^{k}\|,\|\Delta_{h}{\tilde{u}}^{k}\|\\}\leq$
$\displaystyle\max\\{1,C_{\Omega}\\}\left(\frac{1}{\nu^{2}}\|{\tilde{w}}^{k}\|+\|{\cal
S}^{k}\|\right)$ $\displaystyle\leq$
$\displaystyle\max\\{1,C_{\Omega}\\}\left(\frac{1}{\nu^{2}}{\cal
E}_{k}{\tilde{\cal T}}^{k-\theta}+C_{u}h^{2}\right)\leq 1,$

for $1\leq k\leq n-1$ and small step sizes. So according to (41), the
numerical solutions satisfy

$\max\\{\|{u}^{k}\|_{\infty},\|\nabla_{h}{u}^{k}\|,\|\Delta_{h}{u}^{k}\|\\}\leq
C_{0}+1.$

Now for $k=n$, suppose $f\in
C^{(4,2,0)}({\mathbb{R}}\times\Omega\times[0,T])$. By Lemma 3.4 there exists a
positive constant $C_{2}$ such that

$\displaystyle\|\Delta_{h}{\cal N}^{n-\theta}\|\leq$
$\displaystyle(1-\theta)\bigg{(}\|\Delta_{h}[f_{u}^{\prime}(u^{n-1},{\bf
x},t_{n-\theta})\nabla_{\tau}{\tilde{u}}^{n}]\|$
$\displaystyle+\int_{0}^{1}\|\Delta_{h}[f_{u}^{\prime\prime}(sU^{n-1}+(1-s)u^{n-1},{\bf
x},t_{n-\theta}){\tilde{u}}^{n-1}\nabla_{\tau}U^{n}]\|\,\mathrm{d}s\bigg{)}$
$\displaystyle+\int_{0}^{1}\|\Delta_{h}[f_{u}^{\prime}(sU^{n-1}+(1-s)u^{n-1},{\bf
x},t_{n-\theta}){\tilde{u}}^{n-1}]\|\,\mathrm{d}s$ $\displaystyle\leq$
$\displaystyle(1-\theta)C_{2}\left(\|\Delta_{h}(\nabla_{\tau}{\tilde{u}}^{n})\|+C_{2}\|\Delta_{h}{\tilde{u}}^{n-1}\|\right)+C_{2}\|\Delta_{h}{\tilde{u}}^{n-1}\|$
$\displaystyle\leq$ $\displaystyle
C_{f}\left[(1-\theta)\|\Delta_{h}{\tilde{u}}^{n}\|+\theta\|\Delta_{h}{\tilde{u}}^{n-1}\|\right]$
(44) $\displaystyle\leq$ $\displaystyle
C_{f}\|{\tilde{w}}^{(n-\theta)}\|+C_{f}C_{u}h^{2}.$

So (39) and (44) lead to

$\displaystyle\sum_{k=1}^{n}A_{n-k}^{(n)}\nabla_{\tau}(\|{\tilde{w}}^{n}\|^{2}+\nu^{2}\|\nabla_{h}{\tilde{v}}^{n}\|^{2})\leq
2\nu(1+\nu)C_{f}(1+C_{\Omega})\left(\|{\tilde{w}}^{(n-\theta)}\|+\nu\|\nabla_{h}{\tilde{v}}^{(n-\theta)}\|\right)^{2}$
$\displaystyle+2\nu(1+\nu)(1+C_{\Omega})\left({\cal
T}^{n-\theta}+C_{f}C_{u}h^{2}\right)\left(\|{\tilde{w}}^{(n-\theta)}\|+\nu\|\nabla_{h}{\tilde{v}}^{(n-\theta)}\|\right).$

Applying Lemma 3.2 and utilizing (40) again, it yields

$\displaystyle\|{\tilde{w}}^{n}\|+\nu\|\nabla_{h}{\tilde{v}}^{n}\|\leq{\cal
E}_{n}\left[2\nu(1+\nu)(1+C_{\Omega})\max_{1\leq k\leq
n}\sum_{j=1}^{k}P_{k-j}^{(k)}\left({\cal
T}^{j-\theta}+C_{f}C_{u}h^{2}\right)\right]\leq{\cal E}_{n}{\tilde{\cal
T}}^{n-\theta}.$

Therefore (42) is verified.

Finally, the desired result (36) is reached by (30) and unifying the
constants.

######  Remark 3.8.

A memory and computational storage saving technique (called SOE approximation)
investigated in [8] (see also [8, Theorem 2.5] or [19, Lemma 5.1]) to compute
the discrete Caputo derivative can be directly employed to the nonuniform L1
and Alikhanov formulas, the corresponding coefficients of fast L1 and fast
Alikhanov formulas preserve the properties A1–A2 [17, 19] which further ensure
the theoretical analysis of the associated fast schemes. Therefore, in our
later implementation of adaptive time stepping strategy and numerical tests,
we will always utilize the fast L1 formula [17, Example 2] and the fast
Alikhanov formula [19, eq. (5.3)] while applying the proposed algorithms
(22)–(24) with $\theta=0$ and $\theta=\beta/2$, respectively.

###  3.4 Adaptive time-stepping strategy

The time mesh assumption in Theorems 3.6 permits us to establish adaptive
time-stepping strategy based on the fast L1 and fast Alikhanov algorithms to
reduce the computational costs while solving the semilinear diffusion-wave
equation, especially when the solution of the governing problem may possess
highly oscillatory feature in time. In the following, we refer to [7, 19] for
designing an adaptive time-stepping algorithm of the semilinear diffusion-wave
equation (1), the strategy is presented in Algorithm 1.

Algorithm 1 Adaptive time-stepping strategy

Given: $u^{n}$, $v^{n}$ and time step $\tau_{n+1}$

1: Compute $u_{1}^{n+1}$ by the (fast) L1 scheme ((22)–(24) for $\theta=0$)
with time step $\tau_{n+1}$;

2: Compute $u_{2}^{n+1}$ by the (fast) Alikhanov scheme ((22)–(24) for
$\theta=\beta/2$) with time step $\tau_{n+1}$;

3: Calculate $e^{n+1}=\|u_{2}^{n+1}-u_{1}^{n+1}\|/\|u_{2}^{n+1}\|$;

4: if $e^{n+1}<tol$ or $\tau_{n+1}=\frac{2}{3}\tau_{n}$ then

5: Update time step size
$\tau_{n+2}\leftarrow\min\\{\max\\{\tau_{\min},\tau_{ada}\\},\tau_{\max}\\}$;

6: else

7: Reformulate the time-step size
$\tau_{n+1}\leftarrow\max\\{\min\\{\max\\{\tau_{\min},\tau_{ada}\\},\tau_{\max}\\},\frac{2}{3}\tau_{n+1}\\}$;

8: Goto 1

9: end if

The adaptive time step size in Algorithm 1 is updated by

$\tau_{ada}(e,\tau)=S\left(\frac{tol}{e}\right)^{\frac{1}{2}}\tau,$

where $S$, $tol$ denote the safety coefficient and the tolerance,
respectively.

##  4 Numerical experiments

Numerical examples are carried out in this section to show the accuracy and
efficiency of proposed algorithms. The absolute tolerance error $\epsilon$ and
the cut-off time $\Delta t$ of fast L1 formula [17, Example 2] and the fast
Alikhanov formula [19, Lemma 5.1] are set as $\epsilon=10^{-12}$ and $\Delta
t=\tau_{1}$ in all of the following tests.

######  Example 4.1.

We first consider the problem (1) with $\Omega=(0,1)^{2}$, $T=1$, $\nu=1$ and

$f(u,{\bf x},t)=-u^{3}+[\sin(\pi x)\sin(\pi y)(1+t+t^{\alpha})]^{3}+\sin(\pi
x)\sin(\pi y)\left[\Gamma(\alpha+1)+2\pi^{2}(1+t+t^{\alpha})\right].$

In this situation, the exact solution is $u=\sin(\pi x)\sin(\pi
y)(1+t+t^{\alpha})$.

One may notice that the regularity parameters in (5) are $\sigma_{1}=\alpha$
and $\sigma_{2}=\alpha/2$ for Example 4.1. Therefore, according to Theorem
3.6, the optimal mesh parameter is $\gamma_{opt}=(4-\alpha)/\alpha$ for the
nonuniform L1 scheme and takes the value $\gamma_{opt}=4/\alpha$ for the
nonuniform Alikhanov scheme. They are all bounded for $\alpha\in(1,2)$. These
bounded mesh parameters keep the robustness of the algorithms in practical
implementation if the graded mesh $t_{k}=T(k/N)^{\gamma}$ is imposed to deal
with the weak initial singularity. On the other hand, the optimal grading
parameters of the nonuniform schemes in [12, 15, 16, 19, 20, 35] will grow
without bound while the fractional order becomes small as they are all
possessing the form $\gamma_{opt}=r/\alpha$ where $\alpha\in(0,1)$ for the
sub-diffusion problems and $r$ should be the optimal time rate, this generally
lead to practical limitations.

Since the spatial error ${\cal O}(h^{2})$ is standard, we only display the
temporal accuracy of the fast L1 and fast Alikhanov schemes. For Example 4.1,
we fixed a fine spatial grid mesh with $M=1000$ such that the temporal errors
dominate the spatial errors. In each tests, the time interval $[0,T]$ is
divided into two parts $[0,T_{0}]$ and $(T_{0},T]$ with total $N$ time nodes.
A graded mesh with $t_{k}=T_{0}(k/N_{0})^{\gamma}$ in the first interval
$[0,T_{0}]$ is utilized to resolve the weak initial singularity, where
$T_{0}=\min\\{1/\gamma,T\\}$. For the second interval $(T_{0},T]$, we use
random time-step sizes
$\tau_{N_{0}+k}=(T-T_{0})\epsilon_{k}/\sum_{k=1}^{N_{1}}\epsilon_{k}$ for
$N_{1}=N-N_{0}$, where $\epsilon_{k}$ take values in $(0,1)$ randomly. For
this example, we take $N_{0}=\lceil\frac{N}{T+1-\gamma^{-1}}\rceil$. The
discrete $H^{2}$-norm errors $e_{H^{2}}(N)=\max_{1\leq n\leq
N}\|U^{n}-u^{n}\|_{H^{2}}$ are recorded in each run, and the temporal
convergence order is given by

$\mbox{Order}=\log_{2}\left[\frac{e_{H^{2}}(N/2)}{e_{H^{2}}(N)}\right].$

Table 1: Numerical accuracy in temporal direction of fast L1 scheme for
Example 4.1, where $\alpha=1.1$.

$\gamma=1$ $\gamma_{opt}=(4-\alpha)/\alpha\approx 2.64$
$\gamma=\frac{9}{8}\gamma_{opt}\approx 2.97$ $N$ $e_{H^{2}}(N)$ Order
$e_{H^{2}}(N)$ Order $e_{H^{2}}(N)$ Order $16$ 4.2668e-02 $\ast$ 1.4285e-01
$\ast$ 2.4493e-01 $\ast$ $32$ 3.3723e-02 0.34 3.5519e-02 2.01 5.1210e-02 2.26
$64$ 2.2386e-02 0.59 1.0731e-02 1.73 1.2611e-02 2.02 $128$ 1.3688e-02 0.71
2.4621e-03 2.12 3.8861e-03 1.70 Theoretical Order 0.55 1.45 1.45

Table 2: Numerical accuracy in temporal direction of fast L1 scheme for
Example 4.1, where $\alpha=1.5$.

$\gamma=1$ $\gamma_{opt}=(4-\alpha)/\alpha\approx 1.67$
$\gamma=\frac{9}{8}\gamma_{opt}\approx 1.88$ $N$ $e_{H^{2}}(N)$ Order
$e_{H^{2}}(N)$ Order $e_{H^{2}}(N)$ Order $16$ 3.4875e-02 $\ast$ 2.2507e-01
$\ast$ 2.5657e-01 $\ast$ $32$ 1.2196e-02 1.52 3.6991e-02 2.61 3.4006e-02 2.92
$64$ 8.7566e-03 0.48 1.2921e-02 1.52 1.3962e-02 1.28 $128$ 5.5637e-03 0.65
4.2362e-03 1.61 3.8805e-03 1.85 Theoretical Order 0.75 1.25 1.25

Table 3: Numerical accuracy in temporal direction of fast L1 scheme for
Example 4.1, where $\alpha=1.9$.

$\gamma=1$ $\gamma_{opt}=(4-\alpha)/\alpha\approx 1.11$
$\gamma=\frac{9}{8}\gamma_{opt}\approx 1.24$ $N$ $e_{H^{2}}(N)$ Order
$e_{H^{2}}(N)$ Order $e_{H^{2}}(N)$ Order $16$ 7.4641e-02 $\ast$ 7.4027e-02
$\ast$ 1.3225e-01 $\ast$ $32$ 3.7415e-02 1.00 3.4234e-02 1.11 3.4250e-02 1.95
$64$ 1.7990e-02 1.06 1.6886e-02 1.02 1.6737e-02 1.03 $128$ 8.4156e-03 1.10
8.0910e-03 1.06 8.0567e-03 1.05 Theoretical Order 0.95 1.05 1.05

Table 4: Numerical accuracy in temporal direction of fast Alikhanov scheme for
Example 4.1, where $\alpha=1.2$.

$\gamma=1$ $\gamma_{opt}=4/\alpha\approx 3.33$
$\gamma=\frac{9}{8}\gamma_{opt}\approx 3.75$ $N$ $e_{H^{2}}(N)$ Order
$e_{H^{2}}(N)$ Order $e_{H^{2}}(N)$ Order $16$ 5.2656e-02 $\ast$ 1.2494e-01
$\ast$ .5496e-01 $\ast$ $32$ 3.2671e-02 0.69 3.3236e-02 1.91 4.1575e-02 1.90
$64$ 2.0683e-02 0.66 8.5962e-03 1.95 1.0801e-02 1.94 $128$ 1.1645e-02 0.83
2.1990e-03 1.97 2.8352e-03 1.93 Theoretical Order 0.60 2.00 2.00

Table 5: Numerical accuracy in temporal direction of fast Alikhanov scheme for
Example 4.1, where $\alpha=1.5$.

$\gamma=1$ $\gamma_{opt}=4/\alpha\approx 2.67$
$\gamma=\frac{9}{8}\gamma_{opt}=3$ $N$ $e_{H^{2}}(N)$ Order $e_{H^{2}}(N)$
Order $e_{H^{2}}(N)$ Order $16$ 3.0823e-02 $\ast$ 7.0440e-02 $\ast$ 8.7416e-02
$\ast$ $32$ 1.3857e-02 1.15 1.8560e-02 1.92 2.3212e-02 1.91 $64$ 6.2024e-03
1.16 4.7736e-03 1.96 5.9919e-03 1.95 $128$ 2.6236e-03 1.24 1.2150e-03 1.97
1.5269e-03 1.97 Theoretical Order 0.75 2.00 2.00

Table 6: Numerical accuracy in temporal direction of fast Alikhanov scheme for
Example 4.1, where $\alpha=1.8$.

$\gamma=1$ $\gamma_{opt}=4/\alpha\approx 2.22$
$\gamma=\frac{9}{8}\gamma_{opt}=2.50$ $N$ $e_{H^{2}}(N)$ Order $e_{H^{2}}(N)$
Order $e_{H^{2}}(N)$ Order $16$ 1.9521e-02 $\ast$ 3.5560e-02 $\ast$ 4.3938e-02
$\ast$ $32$ 6.7203e-03 1.54 9.3089e-03 1.93 1.1590e-02 1.92 $64$ 2.6309e-03
1.35 2.3828e-03 1.97 2.9755e-03 1.96 $128$ 1.1487e-03 1.20 6.0470e-04 1.98
7.5559e-04 1.98 Theoretical Order 0.90 2.00 2.00

Tables 1–3 record the numerical results of the proposed fast L1 scheme with
different grading parameters when solving the example for different $\alpha$.
One can observe that the L1 scheme works accurately with the optimal temporal
convergence of ${\cal
O}(\tau^{\min\\{2-\frac{\alpha}{2},\gamma\frac{\alpha}{2}\\}})$. Similar
numerical tests of the fast Alikhanov scheme are carried out for the example,
and the results are listed in Tables 4–6. The temporal convergence of ${\cal
O}(\tau^{\min\\{2,\gamma\frac{\alpha}{2}\\}})$ is well reflected and the
optimal second-order convergence is apparent while
$\gamma\geq\gamma_{opt}={4}/{\alpha}$.

######  Example 4.2.

Consider the semilinear problem (1) with $\Omega=(-1,1)^{2}$, $\nu=1$ and
$f(u,{\bf x},t)=-u^{3}$. The initial data are given as

$\varphi({\bf
x})=(x^{2}-1)(y^{2}-1)\\{\exp\\{-10((x+0.4)^{2}+y^{2})]+\exp[-10((x-0.4)^{2}+y^{2})]\\},\quad{\tilde{\varphi}}({\bf
x})=0.$

Figure 1: The numerical solution in maximum-norm of the Algorithm 1 and the
Graded-Uniform scheme for Example 4.2 with $\alpha=1.5$.

(a) t=0

(b) t=0.5

(c) t=2

(d) t=10

Figure 2: Contour plots of the solutions of Algorithm 1 for Example 4.2 at
different time with $\alpha=1.5$. Figure 3: The variation of time step sizes
of the Algorithm 1 and the Graded-Uniform scheme for Example 4.2 with
$\alpha=1.5$.

For Example 4.2, we choose the spatial node $M=100$ and also divide the time
interval $[0,T]$ into two parts $[0,T_{0}]$ and $(T_{0},T]$ with $T_{0}=0.02$.
The Alikhanov algorithm on graded mesh with $t_{k}=T_{0}(k/N_{0})^{\gamma}$
($\gamma=\alpha/4$) in the first interval $[0,T_{0}]$ is utilized to resolve
the possible weak initial singularity, where $N_{0}=30$. For the remaining
interval $(T_{0},T]$, we employ the proposed adaptive time stepping strategy
(Algorithm 1) to compute the numerical solution until $T=10$. The parameters
of the adaptive algorithm for solving this example are

$tol=10^{-3},~{}S=0.9,~{}\tau_{\min}=10^{-3},~{}\tau_{\max}=10^{-1},~{}\tau_{N_{0}+1}=\tau_{N_{0}}.$

In order to show the efficiency of the adaptive algorithm, the fast linearized
Alikhanov scheme is applied at the same time to find the solution in the
interval $(T_{0},T]$. Its temporal mesh is graded (with $\gamma=\alpha/4$) in
$[0,T_{0}]$ and is uniform in $(T_{0},T]$. In the following, we use ‘Graded-
Uniform’ to represent this scheme.

Figure 1 displays the numerical solution in maximum-norm of the Algorithm 1
and the Graded-Uniform scheme for $\alpha=1.5$. It implies that the adaptive
mesh suits well with a dense uniform mesh in $(T_{0},T]$, provided that the
adaptive mesh requires 277 time nodes in the remain interval $(T_{0},T]$
whereas the uniform mesh needs 970 time nodes. Figure 2 gives the solution
contour plots of the solutions by the adaptive strategy, which simply shows
the wave interactions of the example at different time. The variation of the
temporal step sizes of the adaptive strategy with its comparison to those of
the Graded-Uniform scheme are presented in Figure 3. The results indicate that
the adaptive time-stepping strategy should be efficient and robust in the long
time simulation of the semilinear diffusion-wave equations especially when the
solution may exhibit high oscillations in time.

##  5 Concluding remarks

We proposed a novel order reduction method to equivalently rewrite the
semilinear diffusion-wave equation into coupled equations, where the explicit
time-fractional derivative orders are all $\alpha/2$. The L1 and Alikhanov
schemes combining with linearized approximations have been constructed for the
equivalent problem. By using $H^{2}$ energy method, unconditional convergences
(Theorem 3.6) were obtained for the two proposed algorithms under reasonable
regularity assumptions and weak mesh restrictions. An adaptive time-stepping
strategy was then designed for the semilinear problem to deal with possible
temporal oscillations of the solution. The theoretical results were well
demonstrated by our numerical experiments.

We finally point out several relevant issues that deserve for further study:
(i) deriving the regularity of the linear and semilinear diffusion-wave
equations for the difference schemes; (ii) studying the energy properties of
the nonlinear diffusion-wave equations for both the continuous and discrete
versions, noting that corresponding properties were investigated recently for
the nonlinear sub-diffusion problems [37]; (iii) extending the proposed
methods to some related problems, such as the multi-term time-fractional wave
equation [23].

##  6 Acknowledgement

The authors are very grateful to Prof. Hong-lin Liao for his great help on the
design of the SFOR method and valuable suggestions on other parts of the whole
paper.

##  7 Appendix: Truncation error analysis

The truncation errors in (29)–(31) are defined as

$\displaystyle({\cal T}_{f})_{h}^{n-\theta}:=f(U_{h}^{n-\theta},{\bf
x}_{h},t_{n-\theta})-\left[f(U_{h}^{n-1},{\bf
x}_{h},t_{n-\theta})+(1-\theta)f^{\prime}(U_{h}^{n-1},{\bf
x}_{h},t_{n-\theta})\nabla_{\tau}U_{h}^{n}\right],$ $\displaystyle({\cal
T}_{v1})_{h}^{n-\theta}:={\cal D}_{t}^{\beta}{\bf v}({\bf
x}_{h},t_{n-\theta})-({\cal D}_{\tau}^{\beta}{\bf v}_{h})^{n-\theta},$
$\displaystyle({\cal T}_{w})_{h}^{n-\theta}:=w({\bf
x}_{h},t_{n-\theta})-w_{h}^{n-\theta},$ $\displaystyle({\cal
T}_{u})_{h}^{n-\theta}:={\cal D}_{t}^{\beta}{\bf u}({\bf
x}_{h},t_{n-\theta})-({\cal D}_{\tau}^{\beta}{\bf u}_{h})^{n-\theta},$
$\displaystyle({\cal T}_{v2})_{h}^{n-\theta}:={\bf v}({\bf
x}_{h},t_{n-\theta})-{\bf v}_{h}^{n-\theta},$ $\displaystyle{\cal
S}_{h}^{n}:=\Delta u({\bf x}_{h},t_{n})-\Delta_{h}u_{h}^{n}.$

According to [16, Lemma 3.8 and Theorem 3.9], we have the follow lemma on
estimating the time weighted approximation.

######  Lemma 7.1.

Assume that $g\in C^{2}((0,T])$ and there exists a constant $C_{g}>0$ such
that

$|g^{\prime\prime}(t)|\leq C_{g}(1+t^{\sigma-2}),\quad 0<t\leq T,$

where $\sigma\in(0,1)\cup(1,2)$ is a regularity parameter. Denote the local
truncation error of $g^{n-\vartheta}$ (here $\vartheta=\beta/2$) by

${\tilde{\cal R}}^{n-\vartheta}=g(t_{n-\vartheta})-g^{n-\vartheta},\quad 1\leq
n\leq N.$

If the mesh assumption MA holds, then

$\displaystyle\sum_{j=1}^{n}P_{n-j}^{(n)}|{\tilde{\cal R}}^{j-\vartheta}|\leq
C_{g}\left(\tau_{1}^{\sigma+\beta}/\sigma+t_{n}^{\beta}\max_{2\leq k\leq
n}t_{k-1}^{\sigma-2}\tau_{k}^{2}\right)\leq C\tau^{\min\\{\gamma\sigma,2\\}}.$

The following lemma is provided to analyze $({\cal T}_{f})_{h}^{n-\theta}$,
which is analogous to Lemma 3.4 in [20].

######  Lemma 7.2.

Assume that $\eta\in C([0,T])\cap C^{2}((0,T])$ and there exists a constant
$C_{u}>0$ such that

$|\eta^{(k)}(t)|\leq C_{u}(1+t^{\sigma-k}),\quad 0<t\leq T,\quad k=1,2,$

where $\sigma\in(0,1)\cup(1,2)$ is a regularity parameter. Assume further that
the nonlinear function $f(u,x,t)\in C^{4}({\mathbb{R}})$ with respect to $u$.
Denote $\eta^{n}=\eta(t_{n})$ and the local truncation error

$\displaystyle{\cal
R}_{f}^{n-\theta}:=f(\eta(t_{n-\theta}),x,t)-\left[f(\eta^{n-1},x,t)+(1-\theta)f^{\prime}_{u}(\eta^{n-1},x,t)\nabla_{\tau}\eta^{n}\right].$

If the assumption MA holds, then

$\sum_{j=1}^{n}P_{n-j}^{(n)}|{\cal
R}_{f}^{j-\theta}|\leq\left\\{\begin{array}[]{ll}C\tau^{\min\\{2\gamma\sigma,2\\}},\quad\theta=0;\\\
C\tau^{\min\\{\gamma\sigma,2\\}},\quad\theta=\frac{\beta}{2};\end{array}\quad
1\leq n\leq N.\right.$

######  Proof 7.3.

Denote ${\cal R}_{\eta}^{n-\theta}:=\eta(t_{n-\theta})-\eta^{n-\theta}$. We
have ${\cal R}_{\eta}^{n-\theta}=0$ while $\theta=0$. By the Taylor expansion,
we have

$\displaystyle{\cal R}_{f}^{n-\theta}=$ $\displaystyle
f_{u}^{\prime}(\eta^{n-1}){\cal R}_{\eta}^{n-\theta}$
$\displaystyle+\left((1-\theta)\nabla_{\tau}\eta^{n}+{\cal
R}_{\eta}^{n-\theta}\right)^{2}\int_{0}^{1}f^{\prime\prime}_{u}\left(\eta^{n-1}+s(\eta(t_{n-\theta})-\eta^{n-1}),x,t\right)(1-s)\,\mathrm{d}s.$

Following the proof of [20, Lemma 3.4] and using Lemma 7.1, the desired result
holds immediately.

For ${\bf x}\in\Omega$, let $\xi^{n}({\bf x})$ be a spatially continues
function and denote $\xi^{n}_{h}:=\xi^{n}({\bf x}_{h})$. One may apply the
Taylor expansion to get

$\displaystyle\Delta_{h}\xi^{n}_{h}=$
$\displaystyle\int_{0}^{1}\left[\partial_{xx}\xi^{n}(x_{i}-sh_{x},y_{j})+\partial_{xx}\xi^{n}(x_{i}+sh_{x},y_{j})\right](1-s)\,\mathrm{d}s$
$\displaystyle+\int_{0}^{1}\left[\partial_{yy}\xi^{n}(x_{i},y_{j}-sh_{y})+\partial_{yy}\xi^{n}(x_{i},y_{j}+sh_{y})\right](1-s)\,\mathrm{d}s,\quad
1\leq n\leq N.$

Then we define a function ${\cal T}_{f}^{n}({\bf x})$ by $({\cal
T}_{f})_{h}^{n}={\cal T}_{f}^{n}({\bf x}_{h})$. If the assumptions in (5) and
MA are satisfied, by Lemma 7.2 and the differential formula of composite
function, we can obtain

(45) $\sum_{j=1}^{n}P_{n-j}^{(n)}\|\Delta_{h}({\cal
T}_{f})^{n-\theta}\|\leq\left\\{\begin{array}[]{ll}C\tau^{\min\\{2\gamma\sigma_{1},2\\}},\quad\theta=0;\\\
C\tau^{\min\\{\gamma\sigma_{1},2\\}},\quad\theta=\frac{\beta}{2};\end{array}\quad
1\leq n\leq N.\right.$

For the spatial error, based on the regularity condition, it is easy to know
that

(46) $\displaystyle\|{\cal S}^{n}\|\leq C_{u}h^{2},\quad 1\leq n\leq N.$

Then

(47) $\displaystyle\sum_{j=1}^{n}P_{n-j}^{(n)}\|({\cal D}_{\tau}^{\beta}{\cal
S})^{j}\|\leq\sum_{j=1}^{n}P_{n-j}^{(n)}\sum_{k=1}^{j}A_{j-k}^{(j)}\|\nabla_{\tau}{\cal
S}^{k}\|=\sum_{k=1}^{n}\|\nabla_{\tau}{\cal S}^{k}\|\leq
C_{u}(1+t_{n}^{\sigma_{1}-1})h^{2}.$

We now consider the temporal truncation errors $({\cal
T}_{v1})_{h}^{n-\theta}$, $({\cal T}_{w})_{h}^{n-\theta}$, $({\cal
T}_{u})_{h}^{n-\theta}$ and $({\cal T}_{v2})_{h}^{n-\theta}$ in two
situations: $\theta=0$ and $\theta=\beta/2$.

For a function $g(t)$, define the global error

${\cal R}^{n-\theta}:=({\cal D}_{t}^{\beta}g)(t_{n-\theta})-({\cal
D}_{\tau}^{\beta}g)^{n-\theta},\quad 1\leq n\leq N.$

$\bullet$ For L1 approximation ($\theta=0$): We have $({\cal
T}_{w})_{h}^{n}=({\cal T}_{v2})_{h}^{n}=0$ in this situation.

According to [15, Lemma 3.3] and [20, Lemma 3.3], the global consistency error
of the L1 approximation can be presented in the following lemma.

######  Lemma 7.4.

Assume that $g\in C^{2}((0,T])$ and there exists a constant $C_{g}>0$ such
that

$|g^{\prime\prime}(t)|\leq C_{g}(1+t^{\sigma-2}),\quad 0<t\leq T,$

where $\sigma\in(0,1)\cup(1,2)$ is a regularity parameter. If the assumption
MA holds, it follows that

$\displaystyle\sum_{j=1}^{n}P_{n-j}^{(n)}|{\cal R}^{j}|\leq
C_{g}\left(\tau_{1}^{\sigma}/\sigma+\frac{1}{1-\beta}\max_{2\leq k\leq
n}(t_{k}-t_{1})^{\beta}t_{k-1}^{\sigma-2}\tau_{k}^{2-\beta}\right)\leq
C\tau^{\min\\{2-\beta,\gamma\sigma\\}}.$

Define the functions ${\cal T}_{v1}^{n}({\bf x})$ and ${\cal T}_{u}^{n}({\bf
x})$ by $({\cal T}_{v1})_{h}^{n}:={\cal T}_{v1}^{n}({\bf x}_{h})$ and $({\cal
T}_{u})_{h}^{n}:={\cal T}_{u}^{n}({\bf x}_{h})$ respectively. Using similar
techniques for (45), and Lemma 7.4 with the in assumptions (5) and MA, we have

(48) $\displaystyle\sum_{j=1}^{n}P_{n-j}^{(n)}\|\Delta_{h}{\cal
T}_{v1}^{n}\|\leq
C\tau^{\min\\{2-\beta,\gamma\sigma_{2}\\}}\quad\mbox{and}\quad\sum_{j=1}^{n}P_{n-j}^{(n)}\|\Delta_{h}{\cal
T}_{u}^{n}\|\leq C\tau^{\min\\{2-\beta,\gamma\sigma_{1}\\}}.$

$\bullet$ For Alikhanov approximation ($\theta=\beta/2$):

The global consistency error estimate of the Alikhanov approximation is
estimated in the next lemma.

######  Lemma 7.5.

([16, Lemma 3.6]) Assume that $g\in C^{3}((0,T])$ and there exists a constant
$C_{g}>0$ such that

$|g^{\prime\prime\prime}(t)|\leq C_{g}(1+t^{\sigma-3}),\quad 0<t\leq T,$

where $\sigma\in(0,1)\cup(1,2)$ is a regularity parameter. Then

$\displaystyle\sum_{j=1}^{n}P_{n-j}^{(n)}|{\cal R}^{j-\theta}|\leq
C_{g}\left(\tau_{1}^{\sigma}/\sigma+t_{1}^{\sigma-3}\tau_{2}^{3}+\frac{1}{1-\beta}\max_{2\leq
k\leq
n}t_{k}^{\beta}t_{k-1}^{\sigma-3}\tau_{k}^{3}/\tau_{k-1}^{\beta}\right).$

By Lemma 7.5, Lemma 7.1, the assumptions in (5) and MA, it is easy to get that

(49) $\displaystyle\sum_{j=1}^{n}P_{n-j}^{(n)}\|\Delta_{h}({\cal
T}_{v1})_{h}^{j-\theta}\|\leq
C\left(\tau_{1}^{\sigma_{2}}+\tau_{2}^{3}\tau_{1}^{\sigma_{2}-3}+\max_{2\leq
k\leq
n}(t_{k}-t_{1})^{\beta}t_{k-1}^{\sigma_{2}-3}\tau_{k}^{3-\beta}\right)\leq
C\tau^{\min\\{3-\beta,\gamma\sigma_{2}\\}},$ (50)
$\displaystyle\sum_{j=1}^{n}P_{n-j}^{(n)}\|\Delta_{h}({\cal
T}_{u})_{h}^{j-\theta}\|\leq
C\left(\tau_{1}^{\sigma_{1}}+\tau_{2}^{3}\tau_{1}^{\sigma_{1}-3}+\max_{2\leq
k\leq
n}(t_{k}-t_{1})^{\beta}t_{k-1}^{\sigma_{1}-3}\tau_{k}^{3-\beta}\right)\leq
C\tau^{\min\\{3-\beta,\gamma\sigma_{1}\\}},$ (51)
$\displaystyle\sum_{j=1}^{n}P_{n-j}^{(n)}\|\Delta_{h}({\cal
T}_{w})_{h}^{j-\theta}\|\leq C\left(\tau_{1}^{\sigma_{1}+\beta}+\max_{2\leq
k\leq n}t_{k-1}^{\sigma_{1}-2}\tau_{k}^{2}\right)\leq
C\tau^{\min\\{2,\gamma\sigma_{1}\\}},$ (52)
$\displaystyle\sum_{j=1}^{n}P_{n-j}^{(n)}\|\Delta_{h}({\cal
T}_{v2})_{h}^{j-\theta}\|\leq C\left(\tau_{1}^{\sigma_{2}+\beta}+\max_{2\leq
k\leq n}t_{k-1}^{\sigma_{2}-2}\tau_{k}^{2}\right)\leq
C\tau^{\min\\{2,\gamma\sigma_{2}\\}}.$

## References

  * [1] A. A. Alikhanov, A new difference scheme for the time fractional diffusion equation, J. Comput, Phys., 280 (2015), pp. 424–438, https://doi.org/10.1016/j.jcp.2014.09.031.
  * [2] H. Chen and M. Stynes, Error analysis of a second-order method on fitted meshes for a time-fractional diffusion problem, J. Sci. Comput., 79 (2019), pp. 624–647, https://doi.org/10.1007/s10915-018-0863-y.
  * [3] E. Cuesta, M. Kirane, and S. A. Malik, Image structure preserving denoising using generalized fractional time integrals, Signal Processing, 92 (2012), pp. 553–563, https://doi.org/10.1016/j.sigpro.2011.09.001.
  * [4] E. Cuesta, C. Lubich, and C. Palencia, Convolution quadrature time discretization of fractional diffusion-wave equations, Math. Comput., 75 (2006), pp. 673–696, https://doi.org/10.1090/S0025-5718-06-01788-1.
  * [5] E. Cuesta and C. Palencia, A fractional trapezoidal rule for integro-differential equations of fractional order in banach spaces, Appl. Numer. Math., 45 (2003), pp. 139–159, https://doi.org/10.1016/S0168-9274(02)00186-1.
  * [6] E. Cuesta and C. Palencia, A numerical method for an integro-differential equation with memory in Banach spaces: qualitative properties, SIAM J. Numer. Anal., 41 (2003), pp. 1232–1241, https://doi.org/10.1137/S0036142902402481.
  * [7] H. Gomez and T. J. R. Hughes, Provably unconditionally stable, second-order time-accurate, mixed variational methods for phase-field models, J. Comput. Phys., 230 (2011), pp. 5310–5327, https://doi.org/10.1016/j.jcp.2011.03.033.
  * [8] S. Jiang, J. Zhang, Q. Zhang, and Z. Zhang, Fast evaluation of the Caputo fractional derivative and its applications to fractional diffusion equations, Commun. Comput. Phys., 21 (2017), pp. 650–678, https://doi.org/10.4208/cicp.OA-2016-0136.
  * [9] B. Jin, R. Lazarov, and Z. Zhou, An analysis of the L1 scheme for the subdiffusion equation with nonsmooth data, IMA J. Numer. Anal., 36 (2016), pp. 197–221, https://doi.org/10.1093/imanum/dru063.
  * [10] B. Jin, R. Lazarov, and Z. Zhou, Two fully discrete schemes for fractional diffusion and diffusion-wave equations with nonsmooth data, SIAM J. Sci. Comput., 38 (2016), pp. A146–A170, https://doi.org/10.1137/140979563.
  * [11] B. Jin, B. Li, and Z. Zhou, Discrete maximal regularity of time-stepping schemes for fractional evolution equations, Numer. Math., 138 (2018), pp. 101–131, https://doi.org/10.1007/s00211-017-0904-8.
  * [12] N. Kopteva, Error analysis of the L1 method on graded and uniform meshes for a fractional-derivative problem in two and three dimensions, Math. Comput., 88 (2019), pp. 2135–2155, https://doi.org/10.1090/mcom/3410.
  * [13] B. Li, T. Wang, and X. Xie, Analysis of the L1 scheme for fractional wave equations with nonsmooth data, arXiv: 1908.09145v2 [math.NA].
  * [14] B. Li, T. Wang, and X. Xie, Analysis of a time-stepping discontinuous Galerkin method for fractional diffusion-wave equations with nonsmooth data, J. Sci. Comput., 82 (2020), https://doi.org/10.1007/s10915-019-01118-7.
  * [15] H. L. Liao, D. Li, and J. Zhang, Sharp error estimate of a nonuniform L1 formula for time-fractional reaction-subdiffusion equations, SIAM J. Numer. Anal., 56 (2018), pp. 1112–1133, https://doi.org/10.1137/17M1131829.
  * [16] H. L. Liao, W. McLean, and J. Zhang, A second-order scheme with nonuniform time steps for a linear reaction-subdiffusion problem, arXiv:1803.09873v2 [math.NA].
  * [17] H. L. Liao, W. McLean, and J. Zhang, A discrete Grönwall inequality with applications to numerical schemes for subdiffusion problems, SIAM J. Numer. Anal., 57 (2019), pp. 218–237, https://doi.org/10.1137/16M1175742.
  * [18] H. L. Liao and Z. Z. Sun, Maximum norm error bounds of ADI and compact ADI methods for solving parabolic equations, Numer. Meth. Part Differ. Equ., 26 (2010), pp. 37–60, https://doi.org/10.1002/num.20414.
  * [19] H. L. Liao, T. Tang, and T. Zhou, A second-order and nonuniform time-stepping maximum-principle preserving scheme for time-fractional Allen-Cahn equations, J. Comput. Phys., 414 (2020), https://doi.org/10.1016/j.jcp.2020.109473.
  * [20] H. L. Liao, Y. Yan, and J. Zhang, Unconditional convergence of a two-level linearized fast algorithm for semilinear subdiffusion equations, J. Sci. Comput., 80 (2019), pp. 1–25, https://doi.org/10.1007/s10915-019-00927-0.
  * [21] Y. Lin and C. Xu, Finite difference/spectral approximations for the time-fractional diffusion equation, J. Comput. Phys., 225 (2007), pp. 1533–1552, https://doi.org/10.1016/j.jcp.2007.02.001.
  * [22] H. Luo, B. Li, and X. Xie, Convergence analysis of a Petrov-Galerkin method for fractional wave problems with nonsmooth data, J. Sci. Comput., 80 (2019), pp. 957–992, https://doi.org/10.1007/s10915-019-00962-x.
  * [23] P. Lyu, Y. Liang, and Z. Wang, A fast linearized finite difference method for the nonlinear multi-term time-fractional wave equation, Appl. Numer. Math., 151 (2020), pp. 448–471, https://doi.org/10.1016/j.apnum.2019.11.012.
  * [24] F. Mainardi, Fractional Calculus and Waves in Linear Viscoelasticity, Imperial College Press, London, 2010.
  * [25] F. Mainardi and P. Paradisi, Fractional diffusive waves, J. Comput. Acoust., 9 (2001), pp. 1417–1436, https://doi.org/10.1142/S0218396X01000826.
  * [26] W. McLean and K. Mustapha, A second-order accurate numerical method for a fractional wave equation, Numer. Math., 105 (2007), pp. 481–510, https://doi.org/10.1007/s00211-006-0045-y.
  * [27] W. McLean and V. Thom$\acute{\mbox{e}}$e, Time discretization of an evolution equation via Laplace transforms, IMA J. Numer. Anal., 24 (2004), pp. 439–463, https://doi.org/10.1093/imanum/24.3.439.
  * [28] W. McLean and V. Thom$\acute{\mbox{e}}$e, Maximum-norm error analysis of a numerical solution via Laplace transformation and quadrature of a fractional-order evolution equation, IMA J. Numer. Anal., 30 (2010), pp. 208–230, https://doi.org/10.1093/imanum/drp004.
  * [29] K. Mustapha and W. McLean, Superconvergence of a discontinuous Galerkin method for fractional diffusion and wave equations, SIAM J. Nmer. Anal., 51 (2013), pp. 491–515, https://doi.org/10.1137/120880719.
  * [30] 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), pp. 1426–1446, https://doi.org/10.1093/imanum/drt048.
  * [31] R. R. Nigmatullin, To the theoretical explanation of the “Universal Response”, Phys. Status Solidi B, 123 (1984), pp. 739–745, https://doi.org/10.1002/pssb.2221230241.
  * [32] K. Oldham and J. Spanier, The Fractional Calculus, Academic Press, New York, London, 1974.
  * [33] I. Podnubny, Fractional Differential Equations, Academic Press, San Diego, London, 1999.
  * [34] 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), pp. 426–447, https://doi.org/10.1016/j.jmaa.2011.04.058.
  * [35] 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), pp. 1057–1079, https://doi.org/10.1137/16M1082329.
  * [36] Z. Z. Sun and X. N. Wu, A fully discrete difference scheme for a diffusion-wave system, Appl. Numer. Math., 56 (2006), pp. 193–209, https://doi.org/10.1016/j.apnum.2005.03.003.
  * [37] T. Tang, H. Yu, and T. Zhou, On energy dissipation theory and numerical stability for time-fractional phase-field equations, SIAM J. Sci. Comput., 41 (2019), pp. A3757–A3778, https://doi.org/10.1137/18M1203560.
  * [38] Y. Yan, M. Khan, and N. J. Ford, An analysis of the modified L1 scheme for time-fractional partial differential equations with nonsmooth data, SIAM J. Numer. Anal., 56 (2018), pp. 210–227, https://doi.org/10.1137/16M1094257.
  * [39] W. Zhang, J. Li, and Y. Yang, A fractional diffusion-wave equation with non-local regularization for image denoising, Signal Processing, 103 (2014), pp. 6–15, https://doi.org/10.1016/j.sigpro.2013.10.028.
  * [40] V. A. Zorich, Mathematical Analysis I, Springer, Berlin, 2004.