Compact Difference Method for Time-Fractional Neutral Delay Nonlinear Fourth-Order Equation

Huan Wang, Qing Yang

School of Mathematics and Statistics, Shandong Normal University, Jinan, China.

DOI: 10.4236/eng.2022.1412041                

Abstract

In this paper, we present a compact finite difference method for a class of fourth-order nonlinear neutral delay sub-diffusion equations in two-dimensional space. The fourth-order problem is first transformed into a second-order system by a reduced-order method. Next by using compact operator to approximate the second order space derivatives and L2-1σ formula to approximate the time fractional derivative, the difference scheme which is fourth order in space and second order in time is obtained. Then, the existence and uniqueness of solution, the convergence rate of and the stability of the scheme are proved. Finally, numerical results are given to verify the accuracy and validity of the scheme.

Keywords

Two-Dimensional Nonlinear Sub-Diffusion Equations, Neutral Delay, Compact Difference Scheme, Convergence, Stability

1. Introduction

It has been found that in many experiments and researches the diffusion process of many complex systems no longer satisfies Fick’s second law. Such process is called anomalous diffusion, one remarkable feature of which is that the mean square displacement of the particle and the time variable have the following power-law dependence:

$\langle x^2\left(t\right)\rangle ~\frac{2K_{\alpha }}{\Gamma \left(1+\alpha \right)}{t}^{\alpha },t\to \infty$, (1.1)

where α is the anomalous diffusion exponent and K_α is the generalized diffusion coefficient. When $0<\alpha <1$, we have sub-diffusion, and when $1<\alpha <2$, we have super-diffusion. Various table text styles are provided. The formatter will need to create these components, incorporating the applicable criteria that follow. Anomalous diffusion process can be described by differential equations involving fractional calculus, that is, fractional order diffusion equations.

In recent years, the application of fractional diffusion equations in mechanics, physics viscoelastic mechanics [1] [2] [3] [4], porous media [5] [6], hydrology [7] and so on has been an important research topic. There are several different methods to solve the analytical solutions of some fractional diffusion equations, such as Laplace method, Fourier method, Green function method, etc. [8] [9] [10] [11]. But only certain types of fractional differential equation can be solved for exact solutions. Therefore, in most cases we need to rely on numerical methods. In the past few years, a number of different numerical methods for fractional diffusion equations have been developed [12] - [19]. In [20] [21] [22] [23], L1 formula is introduced for discretization of fractional diffusion equations, and the precision in time is $O\left({\tau }^{2-\alpha }\right)$. In order to improve the precision in time, Zhao and Sun [24] obtained a second order approximation of time fractional derivative by using Crank-Nicolson method. In [25], Gao and Sun proposed the L1-2 approximation formula. On the basis of L1-2 formula, Alikhanov [26] established the L2-1σ approximation formula, which has the 3 − α order uniform convergence rate.

In many applications, it is necessary to use equations which contain fourth-order derivatives in space, for example, wave propagation of light beam [27], modeling of Plane Grooves [28], the formation of ice [29] [30] and the propagation of intense laser beams through the Quer [31] [32] body. In [33], Agrawa gave the analytic solution for fourth order fractional diffusion-wave equation by means of Laplace and Fourier transform. In [34], Hu and Zhang studied finite difference method for spatial fourth-order fractional diffusion equations, and the convergence order of the scheme is $O\left({\tau }^{2-\alpha }+h^2\right)$. Later, in [35], Guo and Li et al. proposed two numerical schemes for the equations in literature [34], and proved unconditional stability and the convergence order of $O\left({\tau }^2+h^2\right)$ of the two schemes. In [36], Sun and Ji constructed a difference scheme of fourth order fractional diffusion equations with the first Dirichlet boundary condition, and obtained fourth order spatial accuracy. Then in [37], Zhang et al. constructed a compact finite difference scheme for fourth order fractional diffusion equations with the second Dirichlet boundary condition by using the L2-1σ formula to approximate the time fractional derivative and the compact operator to approximate spatial fourth order derivative.

Delay differential equations are widely used in many fields, such as population ecology, cell biology, control theory, economics and so on [38] [39] [40]. In [38] [39] Sarita Ndal et al. discuss finite difference scheme for one-dimensional time fractional fourth-order diffusion equation, which contains a nonlinear source function with time delay and a fourth-order space delay term. The unique solvability, stability and convergence of the scheme are proved.

In this work, we consider the following fourth-order nonlinear sub-diffusion neutral delayed equation in two-dimensional space.

$\begin{array}{l}{}_0{}^{c}D{}_t^{\alpha }u\left(x,y,t\right)+{\Delta }^{2}u\left(x,y,t-s\right)\\ =f\left(x,y,t,u\left(x,y,t\right),u\left(x,y,t-s\right)\right),\left(x,y,t\right)\in \Omega \times \left(0,T\right).\end{array}$ (1.2)

subject to initial and boundary conditions:

$\{\begin{array}{l}u\left(x,y,t\right)=\phi \left(x,y,t\right),\left(x,y,t\right)\in \partial \Omega \times \left[0,T\right],\\ \Delta u\left(x,y,t\right)=\psi \left(x,y,t\right),\left(x,y,t\right)\in \partial \Omega \times \left[0,T\right],\\ u\left(x,y,t\right)=\varphi \left(x,y,t\right),\left(x,y,t\right)\in \Omega \times \left[-s,0\right],\end{array}$ (1.3)

where $s>0$ is the delay, $f\left(x,y,t,u\left(x,y,t\right),u\left(x,y,t-s\right)\right)$ represents a nonlinear source term with time delay, $\Omega =\left(0,L_1\right)\times \left(0,L_2\right)$, $\alpha \in \left(0,1\right)$, $f,\varphi ,\psi ,\phi$ are all given and sufficiently smooth functions.

The fractional derivative is defined in Caputo form:

${}_0{}^{c}D{}_t^{\alpha }u\left(x,y,t\right)\equiv \frac{{\partial }^{\alpha }u\left(x,y,t\right)}{\partial t^{\alpha }}:=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^t{\left(t-\xi \right)}^{-\alpha }}\frac{\partial u\left(x,y,\xi \right)}{\partial \xi }\text{d}\xi ,0<\alpha <1$. (1.4)

Let m be the integer satisfying $ms\le T\le \left(m+1\right)s$. Define $I_r=\left(rs,\left(r+1\right)s\right)$, $r=-1,0,\cdots ,m-1$, $I_m=\left(ms,T\right)$, $I=U_{q=-1}{I}_q$.

Assume that the partial derivatives $f_{\mu }\left(x,y,t,\mu ,\nu \right)$ and $f_{\nu }\left(x,y,t,\mu ,\nu \right)$ are continuous in the ${ϵ}_0$ -neighborhood of $\mu$ and $\nu$ for a positive constant ${ϵ}_0$. Define

$c_1={\sup }_{\begin{array}{c}0<x<L_1,0<y<L_2,0<t<T\\ \left|{ϵ}_1\right|\le {ϵ}_0,\left|{ϵ}_2\right|\le {ϵ}_0\end{array}}\left|f_{\mu }\left(x,y,t,u\left(x,y,t\right)+{ϵ}_1,u\left(x,y,t-s\right)+{ϵ}_2\right)\right|,$ (1.5)

$c_2={\sup }_{\begin{array}{c}0<x<L_1,0<y<L_2,0<t<T\\ \left|{ϵ}_1\right|\le {ϵ}_0,\left|{ϵ}_2\right|\le {ϵ}_0\end{array}}\left|f_{\nu }\left(x,y,t,u\left(x,y,t\right)+{ϵ}_1,u\left(x,y,t-s\right)+{ϵ}_2\right)\right|.$ (1.6)

In this article, we construct a compact difference scheme for the problem (1.2)-(1.3). First, the fourth-order equation is transformed into two second-order equations by introducing $v=\Delta u$ as an auxiliary variable. Next, we use the L2-1σ formula and the compact difference operator to approximate temporal Caputo derivative and spatial derivative. Doing in this way, we can obtain the compact difference scheme. Then the existence and uniqueness of the difference scheme is proved. By using mathematical induction and energy method, we prove the convergence and stability of the scheme. Numerical experiments show that the scheme can achieve convergence order $O\left({\tau }^2+h^4\right)$ in discrete $L^2$ norm and $L^{\infty }$ norm, which verify the theoretical results.

2. Notations and Preliminary

Introducing positive integers $M_1,M_2$ and n, let $h_x=\frac{L_1}{M_1}$, $h_y=\frac{L_2}{M_2}$ and $\tau =\frac{s}{n}$, be the spatial steps in x and y directions and the temporal step respectively. Define $x_i=i{h}_x\left(0\le i\le M_1\right)$, $y_j=j{h}_y\left(0\le j\le M_2\right)$, $t_k=k\tau \left(-n\le k\le N\right)$, where $N=\left[\frac{T}{t}\right]$. Define ${\bar{\Omega }}_h=\left\{\left(x_i,y_j\right)|0\le i\le M_1,0\le j\le M_2\right\}$, ${\Omega }_h={\bar{\Omega }}_h\cap \Omega$, $\partial {\Omega }_h={\bar{\Omega }}_h\cap \partial \Omega$, $\omega =\left\{\left(i,j\right)|\left(x_i,y_j\right)\in {\Omega }_h\right\}$, $\partial \omega =\left\{\left(i,j\right)|\left(x_i,y_j\right)\in \partial {\Omega }_h\right\}$, ${\Omega }_t=\left\{t_k|-n\le k\le N\right\}$ .

In addition, denote $t_{k-1+\sigma }=\left(k-1+\sigma \right)\tau$, where $\sigma =\frac{\alpha }{2}$. Define the grid function space $S_h=\left\{u|u=\left\{u_{i,j}\in R\left(M_1+1\right)\times \left(M_2+1\right),0\le i\le M_1,0\le j\le M_2\right\}\right\}$ and $S_{h0}=\left\{u|u\in S_h,u_{0,j}=u_{M_1,j},0\le j\le M_2;u_{i,0}=u_{i,M_2},0\le i\le M_1\right\}$ on ${\Omega }_h$.

Define the following difference operators for any grid functions $u\in S_h$,

$\begin{array}{l}{\delta }_x{u}_{i-\frac{1}{2},j}=\frac{1}{h_x}\left(u_{i,j}-u_{i-1,j}\right),{\delta }_y{u}_{i,j-\frac{1}{2}}=\frac{1}{h_y}\left(u_{i,j}-u_{i,j-1}\right),\\ {\delta }_x^2{u}_{i,j}=\frac{1}{h_x^2}\left(u_{i+1,j}-2u_{i,j}+u_{i-1,j}\right),{\delta }_y^2{u}_{i,j}=\frac{1}{h_y^2}\left(u_{i,j+1}-2u_{i,j}+u_{i,j-1}\right).\end{array}$ (2.1)

Definition 2.1. For $u\in S_h$, the compact difference operator is defined as follows

$A_x{u}_{i,j}=\{\begin{array}{ll}\left(I+\frac{h_x^2}{12}{\delta }_x^2\right)u_{i,j},\hfill & 1\le i\le M_1-1,0\le j\le M_2,\hfill \\ u_{i,j},\hfill & i=0\text{or}i=M_1,0\le j\le M_2,\hfill \end{array}$ (2.2)

$A_y{u}_{i,j}=\{\begin{array}{ll}\left(I+\frac{h_y^2}{12}{\delta }_y^2\right)u_{i,j},\hfill & 1\le j\le M_2-1,0\le j\le M_1,\hfill \\ u_{i,j},\hfill & j=0\text{or}j=M_2,0\le j\le M_1,\hfill \end{array}$ (2.3)

Denote

$A_h=A_x{A}_y,B_h=A_y{\delta }_x^2+A_x{\delta }_y^2.$ (2.4)

Then define the following discrete inner products and the corresponding norms for $u,v\in S_h$,

$\left(u,v\right)=h_x{h}_y{\displaystyle {\sum }_{i=1}^{M_1-1}{\displaystyle {\sum }_{j=1}^{M_2-1}{u}_{i,j}{v}_{i,j}}},{\Vert u\Vert }^2=\left(u,u\right)$, (2.5)

${\left({\delta }_{x}u,{\delta }_{x}v\right)}_x=h_x{h}_y{\displaystyle {\sum }_{i=1}^{M_1-1}{\displaystyle {\sum }_{j=1}^{M_2-1}{\delta }_x{u}_{i-\frac{1}{2},j}{\delta }_x{v}_{i,j}}},{\Vert {\delta }_{x}u\Vert }_x^2={\left({\delta }_{x}u,{\delta }_{x}u\right)}_x$, (2.6)

${\left({\delta }_x^{2}u,{\delta }_x^{2}v\right)}_x=h_x{h}_y{\displaystyle {\sum }_{i=1}^{M_1-1}{\displaystyle {\sum }_{j=1}^{M_2-1}{\delta }_x^2{u}_{i,j}{\delta }_x^2{v}_{i,j}}},{\Vert {\delta }_x^{2}u\Vert }_x^2={\left({\delta }_x^{2}u,{\delta }_x^{2}u\right)}_x$, (2.7)

${\left({\delta }_{y}u,{\delta }_{y}v\right)}_y,{\left({\delta }_y^{2}u,{\delta }_y^{2}v\right)}_y,{\Vert {\delta }_{y}u\Vert }_y^2,{\Vert {\delta }_y^{2}u\Vert }_y^2$ can be defined similarly. We also using the following norm

${\Vert u\Vert }_{\infty }={\max }_{\begin{array}{l}1\le i\le M_1-1\\ {}_{1\le j\le M_2-1}\end{array}}\left|u_{i,j}\right|.$ (2.8)

Some related lemmas for difference operators ${\delta }_x^2,{\delta }_y^2,A_x,A_y,A_h,B_h$ are given as follows.

Lemma 2.1. [41] Denote $\zeta \left(\lambda \right)={\left(1-\lambda \right)}^3\left[5-3{\left(1-\lambda \right)}^2\right]$ and let $g\left(•,y\right)\in C^6\left[x_i-h_x,x_i+h_x\right]$, $g\left(x,·\right)\in C^6\left[y_j-h_y,y_j+h_y\right]$ for any function $g\left(x,y\right)$. We have

$A_x{g}_{xx}\left(x_i,y_j\right)={\delta }_x^{2}g\left(x_i,y_j\right)+\frac{h_x^4}{360}{\displaystyle {\int }_0^1\left[g_x^{\left(6\right)}\left(x_i-\lambda h,y_j\right)+g_x^{\left(6\right)}\left(x_i+\lambda h,y_j\right)\right]\zeta \left(\lambda \right)\text{d}\lambda },$ (2.9)

$A_y{g}_{yy}\left(x_i,y_j\right)={\delta }_y^{2}g\left(x_i,y_j\right)+\frac{h_y^4}{360}{\displaystyle {\int }_0^1\left[g_y^{\left(6\right)}\left(x_i,y_j-\lambda h\right)+g_y^{\left(6\right)}\left(x_i,y_j+\lambda h\right)\right]\zeta \left(\lambda \right)\text{d}\lambda }.$ (2.10)

Lemma 2.2. [42] For any grid functions $u,v\in S_{h0}$, we have

$\left(A_{x}u,v\right)=\left(u,A_{x}v\right),\left(A_{y}u,v\right)=\left(u,A_{y}v\right);$ (2.11)

$\left(A_{h}u,v\right)=\left(u,A_{h}v\right),\left(A_{h}u,v\right)=\left(u,A_{h}v\right);$ (2.12)

$\left({\delta }_x^{2}u,v\right)=\left(u,{\delta }_x^{2}v\right),\left({\delta }_y^{2}u,v\right)=\left(u,{\delta }_y^{2}v\right);$ (2.13)

$\left(A_{h}u,B_{h}v\right)=\left(B_{h}u,A_{h}v\right).$ (2.14)

Proof. Because the first six equalities have been proven by Q. Li in [42], we only prove the last equality. According to the definition of operator $A_h,B_h$, and applying the first and third equalities above, we get

$\begin{array}{c}\left(A_{h}u,B_{h}v\right)=\left(A_x{A}_{y}u,A_y{\delta }_x^{2}v+A_x{\delta }_y^{2}v\right)\\ =\left(A_x{A}_{y}u,A_y{\delta }_x^{2}v\right)+\left(A_x{A}_{y}u,A_x{\delta }_y^{2}v\right)\\ =\left(A_x{A}_y{\delta }_x^{2}u,A_{y}v\right)+\left(A_x{A}_y{\delta }_y^{2}u,A_{x}v\right)\\ =\left(A_y{\delta }_x^{2}u,A_x{A}_{y}v\right)+\left(A_x{\delta }_y^{2}u,A_x{A}_{y}v\right)\\ =\left(B_{h}u,A_{h}v\right).\end{array}$ (2.15)

The lemma has been proved. ◻

Lemma 2.3. [41] For any grid function $u\in S_h$, we have

$\begin{array}{l}\frac{1}{3}{\Vert u\Vert }^2\le {\Vert A_{x}u\Vert }^2\le {\Vert u\Vert }^2,\frac{1}{3}{\Vert u\Vert }^2\le {\Vert A_{y}u\Vert }^2\le {\Vert u\Vert }^2,\\ \frac{1}{9}{\Vert u\Vert }^2\le {\Vert A_{h}u\Vert }^2\le {\Vert u\Vert }^2.\end{array}$ (2.16)

Lemma 2.4. [42] For any grid function $u\in S_h$ , we have

$\begin{array}{l}\frac{2}{3}{\Vert u\Vert }^2\le \left(A_{x}u,u\right)\le {\Vert u\Vert }^2,\frac{2}{3}{\Vert u\Vert }^2\le \left(A_{y}u,u\right)\le {\Vert u\Vert }^2,\\ \frac{4}{9}{\Vert u\Vert }^2\le \left(A_{h}u,u\right)\le {\Vert u\Vert }^2.\end{array}$ (2.17)

According to [26], the L2-1σ approximation formula is defined as

${\Delta }_{\tau }^{\alpha }{u}_{i,j}^{k-1+\sigma }=\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left[c_0^{\left(k\right)}{u}_{i,j}^k-{\displaystyle {\sum }_{m=1}^{k-1}\left(c_{k-m-1}^{\left(k\right)}-c_{k-m}^{\left(k\right)}\right)u_{i,j}^m}-c_{k-1}^{\left(k\right)}{u}_{i,j}^0\right]$. (2.18)

where,

$a_0={\sigma }^{1-\alpha },a_l={\left(l+\sigma \right)}^{1-\alpha }-{\left(l-1+\sigma \right)}^{1-\alpha },l\ge 1$,

$b_l=\frac{1}{2-\alpha }\left[{\left(l+\sigma \right)}^{2-\alpha }-{\left(l-1+\sigma \right)}^{2-\alpha }\right]-\frac{1}{2}\left[{\left(l+\sigma \right)}^{1-\alpha }+{\left(l-1+\sigma \right)}^{1-\alpha }\right],l\ge 1$.(2.19)

when $k=1,c_0^{\left(k\right)}=a_0$, and when $k\ge 2$,

$c_m^{\left(k\right)}=\{\begin{array}{l}{a}_0+b_1,m=0,\\ a_m+b_{m+1}-b_m,1\le m\le k-2,\\ a_m-b_m,m=k-1.\end{array}$ (2.20)

Lemma 2.5. [26] For $a\left(t\right)\in C^3\left[0,t_k\right]$, we have the following error estimates

$\left|{{}_0{}^{c}D{}_t^{\alpha }a\left(t\right)|}_{t=t_{k-1+\sigma }}-{{\Delta }_{\tau }^{\alpha }a\left(t\right)|}_{t=t_{k-1+\sigma }}\right|\le \frac{\left(4\sigma -1\right){\sigma }^{-\alpha }}{12\Gamma \left(2-\alpha \right)}{\max }_{0\le t\le t_k}\left|a^{\left(3\right)}\left(t\right)\right|{\tau }^{3-\alpha }.$ (2.21)

Lemma 2.6. [26] Let $w=\left\{w^k|-n\le k\le N\right\}$ be a grid function defined on ${\Omega }_{\tau }$, then

$\left(\sigma w^k+\left(1-\sigma \right)w^{k-1}\right){}_0{}^{c}D{}_{t_{k-1+\sigma }}^{\alpha }w\ge \frac{1}{2}{}_0{}^{c}D{}_{t_{k-1+\sigma }}^{\alpha }\left(w^2\right).$ (2.22)

Lemma 2.7. [26] For $\alpha \in \left(0,1\right)$, $\sigma =1-\frac{\alpha }{2}$ and $c_m^{\left(k\right)}\left(0\le m\le k-1,k\ge 1\right)$ defined by (2.4), we have the following inequalities

$c_m^{\left(k\right)}>\frac{1-\alpha }{2}{\left(k+\sigma \right)}^{-\alpha }>0,$ (2.23)

$c_0^{\left(k\right)}>c_1^{\left(k\right)}>c_2^{\left(k\right)}>\cdots >c_{k-2}^{\left(k\right)}>c_{k-1}^{\left(k\right)}.$ (2.24)

3. Construction of the Compact Difference Scheme

In this section, we will construct the compact difference scheme of the problem (1.2)-(1.3). Let $v\left(x,y,t\right)=\Delta u\left(x,y,t\right)$. Then the problem (1.2)-(1.3) is equivalent to

$\begin{array}{l}{}_0{}^{c}D{}_t^{\alpha }u\left(x,y,t\right)+\Delta v\left(x,y,t\right)+\Delta v\left(x,y,t-s\right)\\ =f\left(x,y,t,u\left(x,y,t\right),u\left(x,y,t-s\right)\right),\left(x,y,t\right)\in \Omega \times \left[0,T\right],\end{array}$ (3.1)

$v\left(x,y,t\right)=\Delta u\left(x,y,t\right),\left(x,y,t\right)\in \Omega \times \left[0,T\right],$ (3.2)

$v\left(x,y,t-s\right)=\Delta u\left(x,y,t-s\right),\left(x,y,t\right)\in \Omega \times \left[0,s\right],$ (3.3)

$u\left(x,y,t\right)=\phi \left(x,y,t\right),\left(x,y,t\right)\in \partial \Omega \times \left[0,T\right],$ (3.4)

$v\left(x,y,t\right)=\psi \left(x,y,t\right),\left(x,y,t\right)\in \partial \Omega \times \left[0,T\right],$ (3.5)

$u\left(x,y,t\right)=\varphi \left(x,y,t\right),\left(x,y,t\right)\in \Omega \times \left[-s,0\right].$ (3.6)

In addition, we define the following grid functions for the exact solutions $u$ and $v$

$U_{i,j}^k=u\left(x_i,y_j,t_k\right),V_{i,j}^k=v\left(x_i,y_j,t_k\right),0\le i\le M_1,0\le j\le M_2,-n\le k\le N.$ (3.7)

Considering the Equations (3.1)-(3.3) at the grid point $\left(x_i,y_j,t_{k-1+\sigma }\right)$ and applying operators $A_h$, we have

$\begin{array}{l}{A}_h\frac{{\partial }^{\alpha }u\left(x_i,y_j,t_{k-1+\sigma }\right)}{\partial t^{\alpha }}+A_h\frac{{\partial }^{2}v\left(x_i,y_j,t_{k-1+\sigma }\right)}{\partial x^2}+A_h\frac{{\partial }^{2}v\left(x_i,y_j,t_{k-1+\sigma }\right)}{\partial y^2}\\ +A_h\frac{{\partial }^{2}v\left(x_i,y_j,t_{k-1+\sigma -n}\right)}{\partial x^2}+A_h\frac{{\partial }^{2}v\left(x_i,y_j,t_{k-1+\sigma -n}\right)}{\partial y^2}\\ =A_{h}f\left(x_i,y_j,t_{k-1+\sigma },u\left(x_i,y_j,t_{k-1+\sigma }\right),u\left(x_i,y_j,t_{k-1+\sigma -n}\right)\right),\end{array}$ (3.8)

$A_{h}v\left(x_i,y_j,t_{k-1+\sigma }\right)=A_h\frac{{\partial }^{2}u\left(x_i,y_j,t_{k-1+\sigma }\right)}{\partial x^2}+A_h\frac{{\partial }^{2}u\left(x_i,y_j,t_{k-1+\sigma }\right)}{\partial y^2},$ (3.9)

$A_{h}v\left(x_i,y_j,t_{k-1+\sigma -n}\right)=A_h\frac{{\partial }^{2}u\left(x_i,y_j,t_{k-1+\sigma -n}\right)}{\partial x^2}+A_h\frac{{\partial }^{2}u\left(x_i,y_j,t_{k-1+\sigma -n}\right)}{\partial y^2}.$ (3.10)

From Taylor’s series expansion, we can get

$U_{i,j}^{k-1+\sigma }=u\left(x_i,y_j,t_{k-1+\sigma }\right)=\left(\sigma +1\right)U_{i,j}^{k-1}-\sigma \left(\sigma +1\right)U_{i,j}^{k-2}+O\left({\tau }^2\right),$ (3.11)

$U_{i,j}^{k-1+\sigma -n}=u\left(x_i,y_j,t_{k-1+\sigma -n}\right)=\sigma U_{i,j}^{k-n}+\left(1-\sigma \right)U_{i,j}^{k-1-n}+O\left({\tau }^2\right).$ (3.12)

Then the nonlinear source term $f\left(x,y,t,u\left(x,y,t\right),u\left(x,y,t-s\right)\right)$ can be approximated by the following formula

$\begin{array}{l}f\left(x_i,y_j,t_{k-1+\sigma },u\left(x_i,y_j,t_{k-1+\sigma }\right),u\left(x_i,y_j,t_{k-1+\sigma -n}\right)\right)\\ =f\left(x_i,y_j,t_{k-1+\sigma },\left(\sigma +1\right)U_{i,j}^{k-1}-\sigma U_{i,j}^{k-2},\sigma U_{i,j}^{k-n}+\left(1-\sigma \right)U_{i,j}^{k-n-1}\right)+O\left({\tau }^2\right).\end{array}$ (3.13)

By using Lemma 2.1, we have

$\begin{array}{l}{A}_h\frac{{\partial }^{2}u\left(x_i,y_j,t_{k-1+\sigma }\right)}{\partial x^2}\\ =\sigma A_h\frac{{\partial }^{2}u\left(x_i,y_j,t_k\right)}{\partial x^2}+\left(1-\sigma \right)A_h\frac{{\partial }^{2}u\left(x_i,y_j,t_{k-1}\right)}{\partial x^2}+O\left({\tau }^2+h_x^4\right)\\ =A_y{\delta }_x^2{U}_{i,j}^{k-1+\sigma }+O\left({\tau }^2+h_x^4\right),\end{array}$ (3.14)

$A_h\frac{{\partial }^{2}u\left(x_i,y_j,t_{k-1+\sigma }\right)}{\partial y^2}=A_x{\delta }_y^2{U}_{i,j}^{k-1+\sigma }+O\left({\tau }^2+h_y^4\right),$ (3.15)

$\begin{array}{l}{A}_h\frac{{\partial }^{2}v\left(x_i,y_j,t_{k-1+\sigma }\right)}{\partial x^2}+A_h\frac{{\partial }^{2}v\left(x_i,y_j,t_{k-1+\sigma }\right)}{\partial y^2}\\ =A_y{\delta }_x^2{V}_{i,j}^{k-1+\sigma }+A_x{\delta }_y^2{V}_{i,j}^{k-1+\sigma }+O\left({\tau }^2+h_x^4+h_y^4\right)\\ =B_h{V}_{i,j}^{k-1+\sigma }+O\left({\tau }^2+h_x^4+h_y^4\right),\end{array}$ (3.16)

$\begin{array}{l}{A}_h\frac{{\partial }^{2}v\left(x_i,y_j,t_{k-1+\sigma -n}\right)}{\partial x^2}+A_h\frac{{\partial }^{2}v\left(x_i,y_j,t_{k-1+\sigma -n}\right)}{\partial y^2}\\ =A_y{\delta }_x^2{V}_{i,j}^{k-1+\sigma -n}+A_x{\delta }_y^2{V}_{i,j}^{k-1+\sigma -n}+O\left({\tau }^2+h_x^4+h_y^4\right)\\ =B_h{V}_{i,j}^{k-1+\sigma -n}+O\left({\tau }^2+h_x^4+h_y^4\right).\end{array}$ (3.17)

Substituting Equations (3.13)-(3.17) into (3.8)-(3.10), and approximating the time fractional derivative by L2-1σ formula (2.18), we can get

$\begin{array}{l}\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left[c_0^{\left(k\right)}{A}_h{U}_{i,j}^k-\underset{m=1}{\overset{k-1}{{\displaystyle \sum }}}\left(c_{k-m-1}^{\left(k\right)}-c_{k-m}^{\left(k\right)}\right)A_h{U}_{i,j}^m-c_{k-1}^{\left(k\right)}{A}_h{U}_{i,j}^0\right]\\ +B_h{V}_{i,j}^{k-1+\sigma }+B_h{V}_{i,j}^{k-1+\sigma -n}\\ =A_{h}f\left(x_i,y_j,t_{k-1+\sigma },\left(\sigma +1\right)U_{i,j}^{k-1}-\sigma U_{i,j}^{k-2},\sigma U_{i,j}^{k-n}+\left(1-\sigma \right)U_{i,j}^{k-n-1}\right)+{\left(R_1\right)}_{i,j}^k,\end{array}$ (3.18)

$A_h{V}_{i,j}^{k-1+\sigma }=B_h{U}_{i,j}^{k-1+\sigma }+{\left(R_2\right)}_{i,j}^k,$ (3.19)

$A_h{V}_{i,j}^{k-1+\sigma -n}=B_h{U}_{i,j}^{k-1+\sigma -n}+{\left(R_3\right)}_{i,j}^k,$ (3.20)

where

$\left|{\left(R_1\right)}_{i,j}^k\right|+\left|{\left(R_2\right)}_{i,j}^k\right|+\left|{\left(R_3\right)}_{i,j}^k\right|\le \hat{c}\left({\tau }^2+h_x^4+h_y^4\right).$ (3.21)

where $\hat{c}$ is a positive constant that does not depend on τ andh.

Omitting $R_1^k,R_2^k$ and $R_3^k$ and in Equations (3.18)-(3.20), and using numerical solution $u_{i,j}^k,v_{i,j}^k$ to replace the exact solution $U_{i,j}^k,V_{i,j}^k$, the following difference scheme for problem (1.2)-(1.3) is obtained

$\begin{array}{l}\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left[c_0^{\left(k\right)}{A}_h{u}_{i,j}^k-\underset{m=1}{\overset{k-1}{{\displaystyle \sum }}}\left(c_{k-m-1}^{\left(k\right)}-c_{k-m}^{\left(k\right)}\right)A_h{u}_{i,j}^m-c_{k-1}^{\left(k\right)}{A}_h{u}_{i,j}^0\right]\\ +B_h{v}_{i,j}^{k-1+\sigma }+B_h{v}_{i,j}^{k-1+\sigma -n}\\ =A_{h}f\left(x_i,y_j,t_{k-1+\sigma },\left(\sigma +1\right)u_{i,j}^{k-1}-\sigma u_{i,j}^{k-2},\sigma u_{i,j}^{k-n}+\left(1-\sigma \right)u_{i,j}^{k-n-1}\right),\\ 1\le i\le M_1-1,1\le j\le M_2-1,1\le k\le N-1,\end{array}$ (3.22)

$1\le i\le M_1-1,1\le j\le M_2-1,1\le k\le N-1,$

$A_h{v}_{i,j}^{k-1+\sigma }=B_h{u}_{i,j}^{k-1+\sigma },$ (3.23)

$A_h{v}_{i,j}^{k-1+\sigma -n}=B_h{u}_{i,j}^{k-1+\sigma -n},$ (3.24)

with the following discrete initial and boundary conditions

$u_{i,j}^k=\phi \left(x_i,y_j,t_k\right),\left(i,j\right)\in \omega ,-n\le k\le 0,$ (3.25)

$v_{i,j}^k=\psi \left(x_i,y_j,t_k\right),\left(i,j\right)\in \partial \omega ,1\le k\le N,$ (3.26)

$u_{i,j}^k=\varphi \left(x_i,y_j,t_k\right),\left(i,j\right)\in \partial \omega ,1\le k\le N.$ (3.27)

4. Analysis of the Compact Difference Scheme

In this section, we analyze the unique solvability, convergence and stability of the difference scheme (3.22)-(3.27).

4.1. Solvability

Theorem 4.1. (Solvability) The compact difference scheme (3.22)-(3.27) is uniquely solvable.

Proof. Noting that the difference scheme (3.19)-(3.24) is linear, we only need to prove that the homogeneous system has only zero solutions. Therefore we can suppose $u^l=0,-n<l<k-1,f=\phi =\psi =\varphi =0$, and consider the following homogeneous problem

$A_h{\Delta }_{\tau }^{\alpha }{u}_{i,j}^{k-1+\sigma }+B_h{v}_{i,j}^{k-1+\sigma }+B_h{v}_{i,j}^{k-1+\sigma -n}=0$, (4.1)

$A_h{v}_{i,j}^{k-1+\sigma }=B_h{u}_{i,j}^{k-1+\sigma },$ (4.2)

$v_{i,j}^{k-1+\sigma -n}=B_h{u}_{i,j}^{k-1+\sigma -n}.$ (4.3)

Taking the inner product $\left(\cdot ,\cdot \right)$ with $u^{k-1+\sigma }$ and $v^{k-1+\sigma }$ on both sides of (4.1) and (4.2) respectively, we have

$\left({\Delta }_{\tau }^{\alpha }{B}_h{v}^{k-1+\sigma },u^{k-1+\sigma }\right)+\left(B_h{v}^{k-1+\sigma },u^{k-1+\sigma }\right)+\left(B_h{v}^{k-1+\sigma -n},u^{k-1+\sigma }\right)=0,$ (4.4)

$\left(A_h{v}^{k-1+\sigma },v^{k-1+\sigma }\right)=\left(B_h{u}^{k-1+\sigma },v^{k-1+\sigma }\right).$ (4.5)

And then by Lemma 2.2 and (4.5), we have

$\left(B_h{v}^{k-1+\sigma },u^{k-1+\sigma }\right)=\left(B_h{u}^{k-1+\sigma },v^{k-1+\sigma }\right)=\left(A_h{v}^{k-1+\sigma },v^{k-1+\sigma }\right).$ (4.6)

Substituting (4.6) into (4.4), we have

$\left({\Delta }_{\tau }^{\alpha }{A}_h{u}^{k-1+\sigma },u^{k-1+\sigma }\right)+\left(A_h{v}^{k-1+\sigma },v^{k-1+\sigma }\right)+\left(B_h{v}^{k-1+\sigma -n},u^{k-1+\sigma }\right)=0.$ (4.7)

Using (2.18) and noticing $u^0=u^1=\cdots =u^{k-1}=0$, we have

$\left({\Delta }_{\tau }^{\alpha }{A}_h{u}^{k-1+\sigma },u^{k-1+\sigma }\right)=\frac{\sigma {\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}{c}_0^{\left(k\right)}\left(A_h{u}^k,u^k\right),$ (4.8)

$\left(A_h{v}^{k-1+\sigma },v^{k-1+\sigma }\right)=\sigma \left(A_h{v}^k,v^k\right),$ (4.9)

$\left(B_h{v}^{k-1+\sigma -n},u^{k-1+\sigma }\right)=0.$ (4.10)

Substituting (4.8)-(4.10) into (4.7) and using Lemma 2.4, we have

$\frac{4\sigma {\tau }^{-\alpha }}{9\Gamma \left(2-\alpha \right)}{c}_0^{\left(k\right)}{\Vert u^k\Vert }^2+\frac{4\sigma }{9}{\Vert v^k\Vert }^2\le \frac{\sigma {\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}{c}_0^{\left(k\right)}\left(A_h{u}^k,u^k\right)+\sigma \left(A_h{v}^k,v^k\right)=0.$ (4.11)

From Lemma 2.4 and noticing $c_0^{\left(k\right)}>0$, we get $\Vert u^k\Vert =\Vert v^k\Vert =0$, which immediately gives $u^k=0$, $v^k=0$.

The proof is completed. ◻

4.2. Convergence

Lemma 4.2. [43] Let $\left\{z_k\right\},\left\{g_k\right\}$ be two nonnegative sequences. If

$z_k\le K+{\displaystyle {\sum }_{i=0}^k{g}_i{z}_i},k\ge 0$, (4.12)

it holds that

$z_k\le K\cdot \exp \left({\displaystyle {\sum }_{i=0}^k{g}_i}\right),k\ge 0.$ (4.13)

where K is a nonnegative constant.

Theorem 4.2. (Convergence) Let $\left\{\left(U_{i,j}^k,V_{i,j}^k\right)|0\le i\le M_1,0\le j\le M_2,-n\le k\le N\right\}$ be the solution of problem (3.1)-(3.6) and $\left\{\left(u_{i,j}^k,v_{i,j}^k\right)|0\le i\le M_1,0\le j\le M_2,-n\le k\le N\right\}$ be the solution of the difference scheme (3.22)-(3.27) have

${\max }_{0\le k\tau \le T}\left\{{\Vert e^k\Vert }^2+{\Vert {\hat{e}}^{k-1+\sigma }\Vert }^2\right\}\le C{\left({\tau }^2+h_x^4+h_y^4\right)}^2.$ (4.14)

where $e_{i,j}^k=u_{i,j}^k-U_{i,j}^k$, ${\hat{e}}_{i,j}^k=v_{i,j}^k-V_{i,j}^k$.

Proof. Subtracting Equations (3.18)-(3.20) from (3.22)-(3.27), we can get the following error equations

$\begin{array}{l}{A}_h{\Delta }_{\tau }^{\alpha }{e}_{i,j}^{k-1+\sigma }+B_h{\hat{e}}_{i,j}^{k-1+\sigma }+B_h{\hat{e}}_{i,j}^{k-1+\sigma -n}\\ =A_h{\tilde{F}}_{i,j}^{k-1+\sigma }+{\left(R_1\right)}_{i,j}^k,1\le i\le M_1-1,1\le j\le M_2-1,1\le k\le N-1,\end{array}$ (4.15)

$A_h{\hat{e}}_{i,j}^{k-1+\sigma }=B_h{e}_{i,j}^{k-1+\sigma }+{\left(R_2\right)}_{i,j}^k,$ (4.16)

$A_h{\hat{e}}_{i,j}^{k-1+\sigma -n}=B_h{e}_{i,j}^{k-1+\sigma -n}+{\left(R_3\right)}_{i,j}^k,$ (4.17)

$e_{i,j}^k=0,0\le i\le M_1,0\le j\le M_2,-n\le k\le 0,$ (4.18)

$e_{i,j}^k=0,{\hat{e}}_{i,j}^k=0,i=0\text{or}{M}_1,j=0\text{or}{M}_2,1\le k\le N,$ (4.19)

$\begin{array}{c}{\tilde{F}}_{i,j}^{k-1+\sigma }=f\left(x_i,y_j,t_{k-1+\sigma },\left(\sigma +1\right)u_{i,j}^{k-1}-\sigma u_{i,j}^{k-2},\sigma u_{i,j}^{k-n}+\left(1-\sigma \right)u_{i,j}^{k-n-1}\right)\\ -f\left(x_i,y_j,t_{k-1+\sigma },\left(\sigma +1\right)U_{i,j}^{k-1}-\sigma U_{i,j}^{k-2},\sigma U_{i,j}^{k-n}+\left(1-\sigma \right)U_{i,j}^{k-n-1}\right).\end{array}$ (4.20)

Taking the inner product $\left(\cdot ,\cdot \right)$ with $A_h{\text{e}}^{k-1+\sigma }$ on both sides of (4.15), and with $A_h{\hat{e}}^{k-1+\sigma }$, $A_h{\hat{e}}^{k-1+\sigma -n}$ on both sides of (4.16) and (4.17) respectively, we have

$\begin{array}{l}\left(A_h{\Delta }_{\tau }^{\alpha }{e}^{k-1+\sigma },A_h{e}^{k-1+\sigma }\right)+\left(B_h{\hat{e}}^{k-1+\sigma },A_h{e}^{k-1+\sigma }\right)+\left(B_h{\hat{e}}^{k-1+\sigma -n},A_h{e}^{k-1+\sigma }\right)\\ =\left(A_h{\tilde{F}}^{k-1+\sigma },A_h{e}^{k-1+\sigma }\right)+\left(R_1^k,A_h{e}^{k-1+\sigma }\right),\end{array}$ (4.21)

$\left(A_h{\hat{e}}^{k-1+\sigma },A_h{\hat{e}}^{k-1+\sigma }\right)=\left(B_h{e}^{k-1+\sigma },A_h{\hat{e}}^{k-1+\sigma }\right)+\left(R_2^k,A_h{\hat{e}}^{k-1+\sigma }\right),$ (4.22)

$\left(A_h{\hat{e}}^{k-1+\sigma -n},A_h{\hat{e}}^{k-1+\sigma -n}\right)=\left(B_h{e}^{k-1+\sigma -n},A_h{\hat{e}}^{k-1+\sigma -n}\right)+\left(R_3^k,A_h{\hat{e}}^{k-1+\sigma -n}\right)$. (4.23)

Using Lemma 2.4 we get

$\left(B_h{\hat{e}}^{k-1+\sigma },A_h{e}^{k-1+\sigma }\right)=\left(A_h{\hat{e}}^{k-1+\sigma },A_h{\hat{e}}^{k-1+\sigma }\right)-\left(R_2^k,A_h{\hat{e}}^{k-1+\sigma }\right)$ (4.24)

and

$\left(B_h{\hat{e}}^{k-1+\sigma -n},A_h{e}^{k-1+\sigma }\right)=\left(A_h{\hat{e}}^{k-1+\sigma -n},A_h{\hat{e}}^{k-1+\sigma -n}\right)-\left(R_3^k,A_h{\hat{e}}^{k-1+\sigma -n}\right).$ (4.25)

Substituting (4.28)-(4.29) into (4.21), we have

$\begin{array}{l}\left({\Delta }_{\tau }^{\alpha }{A}_h{e}^{k-1+\sigma },A_h{e}^{k-1+\sigma }\right)+\left(A_h{\hat{e}}^{k-1+\sigma },A_h{\hat{e}}^{k-1+\sigma }\right)+\left(A_h{\hat{e}}^{k-1+\sigma -n},A_h{\hat{e}}^{k-1+\sigma -n}\right)\\ =\left(A_h{\tilde{F}}^{k-1+\sigma },A_h{e}^{k-1+\sigma }\right)+\left(R_1^k,A_h{e}^{k-1+\sigma }\right)+\left(R_2^k,A_h{\hat{e}}^{k-1+\sigma }\right)+\left(R_3^k,A_h{\hat{e}}^{k-1+\sigma -n}\right).\end{array}$ (4.26)

Using Lemma 2.6, we get

$\left({\Delta }_{\tau }^{\alpha }{A}_h{e}^{k-1+\sigma },A_h{e}^{k-1+\sigma }\right)\ge \frac{1}{2}{\Delta }_{\tau }^{\alpha }{\Vert A_h{e}^{k-1+\sigma }\Vert }^2.$ (4.27)

Then substituting the above inequality into (4.26), we can get the following inequality

$\begin{array}{l}\frac{1}{2}{\Delta }_{\tau }^{\alpha }{\Vert A_h{e}^{k-1+\sigma }\Vert }^2+{\Vert A_h{\hat{e}}^{k-1+\sigma }\Vert }^2+{\Vert \Delta A_h{\hat{e}}^{k-1+\sigma -n}\Vert }^2\\ \le \left(A_h{\tilde{F}}^{k-1+\sigma },A_h{e}^{k-1+\sigma }\right)+\left(R_1^k,A_h{e}^{k-1+\sigma }\right)+\left(R_2^k,A_h{\hat{e}}^{k-1+\sigma }\right)+\left(R_3^k,A_h{\hat{e}}^{k-1+\sigma -n}\right).\end{array}$ (4.28)

Applying Cauchy inequality to terms on the right-hand side of (4.28), we can get

$\left(A_h{\tilde{F}}^{k-1+\sigma },A_h{e}^{k-1+\sigma }\right)\le \frac{1}{2}{\Vert A_h{\tilde{F}}^{k-1+\sigma }\Vert }^2+\frac{1}{2}{\Vert A_h{e}^{k-1+\sigma }\Vert }^2,$ (4.29)

$\left(R_1,A_h{e}^{k-1+\sigma }\right)\le \frac{1}{2}{\Vert R_1^k\Vert }^2+\frac{1}{2}{\Vert A_h{e}^{k-1+\sigma }\Vert }^2,$ (4.30)

$\left(R_2,A_h{\hat{e}}^{k-1+\sigma }\right)\le \frac{1}{2}{\Vert R_2^k\Vert }^2+\frac{1}{2}{\Vert A_h{\hat{e}}^{k-1+\sigma }\Vert }^2,$ (4.31)

$\left(R_3,A_h{\hat{e}}^{k-1+\sigma -n}\right)\le \frac{1}{2}{\Vert R_3^k\Vert }^2+\frac{1}{2}{\Vert A_h{\hat{e}}^{k-1+\sigma -n}\Vert }^2.$ (4.32)

Substituting (4.29)-(4.32) into (4.28), we have

$\begin{array}{l}\frac{1}{2}{\Delta }_{\tau }^{\alpha }{\Vert A_h{e}^{k-1+\sigma }\Vert }^2+\frac{1}{2}{\Vert A_h{\hat{e}}^{k-1+\sigma }\Vert }^2\\ \le \frac{1}{2}{\Vert A_h{\tilde{F}}^{k-1+\sigma }\Vert }^2+{\Vert A_h{e}^{k-1+\sigma }\Vert }^2+\frac{1}{2}{\Vert R_1^k\Vert }^2+\frac{1}{2}{\Vert R_2^k\Vert }^2+\frac{1}{2}{\Vert R_3^k\Vert }^2.\end{array}$ (4.33)

Noticing that ${\Vert e^k\Vert }_{\infty }=0,-n\le k\le 0$, we have ${\Vert e^k\Vert }_{\infty }\le \frac{{ϵ}_0}{2},-n\le k\le 0$. For $0\le k\le l,0\le l\le N-1$, let

${\Vert e^k\Vert }_{\infty }\le \frac{{ϵ}_0}{2}.$ (4.34)

Then we will prove that (4.34) is also true for $k=l+1$.

According to (4.34), we have

$\begin{array}{c}\left|{\tilde{F}}_{i,j}^{k-1+\sigma }\right|=|f\left(x_i,y_j,t_{k-1+\sigma },\left(\sigma +1\right)u_{i,j}^{k-1}-\sigma u_{i,j}^{k-2},\sigma u_{i,j}^{k-n}+\left(1-\sigma \right)u_{i,j}^{k-n-1}\right)\\ -f\left(x_i,y_j,t_{k-1+\sigma },\left(\sigma +1\right)U_{i,j}^{k-1}-\sigma U_{i,j}^{k-2},\sigma U_{i,j}^{k-n}+\left(1-\sigma \right)U_{i,j}^{k-n-1}\right)|\\ \le c_1\left|\left(\sigma +1\right)e_{i,j}^{k-1}-\sigma e_{i,j}^{k-2}\right|+c_2\left|\sigma e_{i,j}^{k-n}+\left(1-\sigma \right)e_{i,j}^{k-n-1}\right|\\ \le c_1\left(\sigma +1\right)\left|e_{i,j}^{k-1}\right|+c_1\sigma \left|e_{i,j}^{k-2}\right|+c_2\sigma \left|e_{i,j}^{k-n}\right|+c_2\left(1-\sigma \right)\left|e_{i,j}^{k-n-1}\right|.\end{array}$ (4.35)

Using Lemma 2.3 and inequality (4.35), we can obtain that

$\begin{array}{c}{\Vert A_h{\tilde{F}}^{k-1+\sigma }\Vert }^2\le 2c_1^2{\left(\sigma +1\right)}^2{\Vert e^{k-1}\Vert }^2+2c_1^2{\sigma }^2{\Vert e^{k-2}\Vert }^2\\ +2c_2^2{\sigma }^2{\Vert e^{k-n}\Vert }^2+2c_2^2{\left(1-\sigma \right)}^2{\Vert e^{k-n-1}\Vert }^2,\end{array}$ (4.36)

and

${\Vert A_h{e}^{k-1+\sigma }\Vert }^2\le {\Vert e^{k-1+\sigma }\Vert }^2\le 2{\left(1+\sigma \right)}^2{\Vert e^{k-1}\Vert }^2+2{\sigma }^2{\Vert e^{k-2}\Vert }^2.$ (4.37)

Using (2.18) and substituting (4.36)-(4.37) into (4.33), we get

$\begin{array}{l}{c}_0^{\left(k\right)}{\Vert A_h{e}^k\Vert }^2+\mu {\Vert A_h{\hat{e}}^{k-1+\sigma }\Vert }^2\\ \le \underset{m=1}{\overset{k-1}{{\displaystyle \sum }}}\left(c_{k-m-1}^{\left(k\right)}-c_{k-m}^{\left(k\right)}\right){\Vert A_h{e}^m\Vert }^2+c_{k-1}^{\left(k\right)}{\Vert A_h{e}^0\Vert }^2+2\mu \left(c_1^2+2\right){\left(\sigma +1\right)}^2{\Vert e^{k-1}\Vert }^2\\ +2\mu \left(c_1^2+2\right){\sigma }^2{\Vert e^{k-2}\Vert }^2+2\mu c_2^2{\sigma }^2{\Vert e^{k-n}\Vert }^2+2\mu c_2^2{\left(1-\sigma \right)}^2{\Vert e^{k-n-1}\Vert }^2\\ +\mu {\Vert R_1^k\Vert }^2+\mu {\Vert R_2^k\Vert }^2+\mu {\Vert R_3^k\Vert }^2,\end{array}$ (4.38)

where

$\begin{array}{c}\mu ={\tau }^{\alpha }\Gamma \left(2-\alpha \right)=T^{\alpha }\Gamma \left(1-\alpha \right)\left(1-\alpha \right)N^{-\alpha }\\ <T^{\alpha }\Gamma \left(1-\alpha \right)\left(1-\alpha \right){\left(k-\frac{\alpha }{2}\right)}^{-\alpha }\\ <2c_{k-1}^k{T}^{\alpha }\Gamma \left(1-\alpha \right).\end{array}$ (4.39)

Substituting (4.39) into (4.38), we get

$\begin{array}{l}{c}_0^{\left(k\right)}{\Vert A_h{e}^k\Vert }^2+\mu {\Vert A_h{\hat{e}}^{k-1+\sigma }\Vert }^2\\ \le \underset{m=1}{\overset{k-1}{{\displaystyle \sum }}}\left(c_{k-m-1}^{\left(k\right)}-c_{k-m}^{\left(k\right)}\right){\Vert A_h{e}^m\Vert }^2+c_{k-1}^{\left(k\right)}[{\Vert A_h{e}^0\Vert }^2\end{array}$

$\begin{array}{l}+2T^{\alpha }\Gamma \left(1-\alpha \right)(2\left(c_1^2+2\right){\left(\sigma +1\right)}^2{\Vert e^{k-1}\Vert }^2+2\left(c_1^2+2\right){\sigma }^2{\Vert e^{k-2}\Vert }^2\\ +2c_2^2{\sigma }^2{\Vert e^{k-n}\Vert }^2+2c_2^2{\left(1-\sigma \right)}^2{\Vert e^{k-n-1}\Vert }^2+{\Vert R_1^k\Vert }^2+{\Vert R_2^k\Vert }^2+{\Vert R_3^k\Vert }^2)].\end{array}$ (4.40)

Using Lemma 2.3 and Lemma 2.7, we can get the following inequality

$\begin{array}{l}\frac{1}{9}{c}_0^{\left(k\right)}{\Vert e^k\Vert }^2+\frac{\mu }{9}{\Vert {\hat{e}}^{k-1+\sigma }\Vert }^2\\ \le \underset{m=1}{\overset{k-1}{{\displaystyle \sum }}}\left(c_{k-m-1}^{\left(k\right)}-c_{k-m}^{\left(k\right)}\right){\Vert e^m\Vert }^2+c_0^{\left(k\right)}[2T^{\alpha }\Gamma \left(1-\alpha \right)(2\left(c_1^2+2\right){\left(\sigma +1\right)}^2{\Vert e^{k-1}\Vert }^2\\ +2\left(c_1^2+2\right){\sigma }^2{\Vert e^{k-2}\Vert }^2+2c_2^2{\sigma }^2{\Vert e^{k-n}\Vert }^2+2c_2^2{\left(1-\sigma \right)}^2{\Vert e^{k-n-1}\Vert }^2\\ +{\Vert R_1^k\Vert }^2+{\Vert R_2^k\Vert }^2+{\Vert R_3^k\Vert }^2]\end{array}$ (4.41)

Dividing (4.41) by $\frac{1}{9}{c}_0^{\left(k\right)}$ on both sides, we have

$\begin{array}{l}{\Vert e^k\Vert }^2+\frac{\mu }{c_0^{\left(k\right)}}{\Vert {\hat{e}}^{k-1+\sigma }\Vert }^2\\ \le \underset{m=1}{\overset{k-1}{{\displaystyle \sum }}}\frac{9}{c_0^{\left(k\right)}}\left(c_{k-m-1}^{\left(k\right)}-c_{k-m}^{\left(k\right)}\right){\Vert e^m\Vert }^2+9[2T^{\alpha }\Gamma \left(1-\alpha \right)(2\left(c_1^2+2\right){\left(\sigma +1\right)}^2{\Vert e^{k-1}\Vert }^2\\ +2\left(c_1^2+2\right){\sigma }^2{\Vert e^{k-2}\Vert }^2+2c_2^2{\sigma }^2{\Vert e^{k-n}\Vert }^2+2c_2^2{\left(1-\sigma \right)}^2{\Vert e^{k-n-1}\Vert }^2\\ +{\Vert R_1^k\Vert }^2+{\Vert R_2^k\Vert }^2+{\Vert R_3^k\Vert }^2].\end{array}$ (4.42)

Then, letting

$C_1=36T^{\alpha }\Gamma \left(1-\alpha \right)\max \left\{\left(c_1^2+2\right){\left(\sigma +1\right)}^2,\left(c_1^2+2\right){\sigma }^2,c_2^2{\sigma }^2,c_2^2{\left(1-\sigma \right)}^2,{\hat{c}}^2/2\right\}$,(4.43)

we can get the following inequality

$\begin{array}{l}{\Vert e^k\Vert }^2+\frac{\mu }{c_0^{\left(k\right)}}{\Vert {\hat{e}}^{k-1+\sigma }\Vert }^2\\ \le \underset{m=1}{\overset{k-1}{{\displaystyle \sum }}}\frac{9}{c_0^{\left(k\right)}}\left(c_{k-m-1}^{\left(k\right)}-c_{k-m}^{\left(k\right)}\right){\Vert e^m\Vert }^2+C_1[{\Vert e^{k-1}\Vert }^2+{\Vert e^{k-2}\Vert }^2\\ +{\Vert e^{k-n}\Vert }^2+{\Vert e^{k-n-1}\Vert }^2+{\left({\tau }^2+h_x^4+h_y^4\right)}^2].\end{array}$ (4.44)

From (2.18), we have $c_0^{\left(k\right)}\le a_0+b_1$. So inequality (4.44) can be arranged as follows

$\begin{array}{l}{\Vert e^k\Vert }^2+\frac{\mu }{a_0+b_1}{\Vert {\hat{e}}^{k-1+\sigma }\Vert }^2\\ \le \underset{m=1}{\overset{k-1}{{\displaystyle \sum }}}\frac{9}{c_0^{\left(k\right)}}\left(c_{k-m-1}^{\left(k\right)}-c_{k-m}^{\left(k\right)}\right){\Vert e^m\Vert }^2+C_1[{\Vert e^{k-1}\Vert }^2+{\Vert e^{k-2}\Vert }^2\\ +{\Vert e^{k-n}\Vert }^2+{\Vert e^{k-n-1}\Vert }^2+{\left({\tau }^2+h_x^4+h_y^4\right)}^2].\end{array}$ (4.45)

Using Lemma 2.7, we know that $\frac{9}{c_0^{\left(k\right)}}\left(c_{k-m-1}^{\left(k\right)}-c_{k-m}^{\left(k\right)}\right)$ is a nonnegative sequence. Then applying Lemma 4.2, we get

$\begin{array}{l}{\Vert e^k\Vert }^2+\frac{\mu }{a_0+b_1}{\Vert {\hat{e}}^{k-1+\sigma }\Vert }^2\\ \le C_1{\left({\tau }^2+h_x^4+h_y^4\right)}^2\exp \left[4C_1+\underset{m=1}{\overset{k-1}{{\displaystyle \sum }}}\frac{9}{c_0^{\left(k\right)}}\left(c_{k-m-1}^{\left(k\right)}-c_{k-m}^{\left(k\right)}\right)\right]\\ =C_1\exp \left[4C_1+9\left(1-\frac{c_{k-1}^{\left(k\right)}}{c_0^{\left(k\right)}}\right)\right]\cdot {\left({\tau }^2+h_x^4+h_y^4\right)}^2\\ \le C_2^2{\left({\tau }^2+h_x^4+h_y^4\right)}^2,\end{array}$ (4.46)

where $C_2=\sqrt{C_1\exp \left(4C_1+9\right)}$. Letting $C_3=\min \left\{1,\frac{\mu }{a_0+b_1}\right\}$ and $C=\frac{C_2^2}{C_3}$, we get

${\Vert e^k\Vert }^2+{\Vert {\hat{e}}^{k-1+\sigma }\Vert }^2\le C{\left({\tau }^2+h_x^4+h_y^4\right)}^2,1\le k\le l+1.$ (4.47)

Let $h=\max \left\{h_x,h_y\right\}$ and assume that $\tau =O\left(h^2\right)$. Then for sufficiently smell h, we obtain

${\Vert e^{l+1}\Vert }_{\infty }\le c{h}^{-1}\Vert e^{l+1}\Vert \le c{h}^{-1}\left({\tau }^2+h^4\right)=O\left(h^3\right)\le \frac{{ϵ}_0}{2}$, (4.48)

where $c=C^{\frac{1}{2}}$. Thus we know that (4.34) holds for $k=l+1$. According to mathematical induction, then (4.34) holds for all $1\le k\le N$ and the theorem is proved. ◻

4.3. Stability

Let $\left\{p_{i,j}^k|0\le i\le M_1,0\le j\le M_2,-n\le k\le N\right\}$ and $\left\{q_{i,j}^k|0\le i\le M_1,0\le j\le M_2,-n\le k\le N\right\}$ are the solution of the following system

$\begin{array}{l}{A}_h{\Delta }_{\tau }^{\alpha }{p}_{i,j}^{k-1+\sigma }+B_h{q}_{i,j}^{k-1+\sigma }+B_h{q}_{i,j}^{k-1+\sigma -n}\\ =A_{h}f\left(x_i,y_j,t_{k-1+\sigma },\left(\sigma +1\right)p_{i,j}^{k-1}-\sigma p_{i,j}^{k-2},\sigma p_{i,j}^{k-n}+\left(1-\sigma \right)p_{i,j}^{k-n-1}\right),\end{array}$ (4.49)

$A_h{q}_{i,j}^{k-1+\sigma }=B_h{p}_{i,j}^{k-1+\sigma },$ (4.50)

$A_h{q}_{i,j}^{k-1+\sigma -n}=B_h{p}_{i,j}^{k-1+\sigma -n},1\le i\le M_1-1,1\le j\le M_2-1,1\le k\le N-1,$ (4.51)

$p_{i,j}^k=\phi \left(x_i,y_j,t_k\right)+{\rho }_{i,j}^k,\left(i,j\right)\in \omega ,-n\le k\le 0,$ (4.52)

$q_{i,j}^k=\psi \left(x_i,y_j,t_k\right),\left(i,j\right)\in \partial \omega ,1\le k\le N,$ (4.53)

$p_{i,j}^k=\varphi \left(x_i,y_j,t_k\right),\left(i,j\right)\in \partial \omega ,1\le k\le N,$ (4.54)

where ${\rho }_{i,j}^k$ is the perturbation of $\phi \left(x_i,y_j,t_k\right)$.

Let ${\theta }_{i,j}^k=u_{i,j}^k-p_{i,j}^k$, ${\hat{\theta }}_{i,j}^k=v_{i,j}^k-q_{i,j}^k$, $0\le i\le M_1,0\le j\le M_2,-n\le k\le N$.

Theorem 4.3. (Stability) There exists a positive integer $c_5$, such that

$\underset{0\le k\tau \le T}{\max }\left\{{\Vert {\theta }^k\Vert }^2+{\Vert {\hat{\theta }}^{k-1+\sigma }\Vert }^2\right\}\le c_5\underset{-n\le k\le 0}{\max }{\Vert {\rho }^k\Vert }^2.$ (4.55)

Proof. We subtract (3.22)-(3.27) from (4.49)-(4.54), we have

$A_h{\Delta }_{\tau }^{\alpha }{\theta }_{i,j}^{k-1+\sigma }+B_h{\hat{\theta }}_{i,j}^{k-1+\sigma }+B_h{\hat{\theta }}_{i,j}^{k-1+\sigma -n}=A_h{\hat{F}}_{i,j}^{k-1+\sigma },$ (4.56)

$A_h{\hat{\theta }}_{i,j}^{k-1+\sigma }=B_h{\theta }_{i,j}^{k-1+\sigma },$ (4.57)

$A_h{\hat{\theta }}_{i,j}^{k-1+\sigma -n}=B_h{\theta }_{i,j}^{k-1+\sigma -n},$ (4.58)

${\theta }_{i,j}^k={\rho }_{i,j}^k,0\le i\le M_1,0\le j\le M_2,-n\le k\le 0$. (4.59)

where

$\begin{array}{c}{\hat{F}}_{i,j}^{k-1+\sigma }=f\left(x_i,y_j,t_{k-1+\sigma },\left(\sigma +1\right)u_{i,j}^{k-1}-\sigma u_{i,j}^{k-2},\sigma u_{i,j}^{k-n}+\left(1-\sigma \right)u_{i,j}^{k-n-1}\right)\\ -f\left(x_i,y_j,t_{k-1+\sigma },\left(\sigma +1\right)p_{i,j}^{k-1}-\sigma p_{i,j}^{k-2},\sigma p_{i,j}^{k-n}+\left(1-\sigma \right)p_{i,j}^{k-n-1}\right).\end{array}$ (4.60)

Taking the inner product $\left(\cdot ,\cdot \right)$ with $A_h{\theta }^{k-1+\sigma }$ on both side of (4.56), and with $A_h{\hat{\theta }}^{k-1+\sigma }$ on both side of (4.57)-(4.58) respectively, we get

$\begin{array}{l}\left(A_h{\Delta }_{\tau }^{\alpha }{\theta }^{k-1+\sigma },A_h{\theta }^{k-1+\sigma }\right)+\left(B_h{\hat{\theta }}^{k-1+\sigma },A_h{\theta }^{k-1+\sigma }\right)+\left(A_h{\hat{\theta }}^{k-1+\sigma -n},A_h{\theta }^{k-1+\sigma }\right)\\ =\left(A_h{\hat{F}}^{k-1+\sigma },A_h{\theta }^{k-1+\sigma }\right),\end{array}$ (4.61)

$\left(A_h{\hat{\theta }}^{k-1+\sigma },A_h{\hat{\theta }}^{k-1+\sigma }\right)=\left(B_h{\theta }^{k-1+\sigma },A_h{\hat{\theta }}^{k-1+\sigma }\right)$, (4.62)

$\left(A_h{\hat{\theta }}^{k-1+\sigma -n},A_h{\hat{\theta }}^{k-1+\sigma }\right)=\left(B_h{\theta }^{k-1+\sigma -n},A_h{\hat{\theta }}^{k-1+\sigma }\right)$. (4.63)

Similar to (4.24)-(4.39) for (4.61)-(4.63), we obtain:

$\begin{array}{l}{c}_0^{\left(k\right)}{\Vert A_h{\theta }^k\Vert }^2+\mu {\Vert A_h{\hat{\theta }}^{k-1+\sigma }\Vert }^2\\ \le \underset{m=1}{\overset{k-1}{{\displaystyle \sum }}}\left(c_{k-m-1}^{\left(k\right)}-c_{k-m}^{\left(k\right)}\right){\Vert A_h{\theta }^m\Vert }^2\\ +c_{k-1}^{\left(k\right)}[{\Vert A_h{\theta }^0\Vert }^2+2T^{\alpha }\Gamma \left(1-\alpha \right)\left(c_1^2+2\right){\left(\sigma +1\right)}^2{\Vert {\theta }^{k-1}\Vert }^2\\ +\left(c_1^2+2\right){\sigma }^2{\Vert {\theta }^{k-2}\Vert }^2+c_2^2{\sigma }^2{\Vert {\theta }^{k-n}\Vert }^2+c_2^2{\left(1-\sigma \right)}^2{\Vert {\theta }^{k-n-1}\Vert }^2].\end{array}$ (4.64)

By equation (4.59), we can get

${\Vert {\theta }^0\Vert }^2\le \underset{-n\le k\le 0}{\max }{\Vert {\rho }^k\Vert }^2.$ (4.65)

Applying Lemma 2.3 and Lemma 2.7, then dividing both sides by $\frac{1}{9}{c}_0^{\left(k\right)}$, we get the following inequality

$\begin{array}{l}{\Vert {\theta }^k\Vert }^2+\frac{\mu }{c_0^{\left(k\right)}}{\Vert {\hat{\theta }}^{k-1+\sigma }\Vert }^2\\ \le \underset{m=1}{\overset{k-1}{{\displaystyle \sum }}}\frac{9}{c_0^{\left(k\right)}}\left(c_{k-m-1}^{\left(k\right)}-c_{k-m}^{\left(k\right)}\right){\Vert {\theta }^m\Vert }^2+9[\underset{-n\le k\le 0}{\max }{\Vert {\rho }^k\Vert }^2\\ +2T^{\alpha }\Gamma \left(1-\alpha \right)(\left(c_1^2+2\right){\left(\sigma +1\right)}^2{\Vert {\theta }^{k-1}\Vert }^2+\left(c_1^2+2\right){\sigma }^2{\Vert {\theta }^{k-2}\Vert }^2\\ +c_2^2{\sigma }^2{\Vert {\theta }^{k-n}\Vert }^2+c_2^2{\left(1-\sigma \right)}^2{\Vert {\theta }^{k-n-1}\Vert }^2)].\end{array}$ (4.66)

According to the inequality $c_0^{\left(k\right)}\le a_0+b_1$. Therefore, inequality (4.66) can be tidied up

$\begin{array}{l}{\Vert {\theta }^k\Vert }^2+\frac{\mu }{a_0+b_1}{\Vert {\hat{\theta }}^{k-1+\sigma }\Vert }^2\\ \le \underset{m=1}{\overset{k-1}{{\displaystyle \sum }}}\frac{9}{c_0^{\left(k\right)}}\left(c_{k-m-1}^{\left(k\right)}-c_{k-m}^{\left(k\right)}\right){\Vert {\theta }^m\Vert }^2\\ +c_3\left(\underset{-n\le k\le 0}{\max }{\Vert {\rho }^k\Vert }^2+{\Vert {\theta }^{k-1}\Vert }^2+{\Vert {\theta }^{k-2}\Vert }^2+{\Vert {\theta }^{k-n}\Vert }^2+{\Vert {\theta }^{k-n-1}\Vert }^2\right).\end{array}$ (4.67)

Letting

$\begin{array}{c}{c}_3=18\max \{T^{\alpha }\Gamma \left(1-\alpha \right)\left(c_1^2+2\right){\left(\sigma +1\right)}^2,T^{\alpha }\Gamma \left(1-\alpha \right)\left(c_1^2+2\right){\sigma }^2,\stackrel{}{}\\ T^{\alpha }\Gamma \left(1-\alpha \right)c_2^2{\sigma }^2,T^{\alpha }\Gamma \left(1-\alpha \right)c_2^2{\left(1-\sigma \right)}^2,\frac{1}{2}\}.\end{array}$ (4.68)

Applying the Lemma 4.1, we can get

$\begin{array}{l}{\Vert {\theta }^k\Vert }^2+\frac{\mu }{a_0+b_1}{\Vert {\hat{\theta }}^{k-1+\sigma }\Vert }^2\\ \le c_3\underset{-n\le k\le 0}{\max }{\Vert {\rho }^k\Vert }^2\exp \left[4c_3+\underset{m=1}{\overset{k-1}{{\displaystyle \sum }}}\frac{9}{c_0^{\left(k\right)}}\left(c_{k-m-1}^{\left(k\right)}-c_{k-m}^{\left(k\right)}\right)\right]\\ =c_3\exp \left[4c_3+9\left(1-\frac{c_{k-1}^{\left(k\right)}}{c_0^{\left(k\right)}}\right)\right]\cdot \underset{-n\le k\le 0}{\max }{\Vert {\rho }^k\Vert }^2\\ \le c_4^2\underset{-n\le k\le 0}{\max }{\Vert {\rho }^k\Vert }^2,\end{array}$ (4.69)

where $c_4=\sqrt{c_3\exp \left(4c_3+9\right)}$, according to $C_3=\min \left\{1,\frac{\mu }{a_0+b_1}\right\}$ and let $c_5=\frac{c_4^2}{C_3}$, so we have

${\Vert {\theta }^k\Vert }^2+{\Vert {\hat{\theta }}^{k-1+\sigma }\Vert }^2\le c_5\underset{-n\le k\le 0}{\max }{\Vert {\rho }^k\Vert }^2$. (4.70)

◻

5. Numerical Experiments

In this section, we verify the validity and accuracy of the compact difference scheme by two numerical examples. We choose $\Omega ={\left(0,1\right)}^2$ and $T=1$ for all examples. The error between exact solution and numerical solution in discrete $L^{\infty }$ norm and discrete $L^2$ norm are given as follows

$E_{\infty }\left(h,\tau \right)={\Vert u^N-U^N\Vert }_{\infty },E_{L^2}\left(h,\tau \right)=\Vert u^N-U^N\Vert .$ (5.1)

The convergence orders are given

$Order\left(h\right)={\log }_2\left(\frac{E_{L^2}\left(2h,\tau \right)}{E_{L^2}\left(h,\tau \right)}\right),Order\left(\tau \right)={\log }_2\left(\frac{E_{L^2}\left(h,2\tau \right)}{E_{L^2}\left(h,\tau \right)}\right).$ (5.2)

Example 1. In this example, we choose

$f\left(x,y,t,u\left(x,y,t\right),u\left(x,y,t-s\right)\right)=u^2\left(x,y,t\right)+u\left(x,y,t-s\right)+G\left(x,y,t\right),$ (5.3)

where

$\begin{array}{c}G\left(x,y,t\right)=\frac{1}{\Gamma \left(5\right)}\Gamma \left(\alpha +5\right)t^4\sin \left(\pi x\right)\sin \left(\pi y\right)+4{\pi }^4{t}^{4+\alpha }\sin \left(\pi x\right)\sin \left(\pi y\right)\\ +4{\pi }^4{\left(t-0.2\right)}^{4+\alpha }\sin \left(\pi x\right)\sin \left(\pi y\right)-{\left(t^{4+\alpha }\sin \left(\pi x\right)\sin \left(\pi y\right)\right)}^2\\ -{\left(t-0.2\right)}^{4+\alpha }\sin \left(\pi x\right)\sin \left(\pi y\right),\end{array}$ (5.4)

$\varphi \left(x,y,t\right)=t^{4+\alpha }\sin \left(\pi x\right)\sin \left(\pi y\right),0\le x\le 1,0\le y\le 1,t\in \left[-0.2,0\right],$ (5.5)

$\phi \left(x,y,t\right)=\psi \left(x,y,t\right)=0,\left(x,y,t\right)\in \partial \Omega \times \left[0,1\right].$ (5.6)

The exact solution of the problem is $u\left(x,y,t\right)=t^{4+\alpha }\sin \left(\pi x\right)\sin \left(\pi y\right)$.

We use the compact difference scheme (3.22)-(3.24) to solve the above problem. The errors and convergence order in discrete $L^2$ norm and $L^{\infty }$ norm are given in Table 1 and Table 2. Obviously, the spatial accuracy of order $O\left(h^4\right)$ is consistent with our theoretical results. In Table 3 and Table 4, the numerical

Table 1. The computational error and convergence order in spatial dimension for U at T = 1 for Example 1 using present scheme.

Table 2. The computational error and convergence order in spatial dimension for V at T = 1 for Example 1 using present scheme.

Table 3. The computational error and convergence order in spatial dimension for U at T = 1 for Example 1 using central difference scheme.

Table 4. The computational error and convergence order in spatial dimension for V at T = 1 for Example 1 using central difference scheme.

(a) (b)

Figure 1. Exact and numerical solution surface with $\alpha =0.3,\tau =h^2=\frac{1}{50}$ for Example 1.

results of the central difference scheme are compared with those of the theoretical scheme. From the following Tables 1-4, it is easy to find that the compact difference scheme can achieve higher accuracy than the central difference scheme. (Figure 1)

Example 2. In this example, we choose

$\begin{array}{l}f\left(x,y,t,u\left(x,y,t\right),u\left(x,y,t-s\right)\right)\\ =4u\left(x,y,t\right)+u^2\left(x,y,t\right)+u\left(x,y,t-s\right)+G\left(x,y,t\right).\end{array}$ (5.7)

Where

$\begin{array}{c}G\left(x,y,t\right)=\left(\frac{\Gamma \left(\alpha +4\right)}{\Gamma \left(4\right)}-t^{3+2\alpha }\exp \left(x+y\right)\right)t^3\exp \left(x+y\right)\\ +3{\left(t-0.2\right)}^{3+\alpha }\exp \left(x+y\right),\end{array}$ (5.8)

$u\left(x,y,t\right)=t^{3+\alpha }\exp \left(x+y\right),0\le x\le 1,0\le y\le 1,t\in \left[-0.2,0\right],$ (5.9)

$u\left(0,y,t\right)=t^{3+\alpha }\exp \left(y\right),u\left(1,y,t\right)=t^{3+\alpha }\exp \left(1+y\right),$ (5.10)

$u\left(x,0,t\right)=t^{3+\alpha }\exp \left(x\right),u\left(x,1,t\right)=t^{3+\alpha }\exp \left(1+x\right),0\le t\le 1.$ (5.11)

The exact solution of the problem is $u\left(x,y,t\right)=t^{3+\alpha }\exp \left(x+y\right)$.

In this example, we use the formula (3.13) to calculate the nonlinear source term. The error of the $L^2$ and $L^{\infty }$ norms for $\alpha =0.3,0.6,0.9$ and convergence orders in the spatial directions are listed in Table 5 and Table 6, from which can easy to see that the convergence order in space of the compact difference scheme (3.22)-(3.24) we proposed reaches the fourth-order accuracy, which is in agreement with our theoretical results. In Table 7 and Table 8, we compare the numerical results of the central difference scheme with those of the theoretical scheme in the spatial direction. From these tables, we can see that the accuracy of the central difference scheme in the spatial direction is $O\left(h^2\right)$ which is less accurate than the compact difference scheme. (Figure 2)

Table 5. The computational error and convergence order in spatial dimension for U at T = 1 for Example 2 using present scheme.

Table 6. The computational error and convergence order in spatial dimension for V at T = 1 for Example 2 using present scheme.

Table 7. The computational error and convergence order in spatial dimension for U at T = 1 for Example 2 using central difference scheme.

Table 8. The computational error and convergence order in spatial dimension for V at T = 1 for Example 2 using central difference scheme.

(a) (b)

Figure 2. Exact and numerical solution surface with $\alpha =0.3,\tau =h^2=\frac{1}{40}$ for Example 2.

6. Conclusion

In this work, we constructed a linearized compact difference scheme for a non-linear sub-diffusion time-delay equation in two-dimensional space. The global convergence order of the scheme is $O\left({\tau }^2+h_x^4+h_y^4\right)$. We linearize the nonlinear term and prove the uniqueness of the scheme by proving that the corresponding homogeneous problem has only zero solutions. The convergence of the scheme under the discrete $L^2$ norm is obtained by the discrete energy method, and the stability of the scheme is also proved. Numerical experiments have been conducted to show the robustness and accuracy of the proposed numerical scheme. In a word, our scheme is very effective to solve a class of fourth-order nonlinear sub-diffusion neutral delayed equation. In further research, we will apply the proposed method to nonlinear neutral time-delay sub-diffusion equations with variable fractional order derivative.

Acknowledgements

The authors would like to thank editor and referees for their valuable advices for the improvement of this article.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Sokolov, I.M. and Klafter, J. (2005) From Diffusion to Anomalous Diffusion: A Century after Einstein’s Brownian Motion. Chaos: An Interdisciplinary Journal of Nonlinear Science, 15, Article ID: 026103. https://doi.org/10.1063/1.1860472
[2]	El-Nabulsi, R.A. (2008) Fractional Field Theories from Multi-Dimensional Fractional Variational Problems. International Journal of Geometric Methods in Modern Physics, 5, 863-892. https://doi.org/10.1142/S0219887808003119
[3]	Baleanu, D., Defterli, O. and Agrawal, O.P. (2009) A Central Difference Numerical Scheme for Fractional Optimal Control Problems. Journal of Vibration and Control, 15, 583-597. https://doi.org/10.1177/1077546308088565
[4]	Saxena, R.K., Mathai A.M. and Haubold, H.J. (2004) On Generalized Fractional Kinetic Equations. Physica A: Statistical Mechanics and Its Applications, 344, 657-664. https://doi.org/10.1016/j.physa.2004.06.048
[5]	Schumer, R., Benson, D.A., Meerschaert, M.M. and Wheatcraft, S.W. (2001) Eulerian Derivation of the Fractional Advection-Dispersion Equation. Journal of Contaminant Hydrology, 48, 69-88. https://doi.org/10.1016/S0169-7722(00)00170-4
[6]	Lu, S., Molz, F.J. and Fix, G.J. (2002) Possible Problems of Scale Dependency in Applications of the Three-Dimensional Fractional Advection-Dispersion Equation to Natural Porous Media. Water Resources Research, 38, 4-1-4-7. https://doi.org/10.1029/2001WR000624
[7]	Benson, D., Schumer, R., Meerschaert, M. and Wheatcraft, S. (2001) Fractional Dispersion, Lévy Motion, and the MADE Tracer Tests. In: Berkowitz, B., Ed., Dispersion in Heterogeneous Geological Formations, Springer, Dordrecht, 211-240. https://doi.org/10.1007/978-94-017-1278-1_11
[8]	Schneider, W.R. and Wyss, W. (1989) Fractional Diffusion and Wave Equations. Journal of Mathematical Physics, 30, 134-144. https://doi.org/10.1063/1.528578
[9]	Meerschaert, M.M., Benson, D.A., Scheffler, H.P. and Baeumer, B. (2002) Stochastic Solution of Space-Time Fractional Diffusion Equations. Physical Review E, 65, 1103-1106. https://doi.org/10.1103/PhysRevE.65.041103
[10]	Chen, C.M., Liu, F., Turner, I. and Anh, V.A. (2007) Fourier Method for the Fractional Diffusion Equation Describing Sub-Diffusion. Journal of Computational Physics, 227, 886-897. https://doi.org/10.1016/j.jcp.2007.05.012
[11]	Tomovski, Z. and Sandev, T. (2013) Exact Solutions for Fractional Diffusion Equation in a Bounded Domain with Different Boundary Conditions. Nonlinear Dynamics, 71, 671-683. https://doi.org/10.1007/s11071-012-0710-x
[12]	Langlands, T.A.M. and Henry, B.I. (2005) The Accuracy and Stability of an Implicit Solution Method for the Fractional Diffusion Equation. Journal of Computational Physics, 205, 719-736. https://doi.org/10.1016/j.jcp.2004.11.025
[13]	Yuste, S.B. and Acedo, L. (2005) An Explicit Finite Difference Method and a New Von-Neumann Type Stability Analysis for Fractional Diffusion Equations. SIAM Journal on Numerical Analysis, 42, 1862-1874. https://doi.org/10.1137/030602666
[14]	Tadjeran, C., Meerschaert, M.M. and Scheffler, H.P. (2006) A Second-Order Accurate Numerical Approximation for the Fractional Diffusion Equation. Journal of Computational Physics, 213, 205-213. https://doi.org/10.1016/j.jcp.2005.08.008
[15]	Zhuang, P., Liu, F., Anh, V. and Turner, I. (2008) New Solution and Analytical Techniques of the Implicit Numerical Method for the Anomalous Sub-Diffusion Equation. SIAM Journal on Numerical Analysis, 46, 1079-1095. https://doi.org/10.1137/060673114
[16]	Liu, F., Yang, C. and Burrage, K. (2009) Numerical Method and Analytical Technique of the Modified Anomalous Sub-Diffusion Equation with a Nonlinear Source Term. Journal of Computational and Applied Mathematics, 231, 160-176. https://doi.org/10.1016/j.cam.2009.02.013
[17]	Chen, J., Liu, F., Anh, V., Shen, S., Liu, Q., et al. (2012) The Analytical Solution and Numerical Solution of the Fractional Diffusion-Wave Equation with Damping. Applied Mathematics and Computation, 219, 1737-1748. https://doi.org/10.1016/j.amc.2012.08.014
[18]	Mohammed, O.H., Fadhel, S.F. and AL-Safi. M.G.S. (2015) Numerical Solution for the Time Fractional Diffusion-Wave Equations by Using Sinc-Legendre Collocation Method. Mathematical Theory and Modeling, 5, 49-57.
[19]	Zhou, F.Y. and Xu, X.Y. (2017) Numerical Solution of Time-Fractional Diffusion-Wave Equations via Chebyshev Wavelets Collocation Method. Advances in Mathematical Physics, 2017, Article ID: 2610804. https://doi.org/10.1155/2017/2610804
[20]	Cui, M. (2009) Compact Finite Difference Method for the Fractional Diffusion Equation. Journal of Computational Physics, 228, 7792-7804. https://doi.org/10.1016/j.jcp.2009.07.021
[21]	Chen, C.M., Liu, F., Anh, V. and Turner, I. (2010) Numerical Schemes with High Spatial Accuracy for a Variable Order Anomalous Sub-Diffusion Equations. SIAM Journal on Scientific Computing, 32, 1740-1760. https://doi.org/10.1137/090771715
[22]	Gao, G.H. and Sun, Z.Z. (2011) A Compact Difference Scheme for the Fractional Subdiffusion Equations. Journal of Computational Physics, 230, 586-595. https://doi.org/10.1016/j.jcp.2010.10.007
[23]	Vong, S., Lyu, P. and Wang, Z. (2015) A Compact Difference Scheme for Fractional Sub-Diffusion Equations with the Spatially Variable Coefficient under Neumann Boundary Conditions. Journal of Scientific Computing, 66, 725-739. https://doi.org/10.1007/s10915-015-0040-5
[24]	Zhao, X., Sun, Z.Z. and Karniadakis, G.E. (2015) Second-Order Approximations for Variable Order Fractional Derivatives: Algorithms and Applications. Journal of Computational Physics, 293, 184-200. https://doi.org/10.1016/j.jcp.2014.08.015
[25]	Gao, G.H., Sun, H.W. and Sun, Z.Z. (2015) Stability and Convergence of Finite Difference Schemes for a Class of Time-Fractional Sub-Diffusion Equations Based on Certain Superconvergence. Journal of Computational Physics, 280, 510-528. https://doi.org/10.1016/j.jcp.2014.09.033
[26]	Alikhanov, A.A. (2015) A New Difference Scheme for the Time Fractional Diffusion Equation. Journal of Computational Physics, 280, 424-438. https://doi.org/10.1016/j.jcp.2014.09.031
[27]	Sneddon, I.N. (1951) Fourier Transforms. McGraw Hill, New York.
[28]	Oldhan, K.B. and Spainer, J. (1974) The Fractional Calculus. Academic Press, New York.
[29]	Myers, T.G., Charpin, J.P.F. and Chapman, S.J. (2002) The Flow and Solidification of a Thin Fluid Film on an Arbitrary Three-Dimensional Surface. Physics of Fluids, 14, 2788-2803. https://doi.org/10.1063/1.1488599
[30]	Myers, T.G. and Charpin, J.P.F. (2004) A Mathematical Model for Atmospheric Ice Accretion and Water Flow on a Cold Surface. International Journal of Heat and Mass Transfer, 47, 5483-5500. https://doi.org/10.1016/j.ijheatmasstransfer.2004.06.037
[31]	Karpman, V.I. (1996) Stabilization of Soliton Instabilities by Higher-Order Dispersion: Fourth Order Nonlinear Schrödinger-Type Equations. Physical Review E, 53, 1336-1339. https://doi.org/10.1103/PhysRevE.53.R1336
[32]	Karpman, V.I. and Shagalov, A.G. (2000) Stability of Soliton Described by Nonlinear Schrödinger-Type Equations with Higher Order Dispersion. Physica D: Nonlinear Phenomena, 144, 194-210. https://doi.org/10.1016/S0167-2789(00)00078-6
[33]	Agrawal, O.P. (2001) A General Solution for a Fourth-Order Fractional Diffusion-Wave Equation Defined in A Bounded Domain. Computers and Structures, 79, 1497-1501. https://doi.org/10.1016/S0045-7949(01)00026-8
[34]	Hu, X.L. and Zhang, L.M. (2012) On Finite Difference Methods for Fourth-Order Fractional Diffusionwave and Sub-Diffusion Systems. Applied Mathematics and Computation, 218, 5019-5034. https://doi.org/10.1016/j.amc.2011.10.069
[35]	Guo, J., Li, C.P. and Ding, H.F. (2014) Finite Difference Methods for Time Subdiffusion Equation with Space Fourth Order. Commun. Applied Mathematics and Computation, 28, 96-108. (In Chinese)
[36]	Ji, C.-C., Sun, Z.-Z. and Hao, Z.-P. (2015) Numerical Algorithms with High Spatial Accuracy for the Fourth-Order Fractional Sub-Diffusion Equations with the First Dirichlet Boundary Conditions. Journal of Scientific Computing, 66, 1148-1174. https://doi.org/10.1007/s10915-015-0059-7
[37]	Zhang, P. and Pu, H. (2017) A Second-Order Compact Difference Scheme for the Fourth-Order Fractional Sub-Diffusion Equation. Numerical Algorithms, 76, 573-598. https://doi.org/10.1007/s11075-017-0271-7
[38]	Nandal, S. and Pandey, D.N. (2019) Numerical Solution of Time Fractional Non-Linear Neutral Delay Differential Equations of Fourth-Order. Malaya Journal of Matematik, 78, 1467-1487. https://doi.org/10.26637/MJM0703/0035
[39]	Nandal, S. and Pandey, D.N. (2020) Numerical Solution of Non-Linear Fourth-Order Fractional Subdiffusion Wave Equation with Time Delay. Applied Mathematics and Computation, 369, Article ID: 124900. https://doi.org/10.1016/j.amc.2019.124900
[40]	Pimenov, V.G., Hendy, A.S. and De Staelen, R.H. (2017) On a Class of Non-Linear Delay Distributed Order Fractional Diffusion Equations. Journal of Computational and Applied Mathematics, 318, 433-443. https://doi.org/10.1016/j.cam.2016.02.039
[41]	Sun, Z.Z. (2012) Numerical Methods for Partial Differential Equations. 2nd Edition, Science Press, Beijing. (In Chinese)
[42]	Li, Q., Yang, Q. and Chen, H.Z. (2020) Compact Difference Scheme for Two-Dimensional Fourth-Order Nonlinear Hyperbolic Equation. Numerical Methods for Partial Differential Equations, 36, 1938-1961. https://doi.org/10.1002/num.22511
[43]	Zhang, Q.F., Ran, M.H. and Xu, D.H. (2017) Analysis of the Compact Difference Scheme for the Semilinear Fractional Partial Differential Equation with Time Delay. Applicable Analysis, 96, 1867-1884. https://doi.org/10.1080/00036811.2016.1197914