Improved Mass- and Energy-Conserving Difference Scheme for Two-Dimensional Nonlinear Space-Fractional Schrödinger Equation

1. Introduction

In this paper, we will devoted to constructing a high-order numerical scheme for the following two-dimensional nonlinear space fractional Schrödinger equation

$$\begin{array}{c}\hfill \mathrm{i}{\displaystyle \frac{\partial u}{\partial t}}+\mu {\mathcal{L}}_{\alpha ,\beta }u+\rho {\left|u\right|}^{2}u=0,(x,y)\in \Omega ,t\in (0,T],\end{array}$$

(1)

subject to the initial condition,

$${\displaystyle u(x,y,0)=u_0(x,y),(x,y)\in \partial \Omega ,}$$

where $\Omega =(a,b)\times (c,d)$, $\mathrm{i}=\sqrt{-1}$, $\mu$ and $\rho$ are real constants, and $u_0(x,y)$ is a given smooth function. The operator ${\mathcal{L}}_{\alpha ,\beta }$ is defined by ${\mathcal{L}}_{\alpha ,\beta }:={\partial }_x^{\alpha }+{\partial }_y^{\beta }$ with $1<\alpha ,\beta \le 2$, and ${\partial }_x^{\alpha }$ and ${\partial }_y^{\beta }$ denote the Riesz fractional derivative operators in x and y direction, respectively, which are defined by [8]

$${\displaystyle {\partial }_x^{\alpha }u=-\frac{1}{2cos\left(\frac{\pi }{2}\alpha \right)\Gamma (2-\alpha )}\frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}{\int }_a^b{|x-s_1|}^{1-\alpha }u\left(s_1,y,t\right)\mathrm{d}{s}_1,1<\alpha \le 2,}$$

and

$${\displaystyle {\partial }_y^{\beta }u=-\frac{1}{2cos\left(\frac{\pi }{2}\beta \right)\Gamma (2-\beta )}\frac{{\mathrm{d}}^2}{\mathrm{d}{y}^2}{\int }_c^d{|y-s_2|}^{1-\beta }u\left(x,s_2,t\right)\mathrm{d}{s}_2,1<\beta \le 2.}$$

The equation (1) own two fundamental conservation properties of the mass and of the energy, i.e.,

$${\displaystyle \mathcal{M}\left(t\right)=\mathcal{M}\left(0\right),\mathcal{E}\left(t\right)=\mathcal{E}\left(0\right),0\le t\le T,}$$

where

$${\displaystyle \mathcal{M}\left(t\right)={\int }_{\Omega }{\left|u(x,y,t)\right|}^2\mathrm{d}x\mathrm{d}y,}$$

and

$${\displaystyle \mathcal{E}\left(t\right)={\int }_{\Omega }-\left[\mu u^{*}(x,y,t){\mathcal{L}}_{\alpha ,\beta }u(x,y,t)+\frac{\rho }{2}{\left|u(x,y,t)\right|}^4\right]\mathrm{d}x\mathrm{d}y,}$$

where $u^{*}(x,y,t)$ denotes the conjugate of $u(x,y,t)$.

In the past decades, some numerical methods have been developed to solve the nonlinear spatial fractional Schrödinger equation in one-dimensional case. For example, Wang and Huang constructed a structure-preserving difference scheme with convergence order of $\mathcal{O}\left({\tau }^2+h^2\right)$ in [8]. Then, a new finite difference scheme with the same convergence order as the previous algorithm was proposed by them in [10], and the conservation properties of this scheme were studied in detail, including the conservation of mass and energy. With the help of the scalar auxiliary variable method proposed by Shen and Xu [6]and the existing fractional center difference formula, Li et al. [5] gave a a fully discrete structure-preserving scheme with convergence order$\mathcal{O}\left({\tau }^2+h^2\right)$. For other related important and meaningful work, please refer to [2,4,7] and their references. Next, let’s review the case of two-dimensional equations. At present, there are actually some effective numerical schemes. For example, Zhao [12] and his collaborators proposed a kind of linearized energy conservation finite difference algorithm. a split-step spectral Galerkin scheme for this equation was proposed by Wang et al.[11], and its convergence was strictly proved. Based on the compact implicit integration factor method, an efficient finite difference scheme was constructed by Zhang et al. [13] and studied carefully.

Up to now, although the numerical methods for nonlinear spatial fractional Schrödinger equation have made great progress, the convergence order of these methods is still relatively low. Up to now, as far as we know, except for [3], the authors have established two linearly finite difference scheme with second-order accuracy in time and spectral accuracy in space, and almost no further work has focused on the scalar auxiliary variable method for the high-dimensional nonlinear space fractional Schrödinger equations. Therefore, inspired by the the scalar auxiliary variable method and the references [3,5], we will establish a higher-order numerical method to solve the two-dimensional nonlinear space fractional Schrödinger equations, whose convergence order can reach $\mathcal{O}\left({\tau }^2+h_x^4+h_y^4\right)$, and the established difference scheme maintains mass conservation and energy conservation.

The outline of this paper is arranged as follows: In Section 2, we proposed a high-order implicit conservative scheme for the two-dimensional nonlinear space fractional Schrödinger equation. In Section 3, we study the conservation properties of the established numerical algorithm.

2. Construction of the Conservative Scheme

Let $x_i=i{h}_1,i=0,1,...,M_1$, $y_j=j{h}_2,j=0,1,...,M_2$, and $t_k=k\tau ,k=0,1,...,N$, where $h_1={\displaystyle \frac{b-a}{M_1}},h_2={\displaystyle \frac{d-c}{M_2}},\tau ={\displaystyle \frac{T}{N}}$ are spatial and temporal step sizes, respectively. Denote any given grid function

$$u=\left\{u_{ij}^{k}0\le i\le M_1,0\le j\le M_2,0\le k\le N\right\},$$

we introduce the following notations:

$${\delta }_t{u}_{ij}^k={\displaystyle \frac{u_{ij}^{k+1}-u_{ij}^k}{\tau }},f_{ij}^{k+{\displaystyle \frac{1}{2}}}={\displaystyle \frac{1}{2}}\left(f_{ij}^{k+1}+f_{ij}^k\right),$$

$${\delta }_x{u}_{i-{\displaystyle \frac{1}{2}}}^n={\displaystyle \frac{1}{h}}\left(u_i^n-u_{i-1}^n\right),{\delta }_x^2{u}_i^n={\displaystyle \frac{1}{h}}\left({\delta }_x{u}_{i+{\displaystyle \frac{1}{2}}}^n-{\delta }_x{u}_{i-{\displaystyle \frac{1}{2}}}^n\right).$$

For any grid functions $u=\left\{u_{ij}\right\}$ and $v=\left\{v_{ij}\right\}$, we define the discrete inner product and its corresponding norm as

$$\begin{array}{c}{\displaystyle \left(u,v\right)=h_1{h}_2\sum _{i=1}^{M_1-1}\sum _{j=1}^{M_2-1}{u}_{ij}{v}_{ij}^{*},{∥u∥}^2=\left(u,v\right).}\hfill \end{array}$$

First, it is assumed that there is a constant $C_0$ such that ${\int }_{\Omega }F\left(\left|u\right(x,y,t{|}^2\right)\mathrm{d}x\mathrm{d}y+C_0>0$, where $F^{\prime }\left(u\right)=f\left(u\right)=u$, we introduce a scalar auxiliary variable

$${\displaystyle r\left(t\right):=\sqrt{{\int }_{\Omega }F\left({\left|u(x,y,t)\right|}^2\right)\mathrm{d}x\mathrm{d}y+C_0},}$$

which was introduced by Shen et.al. [6] to solve gradient flows, then there holds that

$${\displaystyle \frac{\mathrm{d}r\left(t\right)}{\mathrm{d}t}=\Re {\int }_{\Omega }Q\left(u\right)u_t^{*}\mathrm{d}x\mathrm{d}y,}$$

where $\Re \left(z\right)$ means to take the real part of a complex number z, and

$${\displaystyle Q\left(u\right):=\frac{{\left|u(x,y,t)\right|}^{2}u(x,y,t)}{\sqrt{{\int }_{\Omega }F\left(\left|u\right(x,y,t{|}^2\right)\mathrm{d}x\mathrm{d}y+C_0}}.}$$

At this point, equation (1) can be written in the following equivalent form

$$\begin{array}{c}\hfill \mathrm{i}{\displaystyle \frac{\partial u}{\partial t}}+\mu {\mathcal{L}}_{\alpha ,\beta }u+\rho r\left(t\right)Q\left(u\right)=0,(x,y)\in \Omega ,t\in (0,T],\end{array}$$

(2)

$$\begin{array}{c}\hfill {\displaystyle \frac{\mathrm{d}r\left(t\right)}{\mathrm{d}t}}=\Re {\int }_{\Omega }Q\left(u\right)u_t^{*}\mathrm{d}x\mathrm{d}y,t\in (0,T].\end{array}$$

(3)

Moreover, by applying analytical methods similar to the one-dimensional equation case in [5], it is not difficult to claim that the new systems (2)-(3) have the same mass and energy as the original equation (1).

Next, we will focus on developing a conservative scheme for the systems (2)-(3). For the space Riesz fractional derivatives, we employ a fourth-order numerical differential formula proposed by Ding and Li in [1], that is

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle {\mathcal{H}}_z^{\gamma }{\partial }_z^{\alpha }u(z,t)={\delta }_z^{\gamma }u(z,t)+\mathcal{O}\left(h_z^4\right),(z,\gamma )=(x,\alpha )\mathrm{or}(y,\beta ),}\hfill \end{array}\end{array}$$

(4)

where the operators ${\mathcal{H}}_z^{\gamma }$ and ${\delta }_z^{\gamma }u(z,t)$ are defined by

$$\begin{array}{c}{\displaystyle {\mathcal{H}}_z^{\gamma }u(z,t)=\left(1+\frac{h_z^2}{6}{\delta }_z^2\right)u(z,t),{\delta }_z^{\gamma }u(z,t)=-\frac{1}{2cos\left(\frac{\pi }{2}\gamma \right)}\left.{(}^L{\mathcal{A}}_{h_z}^{\gamma }{+}^R{\mathcal{A}}_{h_z}^{\gamma }\right)u(z,t),}\hfill \end{array}$$

where ^L

$$\begin{array}{c}{\displaystyle {\mathcal{A}}_{h_z}^{\gamma }u(z,t)=\frac{1}{h_z^{\gamma }}\sum _{m=0}^{\left[\frac{z-S_1}{h_z}\right]+1}{\kappa }_m^{\left(\gamma \right)}u\left(z-(m-1)h_z,t\right){,}^R{\mathcal{A}}_{h_z}^{\gamma }u(z,t)=\frac{1}{h_z^{\gamma }}\sum _{m=0}^{\left[\frac{S_2-z}{h_z}\right]+1}{\kappa }_m^{\left(\gamma \right)}u\left(z+(m-1)h_z,t\right).}\hfill \end{array}$$

in which $\left(S_1,S_2\right)=(a,b)$ for $z=x$, and $\left(S_1,S_2\right)=(c,d)$ for $z=y$. Besides, the coefficient ${\kappa }_m^{\left(\gamma \right)}$ can be also computed by the following recursion formulas [1],

$$\left\{\begin{array}{c}{\displaystyle {\kappa }_0^{\left(\gamma \right)}=b_0^{\gamma },}\hfill \\ {\displaystyle {\kappa }_1^{\left(\gamma \right)}=\frac{b_1\gamma }{b_0}{\kappa }_0^{\left(\gamma \right)},}\hfill \\ {\displaystyle {\kappa }_2^{\left(\gamma \right)}=\frac{1}{2b_0}\left[b_1(\gamma -1){\kappa }_1^{\left(\gamma \right)}+2b_2\gamma {\kappa }_0^{\left(\gamma \right)}\right],}\hfill \\ {\displaystyle {\kappa }_3^{\left(\gamma \right)}=\frac{1}{3b_0}\left[b_1(\gamma -2){\kappa }_2^{\left(\gamma \right)}+b_2(2\gamma -1){\kappa }_1^{\left(\gamma \right)}+3b_3\gamma {\kappa }_0^{\left(\gamma \right)}\right],}\hfill \\ {\displaystyle {\kappa }_m^{\left(\gamma \right)}=\frac{1}{m{b}_0}\left[b_1(\gamma +1-m){\kappa }_{m-1}^{\left(\gamma \right)}+b_2(2\gamma +2-m){\kappa }_{m-2}^{\left(\gamma \right)}\right.}\hfill \\ \left.+b_3(3\gamma +3-m){\kappa }_{m-3}^{\left(\gamma \right)}+b_4(4\gamma +4-m){\kappa }_{m-4}^{\left(\gamma \right)}\right],m\ge 4.\hfill \end{array}\right.$$

where

$$\begin{array}{c}{\displaystyle b_0=\frac{25{\gamma }^3-30{\gamma }^2+13\gamma -2}{12{\gamma }^3},b_1=\frac{-\left(24{\gamma }^3-43{\gamma }^2+23\gamma -4\right)}{6{\gamma }^3},}\hfill \\ {\displaystyle b_2=\frac{6{\gamma }^3-15{\gamma }^2+10\gamma -2}{2{\gamma }^3},b_3=\frac{-\left(8{\gamma }^3-21{\gamma }^2+17\gamma -4\right)}{6{\gamma }^3},b_4=\frac{3{\gamma }^3-8{\gamma }^2+7\gamma -2}{12{\gamma }^3}.}\hfill \end{array}$$

Hence, it follows from (4), we have

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle {\mathcal{H}}_x^{\alpha }{\mathcal{H}}_y^{\beta }{\mathcal{L}}_{\alpha ,\beta }u(x,y,t)={\mathcal{H}}_y^{\beta }{\delta }_x^{\alpha }u(x,y,t)+{\mathcal{H}}_x^{\alpha }{\delta }_y^{\beta }u(x,y,t)+\mathcal{O}\left(h_x^4+h_y^4\right).}\hfill \end{array}\end{array}$$

(5)

Considering the systems (2)-(3) at the grid point $\left(x_i,y_j,t_{k+{\displaystyle \frac{1}{2}}}\right)$. For the integral ${\int }_{\Omega }F\left(\left|u\right(x,y,t{|}^2\right)\mathrm{d}x\mathrm{d}y$, it follows from the composite Simpson’s formula, we get

$$\begin{array}{c}\hfill {\int }_{\Omega }F\left({\left|u\left(x_i,y_j,t_{k+{\displaystyle \frac{1}{2}}}\right)\right|}^2\right)\mathrm{d}x\mathrm{d}y={\displaystyle \frac{h_x{h}_y}{36}}\sum _{i=0}^{M_1-1}\sum _{j=0}^{M_2-1}S\left(x_i,y_j,t_{k+{\displaystyle \frac{1}{2}}}\right)+\mathcal{O}\left(h_x^4+h_y^4\right),\end{array}$$

where

$$\begin{array}{cc}\hfill S\left(x_i,y_j,t_{k+{\displaystyle \frac{1}{2}}}\right)=& F\left({\left|u\left(x_i,y_j,t_{k+{\displaystyle \frac{1}{2}}}\right)\right|}^2\right)+F\left({\left|u\left(x_{i+1},y_j,t_{k+{\displaystyle \frac{1}{2}}}\right)\right|}^2\right)+F\left({\left|u\left(x_i,y_{j+1},t_{k+{\displaystyle \frac{1}{2}}}\right)\right|}^2\right)\hfill \\ \hfill & +F\left({\left|u\left(x_{i+1},y_{j+1},t_{k+{\displaystyle \frac{1}{2}}}\right)\right|}^2\right)+4\left[F\left({\left|u\left(x_{i+{\displaystyle \frac{1}{2}}},y_j,t_{k+{\displaystyle \frac{1}{2}}}\right)\right|}^2\right)+F\left({\left|u\left(x_i,y_{j+{\displaystyle \frac{1}{2}}},t_{k+{\displaystyle \frac{1}{2}}}\right)\right|}^2\right)\right.\hfill \\ \hfill & \left.+F\left({\left|u\left(x_{i+1},y_{j+{\displaystyle \frac{1}{2}}},t_{k+{\displaystyle \frac{1}{2}}}\right)\right|}^2\right)+F\left({\left|u\left(x_{i+{\displaystyle \frac{1}{2}}},y_{j+1},t_{k+{\displaystyle \frac{1}{2}}}\right)\right|}^2\right)\right]+16F\left({\left|u\left(x_{i+{\displaystyle \frac{1}{2}}},y_{j+{\displaystyle \frac{1}{2}}},t_{k+{\displaystyle \frac{1}{2}}}\right)\right|}^2\right).\hfill \end{array}$$

By using (5) to spatial discretization, and applying the Crank-Nicolson method in time, then we get

$$\begin{array}{cc}\hfill & \mathrm{i}{\mathcal{H}}_x^{\alpha }{\mathcal{H}}_y^{\beta }{\delta }_{t}u\left(x_i,y_j,t_{k+{\displaystyle \frac{1}{2}}}\right)+\mu \left({\mathcal{H}}_y^{\beta }{\delta }_x^{\alpha }+{\mathcal{H}}_x^{\alpha }{\delta }_y^{\beta }\right)u\left(x_i,y_j,t_{k+{\displaystyle \frac{1}{2}}}\right)+\hfill \\ \hfill & \rho {\mathcal{H}}_x^{\alpha }{\mathcal{H}}_y^{\beta }r\left(t_{k+{\displaystyle \frac{1}{2}}}\right)Q\left(u\left(x_i,y_j,t_{k+{\displaystyle \frac{1}{2}}}\right)\right)=\mathcal{O}\left({\tau }^2+h_x^4+h_y^4\right),\hfill \end{array}$$

(6)

$$\begin{array}{c}\hfill {\delta }_{t}r\left(t_{k+{\displaystyle \frac{1}{2}}}\right)=\Re \left(Q\left(u\left(x_i,y_j,t_{k+{\displaystyle \frac{1}{2}}}\right)\right),{\delta }_{t}u\left(x_i,y_j,t_{k+{\displaystyle \frac{1}{2}}}\right)\right)+\mathcal{O}\left({\tau }^2+h_x^4+h_y^4\right),\end{array}$$

(7)

where

$$\begin{array}{c}\hfill r\left(t_{k+{\displaystyle \frac{1}{2}}}\right)=\sqrt{{\displaystyle \frac{h_x{h}_y}{36}}\sum _{i=0}^{M_1-1}\sum _{j=0}^{M_2-1}S\left(x_i,y_j,t_{k+{\displaystyle \frac{1}{2}}}\right)+C_0},\end{array}$$

and

$$\begin{array}{c}\hfill Q\left(u\left(x_i,y_j,t_{k+{\displaystyle \frac{1}{2}}}\right)\right)={\displaystyle \frac{{\left|u\left(x_i,y_j,t_{k+\frac{1}{2}}\right)\right|}^{2}u\left(x_i,y_j,t_{k+\frac{1}{2}}\right)}{\sqrt{\frac{h_x{h}_y}{36}{\sum }_{i=0}^{M_1-1}{\sum }_{j=0}^{M_2-1}S\left(x_i,y_j,t_{k+\frac{1}{2}}\right)+C_0}}}.\end{array}$$

Omitting the high-order terms in (6) and (7), replacing the exact solutions $u\left(x_i,y_j,t_{k+{\displaystyle \frac{1}{2}}}\right)$, $Q\left(u\left(x_i,y_j,t_{k+{\displaystyle \frac{1}{2}}}\right)\right)$, and $r\left(t_{k+{\displaystyle \frac{1}{2}}}\right)$ with their numerical solutions $u_{ij}^{k+{\displaystyle \frac{1}{2}}}$, $Q\left(u_{ij}^{k+{\displaystyle \frac{1}{2}}}\right)$ and $r^{k+{\displaystyle \frac{1}{2}}}$, we arrive at

$$\begin{array}{c}\hfill \mathrm{i}{\mathcal{H}}_x^{\alpha }{\mathcal{H}}_y^{\beta }{\delta }_t{u}_{ij}^{k+{\displaystyle \frac{1}{2}}}+\mu \left({\mathcal{H}}_y^{\beta }{\delta }_x^{\alpha }+{\mathcal{H}}_x^{\alpha }{\delta }_y^{\beta }\right)u_{ij}^{k+{\displaystyle \frac{1}{2}}}+\rho {\mathcal{H}}_x^{\alpha }{\mathcal{H}}_y^{\beta }{r}^{k+{\displaystyle \frac{1}{2}}}Q\left(u_{ij}^{k+{\displaystyle \frac{1}{2}}}\right)=0,\end{array}$$

(8)

$$\begin{array}{c}\hfill {\delta }_t{r}^{k+{\displaystyle \frac{1}{2}}}=\Re \left(Q\left(u^{k+{\displaystyle \frac{1}{2}}}\right),{\delta }_t{u}^{k+{\displaystyle \frac{1}{2}}}\right).\end{array}$$

(9)

$${\displaystyle u_{ij}^0=u_0(x_i,y_j),(x_i,y_j)\in \partial \Omega ,}$$

3. Theoretical Analysis of the Numerical Scheme

In this section, we study the conservation properties of the schemes (8)-(9). Firstly, we must introduce several lemmas before giving the main results, which are useful in the subsequent analysis.

Lemma 1.

Let ${\mathbf{C}}_p(p=1,2)$ is a symmetric tri-diagonal matrix of $(M_p-1)$-square, denoted as

$${\mathbf{C}}_p={\left(\begin{array}{ccccc}{\displaystyle \frac{2}{3}}& {\displaystyle \frac{1}{6}}& 0& \cdots & 0\\ {\displaystyle \frac{1}{6}}& {\displaystyle \frac{2}{3}}& {\displaystyle \frac{1}{6}}& 0& \cdots \\ ⋮& ⋮& \ddots & \ddots & \ddots \\ 0& \dots & {\displaystyle \frac{1}{6}}& {\displaystyle \frac{2}{3}}& {\displaystyle \frac{1}{6}}\\ 0& 0& \dots & {\displaystyle \frac{1}{6}}& {\displaystyle \frac{2}{3}}\end{array}\right)}_{(M_p-1)\times (M_p-1)},$$

then ${\mathbf{C}}_p$ is positive define.

Proof. 

It is not difficult to know that the eigenvalues of the matrix ${\mathbf{C}}_p$ are

$$\begin{array}{c}{\displaystyle {\lambda }_j\left({\mathbf{C}}_p\right)=\frac{2}{3}+\frac{1}{3}cos\left(\frac{j\pi }{M_p}\right)=\frac{1}{3}+\frac{2}{3}{cos}^2\left(\frac{j\pi }{2M_p}\right)>0,j=1,2,\dots ,M_p-1,}\hfill \end{array}$$

then we can claim that the matrix ${\mathbf{C}}_p$ is positive definite. This ends the proof. □

Lemma 2.

[1] Denote

$${\mathbf{G}}_p^{\left(\gamma \right)}={\displaystyle \frac{1}{h_z^{\gamma }}}{\left(\begin{array}{ccccc}{\kappa }_1^{\left(\gamma \right)}& {\kappa }_0^{\left(\gamma \right)}& 0& \cdots & 0\\ {\kappa }_2^{\left(\gamma \right)}& {\kappa }_1^{\left(\gamma \right)}& {\kappa }_0^{\left(\gamma \right)}& 0& \cdots \\ ⋮& ⋮& \ddots & \ddots & \ddots \\ {\kappa }_{M_p-2}^{\left(\gamma \right)}& \dots & {\kappa }_2^{\left(\gamma \right)}& {\kappa }_1^{\left(\gamma \right)}& {\kappa }_0^{\left(\gamma \right)}\\ {\kappa }_{M_p-1}^{\left(\gamma \right)}& {\kappa }_{M_p-2}^{\left(\gamma \right)}& \dots & {\kappa }_2^{\left(\gamma \right)}& {\kappa }_1^{\left(\gamma \right)}\end{array}\right)}_{(M_p-1)\times (M_p-1)},$$

which is an associate matrix of fractional difference quotient operator ^L${\mathcal{A}}_{h_z}^{\gamma }$. Then the matrix

$$\begin{array}{c}{\displaystyle {\mathbf{D}}_p^{\left(\gamma \right)}=-\frac{1}{2cos\left(\frac{\pi }{2}\gamma \right)}\left({\mathbf{G}}_p^{\left(\gamma \right)}+{{\mathbf{G}}_p^{\left(\gamma \right)}}^T\right),1<\gamma <2,}\hfill \end{array}$$

is the associate matrix of fractional difference quotient operator ${\delta }_z^{\gamma }$ and it is negative semi-definite.

Lemma 3.

For any two grid function u and v, there exists a linear difference operator denoted by ${\delta }^{\alpha ,\beta }$ such that

$$\begin{array}{c}\hfill \left({\mathcal{H}}_y^{\beta }{\delta }_x^{\alpha }+{\mathcal{H}}_x^{\alpha }{\delta }_y^{\beta }u,v\right)=-\left({\delta }^{\alpha ,\beta }u,{\delta }^{\alpha ,\beta }v\right).\end{array}$$

Proof. 

Denote

$$\begin{array}{c}u={(u_{1,1},u_{2,1},\dots ,u_{M_1-1,1},\dots ,u_{1,M_2-1},u_{2,M_2-1},\dots ,u_{M_1-1,M_2-1})}^T\hfill \end{array}$$

and

$$\begin{array}{c}v={(v_{1,1},v_{2,1},\dots ,v_{M_1-1,1},\dots ,v_{1,M_2-1},v_{2,M_2-1},\dots ,v_{M_1-1,M_2-1})}^T.\hfill \end{array}$$

Then we know that there holds that

$$\begin{array}{c}{\displaystyle \left({\mathcal{H}}_y^{\beta }{\delta }_x^{\alpha }+{\mathcal{H}}_x^{\alpha }{\delta }_y^{\beta }u,v\right)=h_1{h}_2{v}^T\mathbf{H}u,}\hfill \end{array}$$

where

$$\begin{array}{c}{\displaystyle \mathbf{H}=\left({\mathbf{C}}_2\otimes {\mathbf{I}}_1\right)\left({\mathbf{I}}_2\otimes {\mathbf{D}}_1^{\left(\alpha \right)}\right)+\left({\mathbf{I}}_2\otimes {\mathbf{C}}_1\right)\left({\mathbf{D}}_2^{\left(\beta \right)}\otimes {\mathbf{I}}_1\right).}\hfill \end{array}$$

Here ${\mathbf{I}}_p$ are the unit matrices of order $M_p-1$$(p=1,2)$, the symbol ⊗ stands for the Kronecker product of any two matrices.

Hence, it follows from the properties of the Kronecker product [9], we get

$$\begin{array}{cc}\hfill {\mathbf{H}}^T=& {\left({\mathbf{I}}_2\otimes {\mathbf{D}}_1^{\left(\alpha \right)}\right)}^T{\left({\mathbf{C}}_2\otimes {\mathbf{I}}_1\right)}^T+{\left({\mathbf{D}}_2^{\left(\beta \right)}\otimes {\mathbf{I}}_1\right)}^T{\left({\mathbf{I}}_2\otimes {\mathbf{C}}_1\right)}^T\hfill \\ \hfill =& \left({\mathbf{I}}_2\otimes {\mathbf{D}}_1^{\left(\alpha \right)}\right)\left({\mathbf{C}}_2\otimes {\mathbf{I}}_1\right)+\left({\mathbf{D}}_2^{\left(\beta \right)}\otimes {\mathbf{I}}_1\right)\left({\mathbf{I}}_2\otimes {\mathbf{C}}_1\right)\hfill \\ \hfill =& \left({\mathbf{C}}_2\otimes {\mathbf{I}}_1\right)\left({\mathbf{I}}_2\otimes {\mathbf{D}}_1\right)+\left({\mathbf{I}}_2\otimes {\mathbf{C}}_1\right)\left({\mathbf{D}}_2\otimes {\mathbf{I}}_1\right)=\mathbf{H},\hfill \end{array}$$

that is to say, the matrix $\mathbf{H}$ is real symmetric. Furthermore, based on the Lemmas 1 and 2, we can also know that matrix $\mathbf{H}$ is semi-negative definite because its eigenvalues are nonpositive. So there exist an orthogonal matrix $\mathbf{L}$ and a diagonal matrix $\Lambda$ such that

$$\begin{array}{c}{\displaystyle \mathbf{H}=-{\mathbf{L}}^T\Lambda \mathbf{L}=-{\left({\Lambda }^{\frac{1}{2}}\mathbf{L}\right)}^T\left({\Lambda }^{\frac{1}{2}}\mathbf{L}\right).}\hfill \end{array}$$

Therefore, we have

$$\begin{array}{c}\left({\mathcal{H}}_y^{\beta }{\delta }_x^{\alpha }+{\mathcal{H}}_x^{\alpha }{\delta }_y^{\beta }u,v\right)=h_1{h}_2{v}^T\mathbf{H}u=-h_1{h}_2{\left({\Lambda }^{{\displaystyle \frac{1}{2}}}\mathbf{L}v\right)}^T{\Lambda }^{{\displaystyle \frac{1}{2}}}\mathbf{L}u=-\left({\delta }^{\alpha ,\beta }u,{\delta }^{\alpha ,\beta }v\right),\hfill \end{array}$$

where ${\Lambda }^{{\displaystyle \frac{1}{2}}}\mathbf{L}$ is the associate matrix of fractional difference quotient operator ${\delta }^{\alpha ,\beta }$, which completes the proof. □

Lemma 4.

For any two grid function u and v, there exists a linear difference operator denoted by ${\mathcal{H}}^{\alpha ,\beta }$ such that

$$\begin{array}{c}\hfill \left({\mathcal{H}}_x^{\alpha }{\mathcal{H}}_y^{\beta }u,v\right)=\left({\mathcal{H}}^{\alpha ,\beta }u,{\mathcal{H}}^{\alpha ,\beta }v\right).\end{array}$$

Proof. 

Using the same method as Lemma 3, we can easily get the results. Therefore, the detailed proof process is omitted. □

Below, we give the main results:

Theorem 1.

The numerical schemes (8)-(9) are mass-conserved, namely

$$\begin{array}{c}\hfill {\mathcal{M}}_{\tau }\left(u^{k+1}\right)={\mathcal{M}}_{\tau }\left(u^k\right),k=0,1,\dots ,N-1,\end{array}$$

(10)

where

$$\begin{array}{c}\hfill {\mathcal{M}}_{\tau }\left(u^k\right)=∥{\mathcal{H}}^{\alpha ,\beta }{u}^k∥\end{array}$$

is the discrete mass of the numerical solution at $t_k$.

Proof. 

Taking the discrete inner product of the equation (8) with $u^{k+{\displaystyle \frac{1}{2}}}$, one has

$$\begin{array}{c}\hfill \mathrm{i}\left({\mathcal{H}}_x^{\alpha }{\mathcal{H}}_y^{\beta }{\delta }_t{u}_{ij}^{k+{\displaystyle \frac{1}{2}}},u^{k+{\displaystyle \frac{1}{2}}}\right)+\mu \left(\left({\mathcal{H}}_y^{\beta }{\delta }_x^{\alpha }+{\mathcal{H}}_x^{\alpha }{\delta }_y^{\beta }\right)u^{k+{\displaystyle \frac{1}{2}}},u^{k+{\displaystyle \frac{1}{2}}}\right)+\rho \left({\mathcal{H}}_x^{\alpha }{\mathcal{H}}_y^{\beta }{r}^{k+{\displaystyle \frac{1}{2}}}Q\left(u^{k+{\displaystyle \frac{1}{2}}}\right),u^{k+{\displaystyle \frac{1}{2}}}\right)=0.\end{array}$$

(11)

Using the Lemmas 3 and 4, and taking the imaginary part of (11) leads to

$$\begin{array}{c}\hfill ∥{\mathcal{H}}^{\alpha ,\beta }{u}^{k+1}∥=∥{\mathcal{H}}^{\alpha ,\beta }{u}^k∥,\end{array}$$

which implies that the (10) holds, the proof is completed. □

Theorem 2.

The numerical schemes (8)-(9) are energy-conserved, that is

$$\begin{array}{c}\hfill {\mathcal{E}}_{\tau }\left(u^{k+1}\right)={\mathcal{E}}_{\tau }\left(u^k\right),k=0,1,\dots ,N-1,\end{array}$$

(12)

where

$$\begin{array}{c}\hfill {\mathcal{E}}_{\tau }\left(u^k\right)=\mu {∥{\delta }^{\alpha ,\beta }{u}^k∥}^2-\rho {\mathcal{H}}_x^{\alpha }{\mathcal{H}}_y^{\beta }{\left|r^k\right|}^2,\end{array}$$

is the discrete energy of the numerical solution at $t_k$.

Proof. 

Taking the discrete inner product of the equation (8) with ${\delta }_t{u}^{k+{\displaystyle \frac{1}{2}}}$, and by follows from the Lemmas 3 and 4 yields

$$\begin{array}{c}\hfill \mathrm{i}∥{\mathcal{H}}^{\alpha ,\beta }{\delta }_t{u}_{ij}^{k+{\displaystyle \frac{1}{2}}}∥-\mu \left({\delta }^{\alpha ,\beta }{u}^{k+{\displaystyle \frac{1}{2}}},{\delta }^{\alpha ,\beta }{\delta }_t{u}^{k+{\displaystyle \frac{1}{2}}}\right)+\rho \left({\mathcal{H}}_x^{\alpha }{\mathcal{H}}_y^{\beta }{r}^{k+{\displaystyle \frac{1}{2}}}Q\left(u^{k+{\displaystyle \frac{1}{2}}}\right),{\delta }_t{u}^{k+{\displaystyle \frac{1}{2}}}\right)=0.\end{array}$$

(13)

Taking the real part of (13) , then one has

$$\begin{array}{c}\hfill {\displaystyle \frac{\mu }{2\tau }}\left[{∥{\delta }^{\alpha ,\beta }{u}^{k+1}∥}^2-{∥{\delta }^{\alpha ,\beta }{u}^k∥}^2\right]=\rho r^{k+{\displaystyle \frac{1}{2}}}\Re \left\{\left({\mathcal{H}}_x^{\alpha }{\mathcal{H}}_y^{\beta }Q\left(u^{k+{\displaystyle \frac{1}{2}}}\right),{\delta }_t{u}^{k+{\displaystyle \frac{1}{2}}}\right)\right\}.\end{array}$$

(14)

Next, if we multiply both sides of equation (9) by ${\mathcal{H}}_x^{\alpha }{\mathcal{H}}_y^{\beta }{r}^{k+{\displaystyle \frac{1}{2}}}$, we get

$$\begin{array}{c}\hfill {\mathcal{H}}_x^{\alpha }{\mathcal{H}}_y^{\beta }{r}^{k+{\displaystyle \frac{1}{2}}}{\delta }_t{r}^{k+{\displaystyle \frac{1}{2}}}=r^{k+{\displaystyle \frac{1}{2}}}\Re \left\{\left({\mathcal{H}}_x^{\alpha }{\mathcal{H}}_y^{\beta }Q\left(u_{ij}^{k+{\displaystyle \frac{1}{2}}}\right),{\delta }_t{u}_{ij}^{k+{\displaystyle \frac{1}{2}}}\right)\right\}.\end{array}$$

(15)

Combine (14) and (15) to get

$$\begin{array}{c}\hfill \mu \left[{∥{\delta }^{\alpha ,\beta }{u}^{k+1}∥}^2-{∥{\delta }^{\alpha ,\beta }{u}^k∥}^2\right]=\rho {\mathcal{H}}_x^{\alpha }{\mathcal{H}}_y^{\beta }\left[{\left|r^{k+1}\right|}^2-{\left|r^k\right|}^2\right],\end{array}$$

in other words, (12) is true, and that completes the proof. □

Remark 1.

In fact, the proposed difference scheme should be both convergent with second-order accuracy in time and fourth-order accuracy in spaces, respectively. However, the detailed convergence analysis and numerical tests will rely on our future studies due to space constraints.