Error Analysis of the Nonuniform Alikhanov Scheme for the Fourth-Order Fractional Diffusion-Wave Equation

Abstract

This paper considers the numerical approximation to the fourth-order fractional diffusion-wave equation. Using a separation of variables, we can construct the exact solution for such a problem and then analyze its regularity. The obtained regularity result indicates that the solution behaves as a weak singularity at the initial time. Using the order reduction method, the fourth-order fractional diffusion-wave equation can be rewritten as a coupled system of low order, which is approximated by the nonuniform Alikhanov scheme in time and the finite difference method in space. Furthermore, the $H^2$-norm stability result is obtained. With the help of this result and a priori bounds of the solution, an $\alpha$-robust error estimate with optimal convergence order is derived. In order to further verify the accuracy of our theoretical analysis, some numerical results are provided.

1. Introduction

As of the last few decades, fractional partial differential equations (FPDEs) have gained widespread acceptance in both physics and engineering [1,2] to describe a variety of phenomena such as the universal response of electromagnetic fields [3], the denoising of images [4], etc. The exact solution to FPDEs can be constructed using tools such as Laplace transforms, Fourier transforms, and Mellin transforms; see [5,6]. In reality, only a few FPDEs can obtain the exact solution. Thus, the numerical simulation of FPDEs has been greatly developed. A number of numerical methods have recently been developed for solving time-fractional diffusion equations with a weakly singular solution. The regularity result for this problem was presented by Stynes et al. [7], and the L1 scheme was constructed on graded meshes. For the nonlinear time-fractional diffusion equation, Liao et al. [8] proposed the nonuniform Alikhanov scheme and constructed the corresponding discrete fractional Gronwall inequality. A nonuniform L2 scheme was designed by Kopteka [9] and the local error of the formulated scheme was estimated. By using the Müntz polynomial, Hou et al. [10] proposed a method for solving the time-fractional diffusion equation using the fractional spectral method.

The purpose of this paper is to develop a numerical method that is efficient in solving the following fourth-order fractional diffusion-wave equation:

$$\begin{array}{cc}\hfill D_t^{\alpha }u+{\kappa }^2{\Delta }^{2}u-c\Delta u& =g(t,\mathbf{x})\forall (t,\mathbf{x})\in \Phi :=(0,T]\times \Omega ,\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill u(0,\mathbf{x})& =\phi \left(\mathbf{x}\right)\mathrm{and}{u}_t(0,\mathbf{x})=\tilde{\phi }\left(\mathbf{x}\right)\mathrm{for}\mathbf{x}\in \Omega ,\hfill \end{array}$$

(1b)

$$\begin{array}{cc}\hfill {u|}_{\partial \Omega }& {=\Delta u|}_{\partial \Omega }=0\mathrm{for}0\le t\le T,\hfill \end{array}$$

(1c)

where $1<\alpha <2$, $\kappa >0$, $c\ge 0$, $\Omega =[x_l,x_r]\times [y_l,y_r]\subset {\mathbb{R}}^2$, $\phi \left(\mathbf{x}\right)\in C^{\infty }\left(\overline{\Omega }\right)$, and $\tilde{\phi }\left(\mathbf{x}\right)\in C^{\infty }\left(\overline{\Omega }\right)$. Here, $D_t^{\alpha }$ is the standard Caputo fractional derivative defined by

$$D_t^{\alpha }v(t,\mathbf{x})={\int }_0^t{\Lambda }_{2-\alpha }(t-s)\frac{{\partial }^{2}v(s,\mathbf{x})}{\partial s^2}ds\mathrm{with}{\Lambda }_{\rho }\left(t\right)=\frac{t^{\rho -1}}{\Gamma \left(\rho \right)}\mathrm{for}\rho >0.$$

(2)

As a matter of fact, Model (1) can be used to describe a wide variety of phenomena in physics, chemistry, engineering, and medicine, such as surfactant and liquid delivery into the lung [11], ice formation [12], forming grooves on flat surfaces [13], removing noise from medical magnetic resonance images [14], and so on. The main advantage of Model (1) is that it is able to be used to describe the evolution processes between diffusion and wave propagation, which is one of its most significant advantages. A number of numerical methods have recently been developed for solving the fractional diffusion-wave equation [15,16]. Despite this, it is relatively rare to simulate Equation (1) with a weakly singular solution numerically. A combination of the L1 scheme and the compact difference scheme is applied in the temporal and spatial directions by Hu and Zhang [17] to resolve Problem (1) by means of order reduction methods. As described in [18], Li and Vong solved Problem (1) by applying the parametric quintic spline in spatial direction and the weighted shifted Grünwald Letnikov scheme in the temporal direction. Bhrawy et al. [19] proposed an efficient and accurate spectral method to simulate Equation (1). To solve Equation (1) with nonlinear terms, Huang et al. [20] used piecewise rectangular formulas and the Crank–Nicolson technique in time, and unconditional convergence was derived using the energy method. A fast Alikhanov scheme was proposed by Ran and Lei [21] to approximate the Caputo derivative of Problem (1) with a variable coefficient, and the unconditional stability as well as the optimal convergence under the maximum norm of the scheme was demonstrated.

It should be noted that the above numerical methods are based on the assumption that the exact solution of Problem (1) is smooth. This assumption, as far as we know, is unrealistic, because at the initial time, the solution of time-fractional equations usually behaves as a weak singularity. Therefore, this paper establishes the regularity of u as well as the intermediate variable for Equation (1). By introducing two intermediate variables, $w=D_t^{\alpha /2}(u-t\tilde{\phi })$ and $p=\Delta (u-t\tilde{\phi })$, we are able to rewrite a fourth-order fractional diffusion-wave Equation (1) and approximate its solution using nonuniform Alikhanov schemes in temporal direction and the finite difference method in the spatial direction. Furthermore, we provide an optimal convergence analysis of the proposed method based on the $H^2$-norm which satisfies $\alpha$-robust.

We organize this paper as follows. The regularity of u and the introduced variable $w=D_t^{\alpha /2}(u-t\tilde{\phi })$ is analyzed in Section 2. In Section 3, an equivalent coupled system for Equation (1) is approximated using nonuniform Alikhanov schemes in time and the finite difference method in space. Moreover, Section 4 provides stability and convergence results in $H^2$-norm for the proposed scheme. In Section 5, two numerical experiments are presented to verify the accuracy of theoretical results. The final section of the paper presents our conclusions.

In the rest of the paper, we note that generic constant C is independent of the mesh, and it takes different values depending on where it is used. Furthermore, this constant is $\alpha$-robust, i.e., the constant’s value remains finite as $\alpha \to 2$.

2. Regularity of the Solution

In this section, the regularity of the exact solution for original Problem (1) is analyzed.

We denote $\tilde{u}(t,\mathbf{x}):=u(t,\mathbf{x})-t\tilde{\phi }\left(\mathbf{x}\right)$. Now, we reformulate the original Problem (1) as

$$\begin{array}{cc}\hfill D_t^{\alpha }\tilde{u}+{\kappa }^2{\Delta }^2\tilde{u}-c\Delta \tilde{u}& =g-t({\kappa }^2{\Delta }^2\tilde{\phi }-c\Delta \tilde{\phi })\forall (t,\mathbf{x})\in \Phi :=(0,T]\times \Omega ,\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill \tilde{u}(0,\mathbf{x})& =\phi \left(\mathbf{x}\right)\mathrm{and}{\tilde{u}}_t(0,\mathbf{x})=0\mathrm{for}\mathbf{x}\in \Omega ,\hfill \end{array}$$

(3b)

$$\begin{array}{cc}\hfill \tilde{u}{|}_{\partial \Omega }& =\Delta \tilde{u}{|}_{\partial \Omega }=0\mathrm{for}0\le t\le T.\hfill \end{array}$$

(3c)

Applying the separation of variables, we construct the solution to Problem (3). We let $\{({\lambda }_i,{\varphi }_i);i=1,2,\dots \}$ be the eigenvalues and eigenfunctions of the following Sturm–Liouville problem:

$$-\Delta {\varphi }_i={\lambda }_i{\varphi }_i\mathrm{with}{\varphi }_i{|}_{\partial \Omega }=0,$$

where ${\varphi }_i$ is the unit orthogonal basis of $L^2(\Omega )$ for all i and $0<{\lambda }_1<{\lambda }_2<\cdots <{\lambda }_i\to \infty$. From [22] (p. 3327), we observe that ${\varphi }_i$ is the eigenfunctions for the eigenvalue problem,

$${\kappa }^2{\Delta }^2{\varphi }_i-c\Delta {\varphi }_i={\tilde{\lambda }}_i{\varphi }_i\mathrm{with}{\varphi }_i{{|}_{\partial \Omega }=\Delta {\varphi }_i|}_{\partial \Omega }=0,$$

where eigenvalue ${\tilde{\lambda }}_i$ is defined by ${\tilde{\lambda }}_i:={\lambda }_i({\kappa }^2{\lambda }_i+c)$ for $i=1,2,\dots ,\infty$.

Now, applying the separation of variables to the equivalent Equation (3) yields

$$\tilde{u}(t,\mathbf{x})=\sum _{i=0}^{\infty }\left[(\phi ,{\varphi }_i)E_{\alpha ,1}(-{\tilde{\lambda }}_i{t}^{\alpha })+J_i\left(t\right)\right]{\varphi }_i\left(\mathbf{x}\right),$$

(4)

where $J_i\left(t\right)$ is defined by

$$J_i\left(t\right):={\int }_0^t{s}^{\alpha -1}{E}_{\alpha ,\alpha }(-{\tilde{\lambda }}_i{s}^{\alpha })g_i(t-s)ds\mathrm{with}{g}_i\left(t\right):=\left(g(t,·)-t({\kappa }^2{\Delta }^2\tilde{\phi }-c\Delta \tilde{\phi }),{\varphi }_i(·)\right).$$

Here, the classic Mittag–Leffler function $E_{\rho ,\eta }$ can be represented as follows:

$$E_{\rho ,\eta }\left(s\right)=\sum _{i=0}^{\infty }\frac{s^i}{\Gamma (\rho i+\eta )}\forall \rho >0,\eta >0.$$

Now, the following fractional power $L^{\gamma }$ is defined for each $\gamma \in \mathbb{R}$ by

$$\mathcal{D}\left(L^{\gamma }\right):=\left\{v\in L^2(\Omega ):\sum _{i=1}^{\infty }{\tilde{\lambda }}_i^{2\gamma }{(v,{\varphi }_i)}^2<\infty \right\},\mathrm{and}\mathrm{we}\mathrm{set}{\parallel v\parallel }_{L^{\gamma }}:={\left(\sum _{i=1}^{\infty }{\tilde{\lambda }}_i^{2\gamma }{(v,{\varphi }_i)}^2\right)}^{1/2}.$$

We set $w=D_t^{\alpha /2}(u-t\tilde{\phi })$. Imitating [23] (Section 3), a priori bounds on $\tilde{u}$ and w of Problem (1) are stated in the following lemma.

Lemma 1. 

We assume that $\phi \left(\mathbf{x}\right)\in \mathcal{D}\left(L^5\right)$, $\tilde{\phi }\left(\mathbf{x}\right)\in \mathcal{D}\left(L^5\right)$, $g(t,·)\in \mathcal{D}\left(L^5\right)$, $g_t(t,·)\in \mathcal{D}\left(L^3\right)$, and $g_{tt}(t,·)\in \mathcal{D}\left(L^3\right)$ for each $t\in (0,T]$ with

$${\parallel g(t,·)\parallel }_{L^5}+\parallel g_t{(t,·)\parallel }_{L^3}+{\parallel g_{tt}(t,·)\parallel }_{L^3}\le C.$$

Then, the solution of (3) and the introduced variable w satisfy

$$\begin{array}{cc}\hfill & {∥\tilde{u}(t,\mathbf{x})∥}_{H^6(\Omega )}+{∥w(t,\mathbf{x})∥}_{H^6(\Omega )}\le Cfort\in (0,T],\hfill \end{array}$$

(5a)

$$\begin{array}{cc}\hfill & {∥{\partial }_t^k\tilde{u}(t,\mathbf{x})∥}_{H^4(\Omega )}\le C(1+t^{\alpha -k})fort\in (0,T]andk=0,1,2,3,\hfill \end{array}$$

(5b)

$$\begin{array}{cc}\hfill & {∥{\partial }_t^{k}w(t,\mathbf{x})∥}_{H^4(\Omega )}\le C(1+t^{\alpha /2-k})fort\in (0,T]andk=0,1,2,3.\hfill \end{array}$$

(5c)

Proof. 

Applying the definition of the Caputo derivative and (4) yields

$$\begin{array}{cc}\hfill w(t,\mathbf{x})& =\frac{1}{\Gamma (1-\beta )}{\int }_0^t{(t-s)}^{-\beta }{\tilde{u}}_s(s,\mathbf{x})ds\hfill \\ \hfill & =\frac{1}{\Gamma (1-\beta )}\sum _{i=0}^{\infty }\left[-{\tilde{\lambda }}_i(\phi ,{\varphi }_i)S_1^i+g_i\left(0\right)S_1^i+S_2^i\right]{\varphi }_i\left(\mathbf{x}\right),\hfill \end{array}$$

(6)

where $S_1^i$ and $S_2^i$ are defined by

$$S_1^i:={\int }_0^t{(t-s)}^{-\beta }{s}^{\alpha -1}{E}_{\alpha ,\alpha }(-{\tilde{\lambda }}_i{s}^{\alpha })ds$$

and

$$S_2^i:={\int }_0^t{(t-s)}^{-\beta }\left[{\int }_0^s{\tau }^{\alpha -1}{E}_{\alpha ,\alpha }(-{\tilde{\lambda }}_i{\tau }^{\alpha })g_i^{\prime }(s-\tau )d\tau \right]ds.$$

Differentiating Series (4) and (6) with respect to t, x, and y, and then imitating the analysis of [23] (3.10) yields (5). □

3. The Fully Discrete Method

In this section, we construct the fully discrete method for Problem (1).

Now, we set $\beta :=\alpha /2$ in the remainder of the paper. From [15] (Lemma 2.1), the following significant property is stated to construct the numerical method.

Lemma 2 

([15] (Lemma 2.1)). We let $\tilde{v}\left(t\right)=v\left(t\right)-t{v}^{\prime }\left(0\right)$. We suppose function $v\left(t\right)\in C^2(0,T]$; we have

$$D_t^{\alpha }v\left(t\right)=D_t^{\alpha }\tilde{v}\left(t\right)=D_t^{\beta }\left(D_t^{\beta }\tilde{v}\left(t\right)\right)\mathrm{for}\alpha \in (1,2).$$

Using Lemma 2, Equation (3) can be rewritten as

$$\left\{\begin{array}{c}{D}_t^{\beta }w-{\kappa }^2\Delta p+cp=g(t,\mathbf{x})-t({\kappa }^2{\Delta }^2\tilde{\phi }-c\Delta \tilde{\phi })\forall (t,\mathbf{x})\in \Phi ,\hfill \\ D_t^{\beta }\tilde{u}-w=0\forall (t,\mathbf{x})\in \Phi ,\hfill \\ p+\Delta \tilde{u}=0\forall (t,\mathbf{x})\in \Phi ,\hfill \\ w(0,\mathbf{x})=0,p(0,\mathbf{x})=0,\mathrm{and}\tilde{u}(0,\mathbf{x})=\phi \left(\mathbf{x}\right)\mathrm{for}\mathbf{x}\in \Omega ,\hfill \\ \tilde{u}{|}_{\partial \Omega }=w{{|}_{\partial \Omega }=p|}_{\partial \Omega }=0\mathrm{for}0<t\le T.\hfill \end{array}\right.$$

(7)

For $n=0,1,\dots ,N$, we set $t_n=T{(n/N)}^r$, where N is a positive integer, and r represents the user’s choice of grading parameter. We denote $t_{n-\sigma }:=\sigma t_{n-1}+(1-\sigma )t_n$ for $n=1,2,\dots ,N$ and $0\le \sigma \le 1$.

For any function $v\in C[0,T]\cap C^3(0,T]$, the fractional derivative $D_t^{\beta }v(t,\mathbf{x})$ of (7) can be calculated at $t=t_{n-\sigma }$ using the following nonuniform Alikhanov method:

$$\begin{array}{cc}\hfill D_t^{\beta }v(\mathbf{x},t_{n-\sigma })\approx {\left(D_N^{\beta }v\right)}^{n-\sigma }& :={\int }_0^{t_{n-\sigma }}\frac{v^{\prime }\left(\eta \right)}{\Gamma (1-\beta ){(t_{n-\sigma }-\eta )}^{t_{n-\sigma }}}d\eta \hfill \\ \hfill & =\sum _{j=1}^{n-1}{\int }_{t_{j-1}}^{t_j}\frac{v^{\prime }\left(\eta \right)}{\Gamma (1-\beta ){(t_{n-\sigma }-\eta )}^{t_{n-\sigma }}}d\eta \hfill \\ \hfill & +{\int }_{t_{n-1}}^{t_{n-\sigma }}\frac{v^{\prime }\left(\eta \right)}{\Gamma (1-\beta ){(t_{n-\sigma }-\eta )}^{t_{n-\sigma }}}d\eta \hfill \\ \hfill & =\sum _{j=1}^{n-1}{\int }_{t_{j-1}}^{t_j}\frac{{\left(I_{2,j}v\left(\eta \right)\right)}^{\prime }}{\Gamma (1-\beta ){(t_{n-\sigma }-\eta )}^{t_{n-\sigma }}}d\eta \hfill \\ \hfill & +{\int }_{t_{n-1}}^{t_{n-\sigma }}\frac{{\delta }_{t}v\left(t_n\right)}{\Gamma (1-\beta ){(t_{n-\sigma }-\eta )}^{t_{n-\sigma }}}d\eta \hfill \\ \hfill & =\sum _{l=1}^n{d}_{n,n-l}^{\sigma }(v^l-v^{l-1}),\hfill \end{array}$$

(8)

in which $I_{2,j}v\left(\eta \right)$ is the quadratic polynomial interpolating at points $t_{j-1}$, $t_j$ and $t_{j+1}$. In (8), coefficients $d_{n,n-l}^{\sigma }$ are defined as follows: if $n=1$, $d_{1,0}^{\sigma }=a_0^1$; otherwise

$$d_{n,n-l}^{\sigma }=\left\{\begin{array}{cc}{a}_0^n+\frac{{\tau }_{n-1}}{{\tau }_n}{b}_1^n\hfill & \mathrm{if}l=n,\hfill \\ a_{n-l}^n+\frac{{\tau }_{l-1}}{{\tau }_k}{b}_{n-l+1}^n-b_{n-l}^n\hfill & \mathrm{if}2\le l\le n-1,\hfill \\ a_{n-1}^n-b_{n-1}^n\hfill & \mathrm{if}l=1,\hfill \end{array}\right.$$

where

$$\begin{array}{cc}\hfill & a_0^n=\frac{1}{{\tau }_n}{\int }_{t_{n-1}}^{t_{n-\sigma }}{\Lambda }_{1-\beta }(t_{n-\sigma }-\eta )d\eta ,a_{n-l}^n=\frac{1}{{\tau }_l}{\int }_{t_{l-1}}^{t_l}{\Lambda }_{1-\beta }(t_{n-\sigma }-\eta )d\eta \mathrm{for}1\le l\le n-1,\hfill \\ \hfill & b_{n-l}^n={\int }_{t_{l-1}}^{t_l}\frac{2(\eta -t_{l-\frac{1}{2}}){\Lambda }_{1-\beta }(t_{n-\sigma }-\eta )}{{\tau }_l({\tau }_l+{\tau }_{l+1})}d\eta \mathrm{for}1\le l\le n-1.\hfill \end{array}$$

Now, the temporal truncation error of the nonuniform Alikhanov scheme is stated in the following two lemmas.

Lemma 3 

([24] (Lemma 3.2)). We let $\sigma =1-\beta /2$. We suppose $\parallel {\partial }_t^l{v(t,\mathbf{x})\parallel }_{L^{\infty }(\Omega )}\le C(1+t^{\delta -l})$ for $l=0,1,2,3$; then, we obtain

$$t_{n-\sigma }^{\beta }\parallel {\left(D_N^{\beta }v\right)}^{n-\sigma }-D_t^{\beta }v(\mathbf{x},t_{n-\sigma }){\parallel }_{L^{\infty }(\Omega )}\le C{N}^{-min\{3-\beta ,r\delta \}}for1\le n\le N.$$

Lemma 4 

([24] (Lemma 3.3)). We let $\sigma =1-\beta /2$. We suppose $\parallel {\partial }_t^l{v(t,\mathbf{x})\parallel }_{L^{\infty }(\Omega )}\le C(1+t^{\delta -l})$ for $l=0,1,2$; then, we have

$${∥v^{n-\sigma }-v^{n,\sigma }∥}_{L^{\infty }(\Omega )}\le C{t}_{n-\sigma }^{-\beta }{N}^{-min\{2,r\delta \}}for1\le n\le N.$$

For spatial meshes, the uniform step sizes in each directions are defined as $h_x:=(x_r-x_l)/M_x$ and $h_y:=(y_r-y_l)/M_y$, where $M_x$ and $M_y$ are positive constants. We let ${\overline{\Omega }}_h=\left\{(x_i,y_j)\right|0\le i\le M_x,0\le j\le M_y\}$, ${\Omega }_h={\overline{\Omega }}_h\cap \Omega$, and $\partial {\Omega }_h={\overline{\Omega }}_h\cap \partial \Omega$, where $x_i:=x_l+i{h}_x$ and $y_j:=y_l+j{h}_y$. We denote $u_{i,j}^n$ or $u_h^n$ as the nodal approximation to $u(t_n,x_i,y_j)$ for all admissible $i,j$ and n.

For each grid function v on ${\overline{\Omega }}_h$, $\Delta v(t_n,x_i,y_j)$ can be calculated by

$$\Delta v(t_n,x_i,y_j)\approx {\Delta }_h{v}_h^n:={\delta }_x^2{v}_{i,j}^n+{\delta }_y^2{v}_{i,j}^n,$$

(9)

where

$${\delta }_x^2{v}_{i,j}^n:=\frac{v_{i+1,j}^n-2v_{i,j}^n+v_{i-1,j}^n}{h_x^2}\mathrm{and}{\delta }_y^2{u}_{i,j}^n:=\frac{v_{i,j+1}^n-2v_{i,j}^n+v_{i,j-1}^n}{h_y^2}.$$

At each point $(t_{n-\sigma },x_i,y_j)$, applying the nonuniform Alikhanov scheme (8) and Scheme (9) to approximate the Caputo derivative and the Laplacian operator of System (7) yields our fully discrete Alikhanov scheme:

$$\left\{\begin{array}{c}{\left(D_N^{\beta }{w}_h\right)}^{n-\sigma }-{\kappa }^2{\Delta }_h{p}_h^{n,\sigma }+c{p}_h^{n,\sigma }=g_h^{n-\sigma }-t_{n-\sigma }({\kappa }^2{\Delta }^2{\tilde{\phi }}_h^{n-\sigma }-c\Delta {\tilde{\phi }}_h^{n-\sigma }),\hfill \\ {\left(D_N^{\beta }{\tilde{u}}_h\right)}^{n-\sigma }-w_h^{n,\sigma }=0,\hfill \\ p_h^{n,\sigma }+{\Delta }_h{\tilde{u}}_h^{n,\sigma }=0\hfill \end{array}\right.$$

(10)

for $n=1,\dots ,N$, where ${\tilde{u}}_h^0=\phi \left(\mathbf{x}\right)$, $w_h^0=0$, and $p_h^0=0$.

4. The Error Estimate of the Fully Discrete Alikhanov Scheme

In this section, the $H^2$-norm stability result of the fully discrete Alikhanov scheme (10) is given. Moreover, we obtain an $\alpha$-robust convergent result.

For any grid function v and $\eta$ that vanish on $\partial \Omega$, we define the discrete inner product, $(v,\eta ):=h_x{h}_y{\sum }_{(x_i,y_j)\in {\Omega }_h}v(x_i,y_j)\eta (x_i,y_j)$, the discrete $L^2$-norm, $\parallel v\parallel :=\sqrt{(v,v)}$, and the discrete $H^1$-norm, $\parallel {\nabla }_{h}v\parallel :=\sqrt{(v,-{\Delta }_{h}v)}=\sqrt{({\nabla }_{h}v,{\nabla }_{h}v)}$. It is easy to check that

$$({\Delta }_h{v}_h,{\eta }_h)=(v_h,{\Delta }_h{\eta }_h)$$

(11)

for any $v_h$ and ${\eta }_h$ that vanish on $\partial \Omega$.

Now, an important property for the Alikhanov scheme (8) is stated.

Lemma 5 

([25] (Corollary 1)). For any grid function $v^n$, we have

$$2\left({\left(D_N^{\beta }v\right)}^{n-\sigma },v^{n,\sigma }\right)\ge {\left(D_N^{\beta }{\parallel v\parallel }^2\right)}^{n-\sigma }for1\le n\le N.$$

For $1\le n\le N$ and $1\le j\le n-1$, we investigate the complementary discrete convolution kernels, ${\Theta }_{n,n-j}$, which are defined by

$${\Theta }_{n,0}:=\frac{1}{d_{n,0}^{\sigma }},{\Theta }_{n,n-s}:=\frac{1}{d_{s,0}^{\sigma }}\sum _{k=s+1}^n\left(d_{k,k-s-1}^{\sigma }-d_{k,k-s}^{\sigma }\right){\Theta }_{n,n-k}.$$

From [26] (5.10), we obtain

$$\sum _{j=1}^n{\Theta }_{n,n-j}{j}^{r(\gamma -\beta )}\le \frac{{\pi }_A\Gamma (1+\gamma -\beta )}{\Gamma (1+\gamma )}{T}^{\beta }{\left(\frac{t_n}{T}\right)}^{\gamma }{N}^{r(\gamma -\beta )}\mathrm{for}\gamma \in (0,1)\mathrm{and}1\le n\le N,$$

(12)

where ${\pi }_A:=11/4$.

We now present the discrete fractional Gronwall inequality below.

Lemma 6 

([8] (Theorem 3.1)). We let the bounded sequences, $\{{\xi }^n\ge 0:n\ge 1\}$ and $\{v^n\ge 0:n\ge 0\}$, satisfy

$${\left(D_N^{\beta }{v}^2\right)}^{n-\sigma }\le {\xi }^n{v}^{n,\sigma }forn\ge 1;$$

(13)

then, we have

$$v^n\le v^0+\underset{1\le k\le n}{max}\sum _{j=1}^k{\Theta }_{k,k-j}{\xi }^{j}for1\le n\le N.$$

(14)

The next step is to consider the following general case, which includes (10) as a special case. For any $1\le n\le N$, we assume the pair of grid functions $({\omega }_h^n,{\nu }_h^n,{\mu }_h^n)$ that vanish on $\partial \Omega$ satisfy

$$\begin{array}{cc}\hfill {\left(D_N^{\beta }{\omega }_h\right)}^{n-\sigma }-{\kappa }^2{\Delta }_h{\nu }_h^{n,\sigma }+c{\nu }_h^{n,\sigma }& ={\zeta }_h^{n-\sigma }\mathrm{for}n=1,\dots ,N,\hfill \end{array}$$

(15a)

$$\begin{array}{cc}\hfill {\left(D_N^{\beta }{\mu }_h\right)}^{n-\sigma }-{\omega }_h^{n,\sigma }& =f_h^{n-\sigma }\mathrm{for}n=1,\dots ,N,\hfill \end{array}$$

(15b)

$$\begin{array}{cc}\hfill {\nu }_h^{n,\sigma }+{\Delta }_h{\mu }_h^{n,\sigma }& ={\varsigma }_h^{n-\sigma }\mathrm{for}n=1,\dots ,N,\hfill \end{array}$$

(15c)

where ${\omega }_h^n$, ${\mu }_h^n$, and ${\nu }_h^n$ vanish on $\partial \Omega$, and ${\zeta }_h^n$, $f_h^n$, and ${\varsigma }_h^n$ are given grid functions on ${\Omega }_h$.

Next, we present a useful bound of (15), which is applied to derive the stability result of the fully discrete Alikhanov scheme (10).

Lemma 7. 

The solution $({\omega }_h^n,{\mu }_h^n,{\nu }_h^n)$ of (15) satisfies

$$\begin{array}{cc}\hfill & \sqrt{\parallel {\omega }_h^n{\parallel }^2+{\kappa }^2\parallel {\Delta }_h{\mu }_h^n{\parallel }^2+c{\parallel {\nabla }_h{\mu }_h^n\parallel }^2}\le \sqrt{\parallel {\omega }_h^0{\parallel }^2+{\kappa }^2\parallel {\Delta }_h{\mu }_h^0{\parallel }^2+c{\parallel {\nabla }_h{\mu }_h^0\parallel }^2}\hfill \\ \hfill & +2\underset{1\le k\le n}{max}\left[\sum _{s=1}^k{\Theta }_{n,k-s}\sqrt{\parallel {\zeta }_h^{s-\sigma }+{\kappa }^2{\Delta }_h{\varsigma }_h^{s-\sigma }-c{\varsigma }_h^{s-\sigma }{\parallel }^2+{\parallel {\Delta }_h{f}_h^{s-\sigma }-\frac{c}{\kappa }{f}_h^{s-\sigma }\parallel }^2}\right]\hfill \end{array}$$

(16)

for $1\le n\le N$.

Proof. 

We fix $n\in \{1,2,\dots ,N\}$. Taking the inner product of Equations (15a)–(15c) by ${\omega }_h^{n,\sigma }$, ${\kappa }^2{\Delta }_h{\mu }_h^{n,\sigma }$, and ${\kappa }^2{\Delta }_h{\omega }_h^{n,\sigma }$ respectively, we have

$$\left\{\begin{array}{c}\left({\left(D_N^{\beta }{\omega }_h\right)}^{n-\sigma },{\omega }_h^{n,\sigma }\right)-{\kappa }^2({\Delta }_h{\nu }_h^{n,\sigma },{\omega }_h^{n,\sigma })+c({\nu }_h^{n,\sigma },{\omega }_h^{n,\sigma })=({\zeta }_h^{n-\sigma },{\omega }_h^{n,\sigma }),\hfill \\ {\kappa }^2\left({\left(D_N^{\beta }{\mu }_h\right)}^{n-\sigma },{\Delta }_h^2{\mu }_h^{n,\sigma }\right)-{\kappa }^2({\omega }_h^{n,\sigma },{\Delta }_h^2{\mu }_h^{n,\sigma })={\kappa }^2(f_h^{n-\sigma },{\Delta }_h^2{\mu }_h^{n,\sigma }),\hfill \\ {\kappa }^2({\nu }_h^{n,\sigma },{\Delta }_h{\omega }_h^{n,\sigma })+{\kappa }^2({\Delta }_h{\mu }_h^{n,\sigma },{\Delta }_h{\omega }_h^{n,\sigma })={\kappa }^2({\varsigma }_h^{n-\sigma },{\Delta }_h{\omega }_h^{n,\sigma }).\hfill \end{array}\right.$$

Adding the above three equations, we have

$$\begin{array}{cc}\hfill \left({\left(D_N^{\beta }{\omega }_h\right)}^{n-\sigma },{\omega }_h^{n,\sigma }\right)+{\kappa }^2\left({\left(D_N^{\beta }{\Delta }_h{\mu }_h\right)}^{n-\sigma },{\Delta }_h{\mu }_h^{n,\sigma }\right)& +c({\nu }_h^{n,\sigma },{\omega }_h^{n,\sigma })=({\zeta }_h^{n-\sigma }+{\kappa }^2{\Delta }_h{\varsigma }_h^{n-\sigma },{\omega }_h^{n,\sigma })\hfill \\ \hfill & +{\kappa }^2({\Delta }_h{f}_h^{n-\sigma },{\Delta }_h{\mu }_h^{n,\sigma }),\hfill \end{array}$$

(17)

where (11) is used. By taking the inner product of Equations (15b) and (15c) with ${\Delta }_h{\mu }_h^{n,\sigma }$ and ${\omega }_h^{n,\sigma }$, we have

$$\left\{\begin{array}{c}\left({\left(D_N^{\beta }{\mu }_h\right)}^{n-\sigma },{\Delta }_h{\mu }_h^{n,\sigma }\right)-({\omega }_h^{n,\sigma },{\Delta }_h{\mu }_h^{n,\sigma })=(f_h^{n-\sigma },{\Delta }_h{\mu }_h^{n,\sigma }),\hfill \\ ({\nu }_h^{n,\sigma },{\omega }_h^{n,\sigma })+({\Delta }_h{\mu }_h^{n,\sigma },{\omega }_h^{n,\sigma })=({\varsigma }_h^{n-\sigma },{\omega }_h^{n,\sigma }).\hfill \end{array}\right.$$

Adding the above two equations and applying (11), we have

$$({\nu }_h^{n,\sigma },{\omega }_h^{n,\sigma })=\left({\left(D_N^{\beta }{\nabla }_h{\mu }_h\right)}^{n-\sigma },{\nabla }_h{\mu }_h^{n,\sigma }\right)+(f_h^{n-\sigma },{\Delta }_h{\mu }_h^{n,\sigma })+({\varsigma }_h^{n-\sigma },{\omega }_h^{n,\sigma }).$$

Substituting the above result into (17) yields

$$\begin{array}{cc}\hfill \left({\left(D_N^{\beta }{\omega }_h\right)}^{n-\sigma },{\omega }_h^{n,\sigma }\right)& +{\kappa }^2\left({\left(D_N^{\beta }{\Delta }_h{\mu }_h\right)}^{n-\sigma },{\Delta }_h{\mu }_h^{n,\sigma }\right)+c\left({\left(D_N^{\beta }{\nabla }_h{\mu }_h\right)}^{n-\sigma },{\nabla }_h{\mu }_h^{n,\sigma }\right)\hfill \\ \hfill & =\left({\zeta }_h^{n-\sigma }+{\kappa }^2{\Delta }_h{\varsigma }_h^{n-\sigma }-c{\varsigma }_h^{n-\sigma },{\omega }_h^{n,\sigma }\right)+({\kappa }^2{\Delta }_h{f}_h^{n-\sigma }-c{f}_h^{n-\sigma },{\Delta }_h{\mu }_h^{n,\sigma }).\hfill \end{array}$$

(18)

Using a Cauchy–Schwarz inequality and Lemma 5, we have

$$\begin{array}{cc}\hfill \frac{1}{2}{\left(D_N^{\beta }{\parallel {\omega }_h\parallel }^2\right)}^{n-\sigma }& +\frac{1}{2}{\kappa }^2{\left(D_N^{\beta }{\parallel {\Delta }_h{\mu }_h\parallel }^2\right)}^{n-\sigma }+\frac{1}{2}c{\left(D_N^{\beta }{\parallel {\nabla }_h{\mu }_h\parallel }^2\right)}^{n-\sigma }\hfill \\ \hfill & \le \parallel {\zeta }_h^{n-\sigma }+{\kappa }^2{\Delta }_h{\varsigma }_h^{n-\sigma }-c{\varsigma }_h^{n-\sigma }\parallel \parallel {\omega }_h^{n,\sigma }\parallel +\parallel {\kappa }^2{\Delta }_h{f}_h^{n-\sigma }-c{f}_h^{n-\sigma }\parallel \parallel {\Delta }_h{\mu }_h^{n,\sigma }\parallel .\hfill \end{array}$$

Applying Hölder inequality produces

$$\begin{array}{cc}\hfill & {\left(D_N^{\beta }\left(\parallel {\omega }_h{\parallel }^2+{\kappa }^2\parallel {\Delta }_h{\mu }_h{\parallel }^2+c{\parallel {\nabla }_h{\mu }_h\parallel }^2\right)\right)}^{n-\sigma }\hfill \\ \hfill & \le 2\sqrt{\parallel {\zeta }_h^{n-\sigma }+{\kappa }^2{\Delta }_h{\varsigma }_h^{n-\sigma }-c{\varsigma }_h^{n-\sigma }{\parallel }^2+{\parallel {\Delta }_h{f}_h^{n-\sigma }-\frac{c}{\kappa }{f}_h^{n-\sigma }\parallel }^2}\hfill \\ \hfill & \times \sqrt{\parallel {\omega }_h^{n,\sigma }{\parallel }^2+{\kappa }^2\parallel {\Delta }_h{\mu }_h^{n,\sigma }{\parallel }^2+c{\parallel {\nabla }_h{\mu }_h^{n,\sigma }\parallel }^2}.\hfill \end{array}$$

Furthermore, applying Lemma 6 yields

$$\begin{array}{cc}\hfill & \sqrt{\parallel {\omega }_h^n{\parallel }^2+{\kappa }^2\parallel {\Delta }_h{\mu }_h^n{\parallel }^2+c{\parallel {\nabla }_h{\mu }_h^n\parallel }^2}\le \sqrt{\parallel {\omega }_h^0{\parallel }^2+{\kappa }^2\parallel {\Delta }_h{\mu }_h^0{\parallel }^2+c{\parallel {\nabla }_h{\mu }_h^0\parallel }^2}\hfill \\ \hfill & +2\underset{1\le k\le n}{max}\left[\sum _{s=1}^k{\Theta }_{n,k-s}\sqrt{\parallel {\zeta }_h^{s-\sigma }+{\kappa }_h^2{\Delta }_h{\varsigma }_h^{s-\sigma }-c{\varsigma }_h^{s-\sigma }{\parallel }^2+{\parallel {\Delta }_h{f}_h^{s-\sigma }-\frac{c}{\kappa }{f}_h^{s-\sigma }\parallel }^2}\right].\hfill \end{array}$$

□

We now present the stability result for the fully discrete Alikhanov scheme (10) in the following theorem.

Theorem 1. 

Solution $(w_h^n,{\tilde{u}}_h^n,p_h^n)$ of (10) satisfies

$$\begin{array}{cc}\hfill \sqrt{\parallel w_h^n{\parallel }^2+{\kappa }^2\parallel {\Delta }_h{\tilde{u}}_h^n{\parallel }^2+c{\parallel {\nabla }_h{\tilde{u}}_h^n\parallel }^2}\le & \sqrt{\parallel w_h^0{\parallel }^2+{\kappa }^2\parallel {\Delta }_h{\tilde{u}}_h^0{\parallel }^2+c{\parallel {\nabla }_h{\tilde{u}}_h^0\parallel }^2}\hfill \\ \hfill & +\frac{2{\pi }_A}{\Gamma (1+\beta )}{t}_n^{\beta }\underset{1\le s\le n}{max}\parallel g_h^{s-\sigma }-t_{s-\sigma }({\kappa }^2{\Delta }^2{\tilde{\phi }}_h-c\Delta {\tilde{\phi }}_h)\parallel \hfill \end{array}$$

(19)

for $1\le n\le N$.

Proof. 

It is clear that (10) is a special case of (15) with ${\omega }_h^n=w_h^n$, ${\nu }_h^n=p_h^n$, ${\mu }_h^n={\tilde{u}}_h^n$, ${\zeta }_h^{n-\sigma }=g_h^{n-\sigma }-t_{n-\sigma }({\kappa }^2{\Delta }^2{\tilde{\phi }}_h-c\Delta {\tilde{\phi }}_h)$, $f_h^{n-\sigma }=0$, and ${\varsigma }_h^{n-\sigma }=0$. Hence, applying Lemma 7 yields

$$\begin{array}{cc}\hfill \sqrt{\parallel w_h^n{\parallel }^2+{\kappa }^2\parallel {\Delta }_h{\tilde{u}}_h^n{\parallel }^2+c{\parallel {\nabla }_h{\tilde{u}}_h^n\parallel }^2}& \le \sqrt{\parallel w_h^0{\parallel }^2+{\kappa }^2\parallel {\Delta }_h{\tilde{u}}_h^0{\parallel }^2+c{\parallel {\nabla }_h{\tilde{u}}_h^0\parallel }^2}\hfill \\ \hfill & +2\underset{1\le k\le n}{max}\sum _{s=1}^k{\theta }_{n,k-s}\parallel g_h^{s-\sigma }-t_{s-\sigma }({\kappa }^2{\Delta }^2{\tilde{\phi }}_h-c\Delta {\tilde{\phi }}_h)\parallel .\hfill \end{array}$$

(20)

From (20), we obtain (19) immediately by invoking (12) with $\gamma =\beta$. □

Next, we state the error analysis of our fully discrete scheme. First, we denote

$$e_w^n:=w_h^n-w(t_n,x_i,y_j),e_u^n:={\tilde{u}}_h^n-\tilde{u}(t_n,x_i,y_j),\mathrm{and}{e}_p^n:=p_h^n-p(t_n,x_i,y_j).$$

Now, we begin to present the error equation of coupled system (10). For $1len\le N$, we define

$$\begin{array}{cc}\hfill & {\phi }_w^{n-\sigma }:=D_t^{\beta }w(·,t_{n-\sigma })-{\left(D_N^{\beta }w\right)}^{n-\sigma },R_p^{n-\sigma }:=p^{n-\sigma }-p^{n,\sigma },\hfill \\ \hfill & {\phi }_u^{n-\sigma }:=D_t^{\beta }u(·,t_{n-\sigma })-{\left(D_N^{\beta }u\right)}^{n-\sigma },R_w^{n-\sigma }:=w^{n-\sigma }-w^{n,\sigma },\hfill \\ \hfill & {\rho }_p^{n,\sigma }:=\Delta p^{n,\sigma }-{\Delta }_h{p}^{n,\sigma },{\rho }_u^{n,\sigma }:=\Delta {\tilde{u}}^{n,\sigma }-{\Delta }_h{u}^{n,\sigma }.\hfill \end{array}$$

From (7) and (10), we obtain the following error equations:

$$\begin{array}{cc}\hfill {\left(D_N^{\beta }{e}_w\right)}^{n-\sigma }-{\kappa }^2{\Delta }_h{e}_p^{n,\sigma }+c{e}_p^{n,\sigma }& ={\chi }^{n-\sigma },\hfill \end{array}$$

(21a)

$$\begin{array}{cc}\hfill {\left(D_N^{\beta }{e}_u\right)}^{n-\sigma }-e_w^{n,\sigma }& ={\phi }_u^{n-\sigma }-R_w^{n-\sigma },\hfill \end{array}$$

(21b)

$$\begin{array}{cc}\hfill e_p^{n,\sigma }+{\Delta }_h{e}_u^{n,\sigma }& ={\rho }_u^{n,\sigma },\hfill \end{array}$$

(21c)

where ${\chi }^{n-\sigma }$ is defined by

$${\chi }^{n-\sigma }:={\phi }_w^{n-\sigma }-{\kappa }^2\Delta R_p^{n-\sigma }-{\kappa }^2{\rho }_p^{n-\sigma }+c{R}_p^{n-\sigma }.$$

The optimal error estimate for Scheme (10) is given by combining the regularity results given in Section 2.

Theorem 2. 

We let $l_N=1/lnN$. At each time level $t_n$, we let $(w^n,{\tilde{u}}^n,p^n)$ and $(w_h^n,{\tilde{u}}_h^n,p_h^n)$ be the solutions of (7) and (10), respectively. Then, we have

$$\parallel {\Delta }_h{\tilde{u}}^n-{\Delta }_h{\tilde{u}}_h^n\parallel +\parallel w^n-w_h^n\parallel \le \frac{C{\pi }_A{t}_n^{\beta }}{\Gamma (1+\beta )}\left(h_x^2+h_y^2\right)+\frac{C{\pi }_A{e}^r\Gamma (1+l_N-\beta )}{\Gamma (1+\beta )}{N}^{-min\{2,r\beta \}}$$

(22)

for $1\le n\le N$.

Proof. 

It is clear that System of error equations (21) is a special case of (15) with ${\omega }_h^n=e_w^n$, ${\nu }_h^n=e_p^n$, ${\mu }_h^n=e_u^n$, ${\zeta }_h^{n-\sigma }={\chi }^{n-\sigma }$, $f_h^{n-\sigma }={\phi }_u^{n-\sigma }-R_w^{n-\sigma }$, and ${\varsigma }_h^{n-\sigma }={\rho }_u^{n,\sigma }$. Hence, invoking Lemma 7 offers

$$\begin{array}{cc}\hfill & \sqrt{\parallel e_w^n{\parallel }^2+{\kappa }^2\parallel {\Delta }_h{e}_u^n{\parallel }^2+c{\parallel {\nabla }_h{e}_u^n\parallel }^2}\le \sqrt{\parallel e_w^0{\parallel }^2+{\kappa }^2\parallel {\Delta }_h{e}_u^0{\parallel }^2+c{\parallel {\nabla }_h{e}_u^0\parallel }^2}\hfill \\ \hfill & +2\underset{1\le k\le n}{max}\sum _{s=1}^k{\Theta }_{k,k-s}\sqrt{\parallel {\chi }^{s-\sigma }+{\kappa }^2{\Delta }_h{\rho }_u^{s,\sigma }-c{\rho }_u^{s,\sigma }{\parallel }^2+{\parallel {\Delta }_h({\phi }_u^{s-\sigma }-R_w^{s-\sigma })-\frac{c}{\kappa }({\phi }_u^{s-\sigma }-R_w^{s-\sigma })\parallel }^2}\hfill \\ \hfill & \le 2\underset{1\le k\le n}{max}\sum _{s=1}^k{\Theta }_{k,k-s}\left(\parallel {\chi }^{s-\sigma }\parallel +(c{C}_{\Omega }+\kappa )\parallel {\Delta }_h{\rho }_u^{s,\sigma }\parallel +(1+\frac{c{C}_{\Omega }}{\kappa })(\parallel {\Delta }_h{\phi }_u^{s-\sigma }\parallel +\parallel {\Delta }_h{R}_w^{s-\sigma }\parallel )\right),\hfill \end{array}$$

(23)

where $e_u^0=0$, $e_w^0=0$, and $\parallel v\parallel \le C_{\Omega }\parallel {\Delta }_{h}v\parallel$ are used.

Applying the triangle inequality and regularity Result (5), we obtain

$$\begin{array}{cc}\hfill \parallel {\chi }^{s-\sigma }\parallel & \le \parallel {\phi }_w^{s-\sigma }\parallel +{\kappa }^2\parallel \Delta R_p^{s-\sigma }\parallel +{\kappa }^2\parallel {\rho }_p^{s-\sigma }\parallel +c\parallel R_p^{s-\sigma }\parallel \hfill \\ \hfill & \le C{t}_{s-\sigma }^{-\beta }{N}^{-min\{3-\beta ,r\beta \}}+C{t}_{s-\sigma }^{-\beta }{N}^{-min\{2,r\alpha \}}+C(h_x^2+h_y^2)\hfill \\ \hfill & \le C\left(h_x^2+h_y^2+t_{s-\sigma }^{-\beta }{N}^{-min\{2,r\beta \}}\right),\hfill \end{array}$$

(24)

where Lemma 3 and Lemma 4 are used. Applying the Taylor expansion, we arrive at

$$\begin{array}{cc}\hfill {\Delta }_h{\rho }_u^{s,\sigma }=& {\int }_0^1\left[{\partial }_{xx}{\rho }_u^{s,\sigma }(x_i-l{h}_x,y_j)+{\partial }_{xx}{\rho }_u^{s,\sigma }(x_i+l{h}_x,y_j)\right](1-l)dl\hfill \\ \hfill & +{\int }_0^1\left[{\partial }_{yy}{\rho }_u^{s,\sigma }(x_i,y_j-l{h}_y)+{\partial }_{yy}{\rho }_u^{s,\sigma }(x_i,y_j-l{h}_y)\right](1-l)dl.\hfill \end{array}$$

Hence, combining the above Taylor expansion with regularity Result (5a) yields

$$\parallel {\Delta }_h{\rho }_u^{s,\sigma }\parallel \le C\left(h_x^2+h_y^2\right).$$

(25)

Similarity to the estimation of $\parallel {\Delta }_h{\rho }_u^{s,\sigma }\parallel$, we have

$$\begin{array}{cc}\hfill \parallel {\Delta }_h{\phi }_u^{s-\sigma }\parallel +\parallel {\Delta }_h{R}_w^{s-\sigma }\parallel & \le C{t}_{s-\sigma }^{-\beta }{N}^{-min\{3-\beta ,r\beta \}}+C{t}_{s-\sigma }^{-\beta }{N}^{-min\{2,r\beta \}}\hfill \\ \hfill & \le C{t}_{s-\sigma }^{-\beta }{N}^{-min\{2,r\beta \}}.\hfill \end{array}$$

(26)

Substituting (24)–(26) into (23) yields

$$\sqrt{\parallel e_w^n{\parallel }^2+{\kappa }^2\parallel {\Delta }_h{e}_u^n{\parallel }^2+c{\parallel {\nabla }_h{e}_u^n\parallel }^2}\le 2C\underset{1\le k\le n}{max}\sum _{s=1}^k{\Theta }_{k,k-s}\left(h_x^2+h_y^2+t_{s-\sigma }^{-\beta }{N}^{-min\{2,r\beta \}}\right).$$

(27)

However,

$$t_{s-\sigma }^{-\beta }\le t_{s-1}^{-\beta }=T^{-\beta }{(s-1)}^{-r\beta }{N}^{r\beta }\le T^{-\beta }{\left(\frac{s}{2}\right)}^{-r\beta }{N}^{r\beta }\le 2^{r\beta }{T}^{-\beta }{N}^{r\beta }{s}^{r(l_N-\beta )}.$$

Utilizing the above bound and using (12) with $\gamma =l_N$ offers

$$\sum _{s=1}^k{\Theta }_{k,k-s}{t}_{s-\sigma }^{-\beta }{N}^{-min\{2,r\beta \}}\le \frac{C{\pi }_A{e}^r\Gamma (1+l_N-\beta )}{\Gamma (1+\beta )}{N}^{-min\{2,r\beta \}}.$$

(28)

Invoking (12) with $\gamma =\beta$ for the term $h_x^2+h_y^2$, we have

$$\sum _{s=1}^k{\Theta }_{k,k-s}\left(h_x^2+h_y^2\right)\le \frac{{\pi }_A}{\Gamma (1+\beta )}{t}_k^{\beta }\left(h_x^2+h_y^2\right).$$

(29)

Substituting (28) and (29) into (27) offers

$$\sqrt{\parallel e_w^n{\parallel }^2+{\kappa }^2\parallel {\Delta }_h{e}_u^n{\parallel }^2+c{\parallel {\nabla }_h{e}_u^n\parallel }^2}\le \frac{C{\pi }_A{t}_n^{\beta }}{\Gamma (1+\beta )}\left(h_x^2+h_y^2\right)+\frac{C{\pi }_A{e}^r\Gamma (1+l_N-\beta )}{\Gamma (1+\beta )}{N}^{-min\{2,r\beta \}}.$$

Finally, the desired Result (22) is obtained immediately. □

Remark 1. 

In [23], the Caputo fractional derivative of Problem (1) is decomposed into $D_t^{\alpha -1}v$ and $u_t$ by investigating the variable $v=u_t$, then they are approximated by the nonuniform L1 scheme and the nonuniform Crank–Nicolson scheme respectively. Furthermore, the proposed numerical method achieves the $min\{r\alpha /2,3-\alpha \}$ order in time, which is suboptimal. However, our convergent result given in Theorem 2 attains the 2 order in time. In addition, Theorem 2 shows that our error analysis is α-robust as $\alpha \to 2^{-}$. Nevertheless, the convergent result given in [15] contains the term $\Gamma (1-\alpha /2)$, which blows up as $\alpha \to 2^{-}$.

Corollary 1. 

We suppose $r\ge 4/\alpha$; then, we obtain that numerical solution $({\tilde{u}}_h^n,w_h^n)$ satisfies

$$\parallel \Delta {\tilde{u}}^n-{\Delta }_h{\tilde{u}}_h^n\parallel +\parallel w^n-w_h^n\parallel \le \frac{C{\pi }_A{t}_n^{\alpha /2}}{\Gamma (1+\alpha /2)}\left(h_x^2+h_y^2\right)+\frac{C{\pi }_A{e}^r\Gamma (1+l_N-\alpha /2)}{\Gamma (1+\alpha /2)}{N}^{-2}$$

for $1\le n\le N$.

5. Numerical Experiments

In this section, two numerical experiments with weakly singular solutions are presented to verify the theoretical result of our fully discrete Alikhanov scheme (10).

Example 1. 

We consider Equation (1) with $\kappa =1$, $c=1$, $\Omega =(0,\pi )\times (0,\pi )$, $T=1$, $\phi (x,y)=0$, $\tilde{\phi }(x,y)=0$, and the source term g is chosen by

$$g(t,x,y)=\left[\Gamma (1+\alpha )+\frac{6}{\Gamma (4-\alpha )}+6(t^{\alpha }+t^3)\right]sinxsiny.$$

The theoretical solution of this numerical example is $u(t,x,y)=(t^{\alpha }+t^3)sinxsiny$, which displays weak singularity at $t=0$.

In our following calculation, we estimate the global $H^2$-norm error ${max}_{1\le n\le N}\parallel \Delta u^n-{\Delta }_h{u}_h^n\parallel$ of the Alikhanov scheme by taking $r=4/\alpha$.

Table 1. Convergent results of $\underset{1\le j\le N}{max}\parallel \Delta u^n-{\Delta }_h{u}_h^n\parallel$ in temporal direction with $r=4/\alpha$.

	N = 20	N = 40	N = 80	N = 160	N = 320
$\alpha =1.2$	3.4823$·{10}^{-2}$	9.1635$·{10}^{-3}$	2.3618$·{10}^{-3}$	5.9906$·{10}^{-4}$	1.5053$·{10}^{-4}$
		1.9261	1.9959	1.9791	1.9925
$\alpha =1.5$	3.2409$·{10}^{-2}$	8.3604$·{10}^{-3}$	2.1279$·{10}^{-3}$	5.3611$·{10}^{-4}$	1.3425$·{10}^{-4}$
		1.9547	1.9740	1.9888	1.9976
$\alpha =1.8$	3.1744$·{10}^{-2}$	8.1407$·{10}^{-3}$	2.0645$·{10}^{-3}$	5.1916$·{10}^{-4}$	1.2969$·{10}^{-4}$
		1.9632	1.9793	1.9915	2.0010

shows the $H^2$-norm error and the rate of the Alikhanov scheme for $\alpha =1.2,1.5,1.8$, where $N=M_x=M_y$ is taken. The results displayed indicate that temporal convergence rates are $N^{-2}$, in agreement with Theorem 2.

Table 2. Convergent results of $\underset{1\le j\le N}{max}\parallel \Delta u^n-{\Delta }_h{u}_h^n\parallel$ in spatial direction.

	${\mathit{M}}_{\mathit{x}}={\mathit{M}}_{\mathit{y}}=10$	${\mathit{M}}_{\mathit{x}}={\mathit{M}}_{\mathit{y}}=20$	${\mathit{M}}_{\mathit{x}}={\mathit{M}}_{\mathit{y}}=40$	${\mathit{M}}_{\mathit{x}}={\mathit{M}}_{\mathit{y}}=80$
$\alpha =1.2$	9.9170$·{10}^{-3}$	2.4773$·{10}^{-3}$	6.1757$·{10}^{-4}$	1.5528$·{10}^{-4}$
		2.0014	2.0040	1.9917
$\alpha =1.5$	9.9080$·{10}^{-3}$	2.4771$·{10}^{-3}$	6.1970$·{10}^{-4}$	1.5699$·{10}^{-4}$
		1.9999	1.9990	1.9836
$\alpha =1.8$	1.0077$·{10}^{-2}$	2.5205$·{10}^{-3}$	6.3258$·{10}^{-4}$	1.6172$·{10}^{-4}$
		1.9993	1.9943	1.9677

shows the errors and their associated spatial orders for different $\alpha$, where $N=200$ and $M_x=M_y$ are taken. Form Table 2, we observe $O\left(h^2\right)$ convergence, again as predicted by Theorem 2.

Example 2. 

Let us solve the original Problem (1) with $\kappa =c=1$, $\Omega =(0,\pi )\times (0,\pi )$, $T=1$, $\phi (x,y)=sinxsiny$, $\tilde{\phi }(x,y)=0$, and $g=0$.

Due to the fact that the solution of this numerical example is unknown, we use the two mesh principles given in [27] to verify the theoretical result of Theorem 2. Figure 1 displays the error $\underset{1\le n\le N}{max}\parallel \Delta u^n-{\Delta }_h{u}_h^n\parallel$ in time for different $\alpha$, where $N=M_x=M_y$ is chosen. It is shown that the rate of convergence attains the 2 order in time, as predicted by Theorem 2. Figure 2 presents the numerical solution at $t_n=1$ for $\alpha =1.1,1.9$. The displayed result of this figure shows that the solution of this problem behaves diffusive as $\alpha \to 1^{+}$, and the characteristic of the initial solution $\phi (x,y)$ is well maintained as $\alpha \to 2^{-}$, which indicates that the propagation speed of waves is finite. Moreover, the attenuation of the numerical solution gradually slows down as $\alpha \to 2^{-}$, verifying the wave feature again.

6. Conclusions

This paper investigates the regularity result of Equation (1), which indicates that the solution behaves as a weak singularity at $t=0$. In order to construct a fully discrete scheme, the nonuniform Alikhanov scheme in time is used along with the finite difference method in space. Then, the stability analysis and an $\alpha$-robust optimal convergent result under $H^2$-norm for the proposed scheme are obtained. Finally, two numerical examples are provided to verify the theoretical result. In the future, we will extend the proposed scheme to the Equation (1) with a nonlinear term and present the unconditional optimal error estimate. In addition, we will extend the proposed scheme to fractional equations with the Caputo-Hadamard fractional derivative.