\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{multirow} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2017 (2017), No. 166, pp. 1--18.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2017 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2017/166\hfil Finite element method for fractional NLS] {Finite element method for time-space-fractional Schr\"odinger equation} \author[X. G. Zhu, Z. B. Yuan, J. G. Wang, Y. F. Nie, Z. Z. Yang \hfil EJDE-2017/166\hfilneg] {Xiaogang Zhu, Zhanbin Yuan, Jungang Wang, Yufeng Nie, Zongze Yang } \address{Xiaogang Zhu \newline Department of Applied Mathematics, Northwestern Polytechnical University, Xi'an 710129, China} \email{zhuxg590@yeah.net} \address{Zhanbin Yuan \newline Department of Applied Mathematics, Northwestern Polytechnical University, Xi'an 710129, China} \email{yzzzb@nwpu.edu.cn} \address{Jungang Wang \newline Department of Applied Mathematics, Northwestern Polytechnical University, Xi'an 710129, China} \email{wangjungang@nwpu.edu.cn} \address{Yufeng Nie \newline Department of Applied Mathematics, Northwestern Polytechnical University, Xi'an 710129, China} \email{yfnie@nwpu.edu.cn} \address{Zongze Yang \newline Department of Applied Mathematics, Northwestern Polytechnical University, Xi'an 710129, China} \email{yangzongze@gmail.com} \dedicatory{Communicated by Mokhtar Kirane} \thanks{Submitted January 30, 2016. Published July 5, 2017.} \subjclass[2010]{35R11, 65M60, 65M12} \keywords{Time-space-fractional NLS; finite element method; convergence} \begin{abstract} In this article, we develop a fully discrete finite element method for the nonlinear Schr\"odinger equation (NLS) with time- and space-fractional derivatives. The time-fractional derivative is described in Caputo's sense and the space-fractional derivative in Riesz's sense. Its stability is well derived; the convergent estimate is discussed by an orthogonal operator. We also extend the method to the two-dimensional time-space-fractional NLS and to avoid the iterative solvers at each time step, a linearized scheme is further conducted. Several numerical examples are implemented finally, which confirm the theoretical results as well as illustrate the accuracy of our methods. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{remark}[theorem]{Remark} \newtheorem{definition}[theorem]{Definition} \newtheorem{example}[theorem]{Example} \allowdisplaybreaks \section{Introduction}\label{intro} The fractional calculus, as a generalization of classical integer calculus, possesses a long history and affluent connotation, which was discovered by mathematicians over three hundred years ago. Recently, due to their excellent properties to characterize the effects of memory and long-range interaction, fractional partial differential equations (FPDEs) have aroused keen interests among academic circles and have also been applied broadly in various applications, examples including complex network, stochastic interfaces, synoptic climatology, certain option pricing mechanism, medical image processing, dielectric polarization, and the chaotic dynamics of nonlinear systems, etc. In virtue of the singular integral form of fractional derivatives, however, solving FPDEs suffers more obstacles than those associated with classical derivatives. In many cases, as we all know, the analytic solutions are unattainable or even unrealistic for most of the mathematical models, so resorting to the efficient numerical approaches to obtain numerical solutions turns into a preferred option. In the past decades, much to our delight, great efforts have been devoted to this area and numerous commonly used methods have been developed, such as finite difference method \cite{Ref004,Ref003,Ref001,Ref005}, finite element method \cite{Ref009,Ref007}, spectral method \cite{Ref006}, adomian decomposition method \cite{Ref0012}, and variational iteration method \cite{Ref0013}. In general, finite difference method and finite element method are the most accepted methods for solving FPDEs. It is noteworthy that, because of the universal mutuality of these models, considering the high-dimensional space-FPDEs efficiently in the numerical aspect appears somewhat challenging \cite{Ref0017,Ref0015,Ref0042,Reh05,Ref0040}. In this work, our aim is to propose a fully discrete finite element method for the one- and two-dimensional time-space-fractional NLS. Due to the similarity, we first focus on the following one-dimensional model \begin{align}\label{eq001} \mathrm{i}\frac{\partial^\alpha \psi}{\partial t^\alpha} +\frac{\partial^{2\beta} \psi}{\partial|x|^{2\beta}} +\lambda f(|\psi|^2)\psi &=0, \quad (x;t) \in\Omega\times(0,T], \end{align} with $\mathrm{i}^2=-1$, real parameters $\lambda$, $0<\alpha\leq 1$, $1/2<\beta\leq 1$, and the initial and boundary conditions given by \begin{gather} \psi(x,0)=\varphi(x), \quad x\in \Omega\cup\partial\Omega,\label{e008}\\ \psi(a,t)=0, \quad \psi(b,t)=0,\quad t \in(0,T], \label{e009} \end{gather} where $\Omega=(a,b)$, $f(s)$ is real on $\mathbb{R}$ and $\varphi(x)$ is a prescribed function. In \eqref{eq001}, the time-fractional derivative is defined in Caputo sense as follow \begin{equation} \label{yy02} \frac{\partial^\alpha \psi(x,t)}{\partial t^\alpha} = {^C_0}D^\alpha_t\psi(x,t) =\frac{1}{\Gamma(1-\alpha)} \int^t_0\frac{\partial \psi(x,\xi)} {\partial \xi}\frac{d\xi}{(t-\xi)^\alpha}, \end{equation} and the Riesz space-fractional derivative is defined as \begin{gather*} \frac{\partial^{2\beta} \psi(x,t)}{\partial|x|^{2\beta}}=\frac{-1}{2\cos(\beta \pi)} \big(D^{2\beta}_L \psi(x,t)+D^{2\beta}_R \psi(x,t)\big),\\ D^{2\beta}_L \psi(x,t) = \frac{1}{\Gamma(2-{2\beta})}\frac{\partial^2}{\partial x^2} \int^x_a(x-\xi)^{1-{2\beta}}\psi(\xi,t)d\xi,\\ D^{2\beta}_R \psi(x,t) = \frac{1}{\Gamma(2-{2\beta})}\frac{\partial^2}{\partial x^2} \int^b_x(\xi-x)^{1-{2\beta}}\psi(\xi,t)d\xi. \end{gather*} In particular, we note \eqref{eq001} degenerates to the classical NLS, while $\alpha,\beta=1$, which has been the subject of intense research in the past few decades \cite{Ref0019,Ref0021,Ref0023,Ref0022}. The time-fractional NLS was generated by Naber \cite{Reh07}, where two kinds of different means were utilized to perform this generalization. The time-fractional NLS, also covering the mixed time-space-fractional NLS, have been investigated by many authors \cite{Reh09,Reh15,Reh16}. For the numerical algorithms to the pure time-fractional NLS, see \cite{Reh13,Reh12,Reh10} for reference. The space-fractional NLS was first raised by Laskin \cite{Ref0024,Ref0025}, via extending the Feynman path integral over the Brownian paths to a new path integral over the L\'evy quantum mechanical paths. Its local and global well-posedness were studied in \cite{Ref0026,Ref0027}. In addition, numerous works are dedicated to researching the theoretical properties for such equation in various regimes \cite{Ref0043,Ref0038,Ref0028,Ref0037}. As to the numerical works for the space- and time-space-fractional NLS in one-dimension, Amore et al. used an effective collocation method to solve the space-fractional NLS \cite{Ref0029}. Herzallah and Gepreel advised an adomian decomposition method for the time-space-fractional NLS \cite{Ref0033}. The variable-order space-fractional NLS was considered in \cite{Ref0032}, using the Crank-Nicholson scheme. The ground state solutions of the semiclassical fractional NLS was investigated in \cite{Ref0046} by a Fourier spectral method. Wang and Huang gave a second-order energy conservative finite difference method for the space-fractional NLS \cite{Ref0035}. Wang et al. proposed two mass conservative difference methods for the space-fractional coupled nonlinear Schr\"odinger equations (CNLS) \cite{Ref0031,Ref0098}. Liu et al. presented the implicit and explicit-implicit finite difference methods for the time-space-fractional NLS and an implicit method for the time-space-fractional CNLS \cite{Ref0096}. In two-dimensional case, Zhao et al. derived a new fourth-order compact operator and applied it to construct a compact alternating direction implicit method for the Riesz space-fractional NLS \cite{Ref0039}, which is linearized and validated to be stable and well convergent. In this context, inspired by these methods in existence, we tend to propose the finite element method to solve the one- and two-dimensional time-space-fractional NLS, and as the character of the space-fractional NLS, some necessary properties of fractional derivatives will be introduced for auxiliary analysis. The layout of this article is organized as follows. In Section \ref{Au}, we recall some fractional derivative spaces and basic properties for the fractional derivatives. In Section \ref{VF}, the weak formulations are described. The fully discrete finite element method is formulated in Section \ref{FD}, and its unconditional stability is analyzed. The convergence of it is discussed in Section \ref{CA}. In Section \ref{NE}, we extend the method to the two-dimensional time-space-fractional NLS and carry out some numerical tests in Section \ref{NY}. \section{Preliminaries} \label{Au} In this section, we present some auxiliary results related to the Riemann-Liouville fractional derivatives. Let $u$, $v$ be the real-valued functions. We employ \[ \|u\|_{L^2(\mathbb{R})}=\langle u,u\rangle^{1/2}_{L^2(\mathbb{R})}, \quad \langle u,v\rangle_{L^2(\mathbb{R})}=\int_{\mathbb{R}}uv dx, \] and denote the norms $\|\cdot\|_{L^2(\Omega)}$ by $\|\cdot\|_0$, $\|\cdot\|_{H^c(\Omega)}$ by $\|\cdot\|_c$, and the inner products $\langle\cdot,\cdot\rangle_{L^2(\Omega)}$, $\langle\cdot,\cdot\rangle_{L^2(\mathbb{R})}$ by $\langle\cdot,\cdot\rangle$, with $c>0$ and a subinterval $\Omega\subset\mathbb{R}$. Referring to \cite{Ref007}, we define the following fractional derivative spaces by the above inner product. \begin{definition}[Left fractional derivative space] \rm For $\mu>0$, we define the left semi-norm \begin{equation*} |u|_{J^\mu_L(\mathbb{R})}=\|D^\mu_Lu\|_{L^2(\mathbb{R})}, \end{equation*} and left norm \begin{equation*} \|u\|_{J^\mu_L(\mathbb{R})}= \big( \|u\|^2_{L^2(\mathbb{R})} +|u|^2_{J^\mu_L(\mathbb{R})} \big)^{1/2} , \end{equation*} and let $J^\mu_L(\mathbb{R})$ be the closure of $C^\infty(\mathbb{R})$ with respect to $\|\cdot\|_{J^\mu_L(\mathbb{R})}$. \end{definition} \begin{definition}[Right fractional derivative space] \rm For $\mu>0$, we define the right semi-norm \begin{equation*} |u|_{J^\mu_R(\mathbb{R})}=\|D^\mu_Ru\|_{L^2(\mathbb{R})}, \end{equation*} and right norm \begin{equation*} \|u\|_{J^\mu_R(\mathbb{R})}=\big( \|u\|^2_{L^2(\mathbb{R})} +|u|^2_{J^\mu_R(\mathbb{R})} \big)^{1/2} , \end{equation*} and let $J^\mu_R(\mathbb{R})$ be the closure of $C^\infty(\mathbb{R})$ with respect to $\|\cdot\|_{J^\mu_R(\mathbb{R})}$. \end{definition} \begin{definition}[Fractional Sobolev space] \label{le002} \rm Let $\mu>0$ and $\mathcal{F}(u)$ be the Fourier transform of a prescribed $u(x)$ defined on $\mathbb{R}$, i.e., $\mathcal{F}(u)=\int_{\mathbb{R}}u(x)e^{-i\omega x}dx$, with the variable $\omega$. Then, we can define the semi-norm \begin{equation*} |u|_{H^\mu(\mathbb{R})}=\|\,|\omega|^\mu \mathcal{F}(u)\|_{L^2(\mathbb{R}_\omega)}, \end{equation*} and norm \begin{equation*} \|u\|_{H^\mu(\mathbb{R})}=\big(\|u\|^2_{L^2(\mathbb{R})}+|u|^2_{H^\mu(\mathbb{R})}\big)^{1/2}, \end{equation*} and denote the closure of $C^\infty(\mathbb{R})$ with respect to $\|\cdot\|_{H^\mu(\mathbb{R})}$ by $H^\mu(\mathbb{R})$. \end{definition} \begin{lemma}[\cite{Ref007}] For $\mu>0$, $J^\mu_L(\mathbb{R})$, $J^\mu_R(\mathbb{R})$, and $H^\mu(\mathbb{R})$ are equivalent with the equivalent semi-norms and norms. \end{lemma} \begin{lemma}[\cite{Ref007}]\label{le003} For $\mu>0$, we have the property in $L^2$-sense \begin{align}\label{s28} \langle D{^{\mu}_L}u,D{^{\mu}_R}u\rangle=\cos(\mu\pi)\|D{^{\mu}_L}u\|^2_{L^2(\mathbb{R})}. \end{align} \end{lemma} Let $J^\mu_{L,0}(\Omega)$, $J^\mu_{R,0}(\Omega)$, $H^\mu_0(\Omega)$ be the closures of $C^\infty_0(\Omega)$ with respect to $\|\cdot\|_{J^\mu_L(\Omega)}$, $\|\cdot\|_{J^\mu_R(\Omega)}$, and $\|\cdot\|_{H^\mu(\Omega)}$, respectively. Then, we have the following lemmas. \begin{lemma}[\cite{Ref0010}] \label{le009} If $0<\beta<1$, $u\in J^{2\beta}_{L,0}(\Omega)$, and $v\in J^{2\beta}_{R,0}(\Omega)$, then \begin{equation}\label{s29} \langle D{^{2\beta}_L}u,v\rangle=\langle D{^{\beta}_L}u,D{^{\beta}_R}v\rangle, \quad \langle D{^{2\beta}_R}u,v\rangle=\langle D{^{\beta}_R}u,D{^{\beta}_L}v\rangle. \end{equation} \end{lemma} \begin{lemma}[\cite{Ref007}]\label{le007} Let $\mu > 0$, $\mu\neq n-1/2$, $n\in \mathbb{N}$. $J^\mu_{L,0}(\Omega)$, $J^\mu_{R,0}(\Omega)$, and $H^\mu_0(\Omega)$ are equivalent with the equivalent semi-norms and norms. \end{lemma} \begin{lemma}[\cite{Ref007}] \label{le006} For $\mu>0$, $u\in J^\mu_{L,0}(\Omega)$, and $0< \gamma< \mu$, we have \begin{equation}\label{s30} \|u\|_0\leq C|u|_{J^\mu_{L}(\Omega)},\quad \ |u|_{J^\gamma_{L}(\Omega)}\leq C|u|_{J^\mu_{L}(\Omega)}, \end{equation} and for $u\in J^\mu_{R,0}(\Omega)$, $0<\gamma< \mu$, we have \begin{equation}\label{s32} \|u\|_0\leq C|u|_{J^\mu_{R}(\Omega)},\quad \ |u|_{J^\gamma_{R}(\Omega)}\leq C|u|_{J^\mu_{R}(\Omega)}. \end{equation} If $u\in H^\mu_0(\Omega)$, $0<\gamma<\mu$, $\mu\neq n-1/2$, $n\in \mathbb{N}$, the analogous result is obtained. \end{lemma} \begin{remark}\label{re002} \rm In two-dimensions, the fractional derivative and Sobolev spaces can also be established and \eqref{s28}-\eqref{s32} still work; see \cite{Ref0017,Reh02} for overall views. \end{remark} \section{Weak problem}\label{VF} At first, decompose the unknown $\psi(x,t)$ into its real and imaginary parts by \begin{equation*} \psi(x,t)=u(x,t)+\mathrm{i}v(x,t). \end{equation*} Inserting it into \eqref{eq001}-\eqref{e009}, the original problem can be recast as a coupled system \begin{gather*} \frac{\partial^\alpha u}{\partial t^\alpha}+\frac{\partial^{2\beta} v} {\partial|x|^{2\beta}} +\lambda f(u^2+v^2)v=0,\\ \frac{\partial^\alpha v}{\partial t^\alpha}-\frac{\partial^{2\beta} u}{\partial|x|^{2\beta}} -\lambda f(u^2+v^2)u=0, \ \ (x;t) \in\Omega\times(0,T], \end{gather*} with the initial and boundary values \begin{align*} & u(x,0)=\operatorname{Re}\varphi(x), \quad v(x,0)=\operatorname{Im}\varphi(x), \quad x\in \Omega\cup\partial\Omega,\\ & u(x,t)=0, \quad v(x,t)=0, \quad (x;t)\in\partial\Omega\times(0,T]. \end{align*} where ``Re", ``Im" mean retaining the real and imaginary parts, respectively. Then, using Lemma \ref{le009}, we can derive the weak problem, i.e., seek $u(\cdot,t)$, $v(\cdot,t)\in H^{\beta}_0(\Omega)$, for any $\chi_1$, $\chi_2\in H^{\beta}_0(\Omega)$, to solve \begin{gather} \big\langle\frac{\partial^\alpha u}{\partial t^\alpha},\chi_1\big\rangle -\Lambda(v,\chi_1) +\lambda\langle f(u^2+v^2)v,\chi_1\rangle=0, \label{s05} \\ \big\langle\frac{\partial^\alpha v}{\partial t^\alpha},\chi_2\big\rangle +\Lambda(u,\chi_2) -\lambda\langle f(u^2+v^2)u,\chi_2\rangle=0, \label{s06}\\ u(x,0)=\operatorname{Re}\varphi(x), \quad v(x,0)=\operatorname{Im}\varphi(x), \label{s07} \end{gather} with zero boundary values, where $\Lambda(\cdot,\cdot)$ takes the form \[ \Lambda(u,v)=\frac{1}{2\cos(\beta\pi)}\langle D^{\beta}_L u,D^{\beta}_R v\rangle +\frac{1}{2\cos(\beta\pi)}\langle D^{\beta}_R u,D^{\beta}_L v\rangle, \] and $\varphi(x)$ denotes the same initial function prescribed before. \begin{theorem}\label{Th3} The bilinear form $\Lambda(\cdot, \cdot)$ is symmetric and enjoys the properties $\Lambda(u,v)\leq C_1\|u\|_{\beta}\|v\|_{\beta}$, $\Lambda(u,u)\geq C_2\|u\|^2_{\beta}$, where $C_1$, $C_2$ are positive constants. \end{theorem} The above theorem was proved in \cite{Ref0010}, by recalling Lemmas \ref{le003}, \ref{le007}, and \ref{le006}. \section{Fully discrete finite element method}\label{FD} Let $t_n=n \tau$, $n=0,1,\ldots,N$, and $T=\tau N $ with a constant $N\in\mathbb{Z}^+$; we discretize the Caputo derivative by a difference approach as follows \begin{align*} {^C_0}D^\alpha_tw(x,t_{n}) &=\frac{1}{\Gamma(1-\alpha)}\sum^{n}_{j=1}\int^{t_{j}}_{t_{j-1}} \frac{\partial w(x,\xi)}{\partial \xi} \frac{d\xi}{(t_n-\xi)^\alpha} \\ &=\frac{1}{\Gamma(1-\alpha)}\sum^{n}_{j=1}\frac{\partial w(x,t_{j-1/2})}{\partial t} \int^{t_{j}}_{t_{j-1}}\frac{d\xi}{(t_n-\xi)^\alpha}+R^{1}_\tau\\ &=\frac{\tau^{1-\alpha}}{\Gamma(2-\alpha)}\sum^{n-1}_{j=0}b_j\frac{\partial w(x,t_{n-j-1/2})}{\partial t} +R^{1}_\tau \\ &=\frac{\tau^{-\alpha}}{\Gamma(2-\alpha)}\sum^{n-1}_{j=0}b_jd_tw(x,t_{n-j}) +R^{1}_\tau + R^{2}_\tau, \end{align*} where $b_j=(j+1)^{1-\alpha}-j^{1-\alpha}$, $d_tw(x,t_{n-j})=w(x,t_{n-j})-w(x,t_{n-j-1})$, and $j=0,1,\ldots,n-1$. Particularly, we assign $0^0=0$ when $\alpha=1$. According to \cite{Re01}, the truncated errors $R^{1}_\tau$, $R^{2}_\tau$ satisfy $R^{2}_\tau= O(\tau^2)$ and \begin{align*} |R^{1}_\tau|&= \Big|\frac{1}{\Gamma(1-\alpha)}\sum^{n}_{j=1}\int^{t_{j}}_{t_{j-1}} \Big\{\frac{\partial w(x,\xi)} {\partial \xi}-\frac{\partial w(x,t_{j-1/2})}{\partial t}\Big\} \frac{d\xi}{(t_n-\xi)^\alpha}\Big| \\ &\leq C\tau^{2-\alpha}\max _{0\leq t\leq t_n}\Big| \frac{\partial^2w(x,t)}{\partial t^2}\Big|, \end{align*} with a bounded constant $C$ independent of $\tau$ for all $\alpha\in(0,1)$ and $n\geq 1$. If we denote the discretized fractional operator by \begin{equation*} \mathcal{D}_t^\alpha w(x,t_n)=\frac{1}{\Gamma(2-\alpha)}\sum^{n-1}_{j=0} b_j\frac{w(x,t_{n-j})-w(x,t_{n-j-1})}{\tau^{\alpha}}, \end{equation*} then $\frac{\partial^\alpha w(x,t_n)}{\partial t^\alpha}$ can be approximated by $\mathcal{D}_t^\alpha w(x,t_n)$, more precisely, \begin{align}\label{eee028} \frac{\partial^\alpha w(x,t_n)}{\partial t^\alpha} =\mathcal{D}_t^\alpha w(x,t_n) + R_\tau, \end{align} where $R_\tau=R^{1}_\tau+R^{2}_\tau$. Let $\Sigma_h$ be a family of subdivisions of $\Omega$, $h$ be their grid parameters, and $X_h\subset H^\beta_0(\Omega)$ be the finite element subspace, denoted as \begin{equation*} X_h=\{v_h\in H^\beta_0(\Omega)\cap C^0(\bar{\Omega}): v_h|_D\in P_{r}(D), \; \forall D\in \Sigma_h \}, \end{equation*} in which, $P_{r}(D)$ is the set of polynomials of degree at most $r$ on $D$. Discretizing the Caputo derivative by \eqref{eee028} and applying finite element method in space, we obtain \begin{gather} \langle\mathcal{D}_{t}^\alpha U^n,\chi_1\rangle-\Lambda(V^n,\chi_1)+\lambda\langle f((U^n)^2+(V^n)^2)V^n,\chi_1\rangle=0, \label{s08} \\ \langle\mathcal{D}_{t}^\alpha V^n,\chi_2\rangle+\Lambda(U^n,\chi_2)-\lambda\langle f((U^n)^2+(V^n)^2)U^n,\chi_2\rangle=0,\label{s09} \end{gather} with any $\chi_1$, $\chi_2\in X_h$, $U^0=\operatorname{Re}\delta\varphi(x)$, $V^0=\operatorname{Im}\delta\varphi(x)$, and $\mathcal{D}^\alpha_{t}w^n$, given by \begin{equation*} \mathcal{D}^\alpha_{t}w^n= \frac{1}{\Gamma(2-\alpha)} \sum^{n-1}_{j=0}b_j\frac{w^{n-j}-w^{n-j-1}}{\tau^{\alpha}}, \end{equation*} where $w^n=U^n$ or $V^n$, $\delta\varphi(x)$ is vital to the starter of \eqref{s08}-\eqref{s09} and thereby shall be properly chosen. Let us introduce the parameter $G_\alpha$: \begin{equation*} G_\alpha:=\tau^\alpha\Gamma(2-\alpha). \end{equation*} Regroup $\mathcal{D}^\alpha_{t}w^n$ as follows \begin{align*} \mathcal{D}_{t}^\alpha w^n &=\frac{\tau^{-\alpha} }{\Gamma(2-\alpha)}\sum^{n-1}_{j=0}b_j(w^{n-j}-w^{n-j-1})\\ &=G_\alpha^{-1}\Big(w^{n}-\sum^{n-1}_{j=1}(b_{j-1}-b_{j})w^{n-j}- b_{n-1}w^0\Big), \end{align*} then \eqref{s08}-\eqref{s09} can be recast, i.e., seek $U^{n}$, $V^{n}\in X_h$, for $\chi_1$, $\chi_2\in X_h$, such that \begin{gather} \begin{aligned} &\langle U^n, \chi_1\rangle-G_\alpha\Lambda(V^n,\chi_1) \\ &= \sum^{n-1}_{j=1} (b_{j-1}-b_{j})\langle U^{n-j},\chi_1\rangle + b_{n-1}\langle U^0,\chi_1 \rangle-\lambda G_\alpha\langle f((U^n)^2+(V^n)^2)V^{n},\chi_1\rangle, \end{aligned} \label{eh003} \\ \begin{aligned} &\langle V^n, \chi_2\rangle+G_\alpha\Lambda(U^n,\chi_2) \\ &= \sum^{n-1}_{j=1} (b_{j-1}-b_{j})\langle V^{n-j},\chi_2\rangle + b_{n-1}\langle V^0,\chi_2 \rangle+\lambda G_\alpha\langle f((U^n)^2+(V^n)^2)U^{n},\chi_2\rangle, \end{aligned} \label{s10} \end{gather} with the initial conditions \begin{equation}\label{e006} U^0=\operatorname{Re}\delta\varphi(x),\quad V^0=\operatorname{Im}\delta\varphi(x), \quad x\in \Omega\cup\partial\Omega, \end{equation} and boundary conditions \begin{equation}\label{e007} U^j=0, \quad V^j=0, \quad 1\leq j\leq N, \quad \text{on } \partial\Omega. \end{equation} In particular, as $n=1$, i.e., at the first step, \eqref{eh003}-\eqref{s10} simply become \begin{gather} \langle U^1, \chi_1\rangle-G_\alpha\Lambda(V^1,\chi_1) = \langle U^0,\chi_1\rangle -\lambda G_\alpha\langle f((U^1)^2+(V^1)^2)V^{1},\chi_1\rangle,\label{e005}\\ \langle V^1, \chi_2\rangle+G_\alpha\Lambda(U^1,\chi_2)= \langle V^0,\chi_2\rangle +\lambda G_\alpha\langle f((U^1)^2+(V^1)^2)U^{1},\chi_2\rangle.\label{s12} \end{gather} Next, we proceed with the stability for \eqref{eh003}-\eqref{e007}. To this end, we recall some requisite properties for $b_j=(j+1)^{1-\alpha}-j^{1-\alpha}$, $j=0,1,\ldots,n$, which read \begin{gather} 1=b_0>b_1>b_2> \cdots >b_n>0,\quad b_n\to 0,\quad \text{ as } n\to +\infty,\label{t01}\\ \sum^{n}_{j=1}(b_{j-1}-b_{j})+b_{n}=(1-b_1)+\sum^{n-1}_{j=2}(b_{j-1}-b_{j})+b_{n-1} =1.\label{t02} \end{gather} \begin{theorem}\label{Th5} The method \eqref{eh003}-\eqref{e007} is stable in the sense that \begin{equation*} \|U^{n}\|^2_{0}+\|V^{n}\|^2_{0}\leq \|U^{0}\|^2_{0}+\|V^{0}\|^2_{0},\quad n=1,2,\ldots,N. \end{equation*} \end{theorem} \begin{proof} Let $\chi_1=U^n$, $\chi_2=V^n$ in \eqref{eh003}-\eqref{s10}. It follows that \begin{gather} \begin{aligned} \langle U^n, U^n\rangle-G_\alpha\Lambda(V^n,U^n) &= \sum^{n-1}_{j=1}(b_{j-1}-b_{j}) \langle U^{n-j},U^n\rangle + b_{n-1}\langle U^0,U^n\rangle\\ &\quad -\lambda G_\alpha \langle f((U^n)^2+(V^n)^2)V^{n},U^n\rangle, \end{aligned} \label{s13} \\ \begin{aligned} \langle V^n, V^n\rangle+G_\alpha\Lambda(U^n,V^n) &=\sum^{n-1}_{j=1}(b_{j-1}-b_{j}) \langle V^{n-j},V^n\rangle + b_{n-1}\langle V^0,V^n\rangle \\ &\quad +\lambda G_\alpha \langle f((U^n)^2+(V^n)^2)U^{n},V^n\rangle. \end{aligned} \label{s15} \end{gather} The sum of \eqref{s13} and \eqref{s15} shows \begin{align*} &\|U^n\|_0^2 + \|V^n\|_0^2 \\ &= \sum^{n-1}_{j=1}( b_{j-1} - b_{j}) (\langle U^{n-j}, U^n\rangle + \langle V^{n-j},V^n\rangle) + b_{n-1}(\langle U^0, U^n\rangle + \langle V^0, V^n\rangle), \end{align*} since $\Lambda(U^n,V^n)=\Lambda(V^n,U^n)$. We go on with the work by the method of induction. In the case of $n=1$, using H\"older's and Young's inequalities, there holds \begin{align*} \|U^1\|_0^2+\|V^1\|_0^2 &=\langle U^0,U^1\rangle+\langle V^0,V^1\rangle \\ &\leq \frac{1}{2}(\|U^0\|_0^2+\|V^0\|_0^2)+\frac{1}{2}(\|U^1\|_0^2+\|V^1\|_0^2), \end{align*} which suggests the result at the first step, i.e., $\|U^1\|_0^2+\|V^1\|_0^2\leq\|U^0\|_0^2+\|V^0\|_0^2$. We now preset the essential hypothesis \begin{equation}\label{s16} \|U^n\|_0^2+\|V^n\|_0^2\leq\|U^0\|_0^2+\|V^0\|_0^2, \quad n=2,3,\ldots,p-1, \end{equation} and begin to consider the case of $n=p$. Using \eqref{s16}, we obtain \begin{align*} \|U^p\|_0^2+\|V^p\|_0^2 &=\sum^{p-1}_{j=1}( b_{j-1}-b_{j})(\|U^{p-j}\|_0 \|U^p\|_0+\|V^{p-j}\|_0\|V^p\|_0) \\ &\quad+b_{p-1}(\|U^0\|_0\|U^p\|_0+\|V^0\|_0\|V^p\|_0)\\ &\leq\frac{1}{2}\sum^{p-1}_{j=1}( b_{j-1} - b_{j})(\|U^{p-j}\|_0^2+\|V^{p-j}\|_0^2) +\frac{1}{2}b_{p-1}(\|U^0\|_0^2+\|V^0\|_0^2) \\ &\quad+\frac{1}{2}\sum^{p-1}_{j=1}( b_{j-1}-b_{j})(\|U^{p}\|_0^2+\|V^{p}\|_0^2) +\frac{1}{2}b_{p-1}(\|U^p\|_0^2+\|V^p\|_0^2). \end{align*} Then, via \eqref{t01}-\eqref{t02}, we easily see that \begin{align*} \|U^p\|_0^2+\|V^p\|_0^2 &\leq \sum^{p-1}_{j=1}( b_{j-1}-b_{j})(\|U^{p-j}\|_0^2 +\|V^{p-j}\|_0^2) +b_{p-1}(\|U^0\|_0^2+\|V^0\|_0^2)\\ &\leq \Big((1-b_1) +\sum^{p-1}_{j=2}(b_{j-1}-b_{j})+b_{p-1}\Big) (\|U^0\|_0^2+\|V^0\|_0^2)\\ &= \|U^0\|_0^2+\|V^0\|_0^2, \end{align*} and hence, the method is unconditionally stable, which proves the theorem. \end{proof} \section{Convergent analysis}\label{CA} In this part, we describe the error estimate for \eqref{eh003}-\eqref{e007}. Define a Ritz projection $\Pi_h:H_0^\beta(\Omega)\mapsto X_h$, via the orthogonal relation \begin{equation}\label{eq0068} \Lambda(u-\Pi_hu,v_h)= 0, \quad \forall v_h\in X_h. \end{equation} Then, relying on \cite{Ref0018,Ref009}, for a $\epsilon\in(0,1/2)$, the following lemma is admitted. \begin{lemma}\label{le008} If $u\in H_0^\beta(\Omega)\cap H^{r+1}(\Omega)$ and $\epsilon\in(0,1/2)$, then \[ \|u-\Pi_h u\|_{0}\leq Ch^{\tilde{r}+1}\|u\|_{r+1}, \] where if $\beta\neq3/4$, $\tilde{r}=r$ and if $\beta=3/4$, $\tilde{r}=r-\epsilon$; $C$ is independent of $h$. \end{lemma} \begin{lemma}\label{le068} Let $\varepsilon^j\geq0$, $R\geq0$, $j=0,1,\ldots,N$, and satisfy \begin{equation*}\label{eee06} \varepsilon^n\leq \sum^{n-1}_{j=1}(b_{j-1}-b_j)\varepsilon^{n-j}+R. \end{equation*} Then, when $0<\alpha<1$, we have \begin{equation}\label{eee01} \varepsilon^n\leq b^{-1}_{n-1} R\leq n^\alpha R/(1-\alpha), \end{equation} and when $\alpha$ closes to $1$, it turns to be \begin{equation}\label{eee02} \varepsilon^n\leq nR. \end{equation} \end{lemma} \begin{proof} In fact, the lemma is implicitly involved in \cite{Ref009,Re01,Reh06}; here, we only underline the process to \eqref{eee02}, since it can not be directly derived from \eqref{eee01} when $\alpha\to 1$. It holds trivially as $n=1$. We now suppose that \begin{equation*} \varepsilon^n\leq nR, \quad n=2,3,\ldots, p-1, \end{equation*} and show the result remains valid at $n=p$. By \eqref{t01}-\eqref{t02}, it is clear that \begin{align*} \varepsilon^p &\leq\sum^{p-1}_{j=1}(b_{j-1}-b_j)\varepsilon^{p-j} +R\leq\sum^{p-1}_{j=1}(b_{j-1}-b_j)(p-j)R+R\\ & \leq \sum^{p-1}_{j=1}(b_{j-1}-b_j)(p-1)R+R \leq pR, \end{align*} in which, $1+ \sum^{p-1}_{j=1}(b_{j-1}-b_{j})(p-1)\leq p$. As a result, \eqref{eee02} is concluded. \end{proof} In the sequel we use the following notation: \[ \partial_tu=\partial u/\partial t, \quad \partial^2_{t}u=\partial^2 u/\partial t^2, \quad \partial^\alpha_{t}u={^C_0}D^\alpha_tu, \] and as stressed before, we select $\delta \varphi(x)=\Pi_h \varphi(x)$ and $\operatorname{Re}\varphi(x)$, $\operatorname{Im}\varphi(x)\in H^{r+1}(\Omega)$ so that the discrete system can be started; $C$ will be regarded as a general constant that may be different at different occasions. In the error analysis, we let $\psi^n$ be the analytic solution to the model \eqref{eq001}-\eqref{e009} at $t=t_n$ with its real and imaginary parts $u^n$, $v^n$, and $\{U^{n}\}_{n=0}^{N}$, $\{V^{n}\}_{n=0}^{N}$ be the numerical solutions obtained by \eqref{eh003}-\eqref{e007}. Also, we define $\psi^n_h=U^n+\mathrm{i}V^n$ and the complex norm \begin{equation}\label{s23} \|\psi^n-\psi^n_h\|_0=\big(\|u^n-U^n\|^2_0+\|v^n-V^n\|^2_0\big)^{1/2}, \quad n=0,1,\ldots,N. \end{equation} Then, with the help of mathematical induction, we state the convergent theorem. \begin{theorem}\label{Th6} Assume $u, v$, $\partial^\alpha_t u, \partial^\alpha_t v$, $\partial_t u, \partial_t v$, $\partial^2_t u, \partial^2_t v \in L^\infty(0,T;H^{r+1}(\Omega))$, and $U^0=\rm{Re}$$\Pi_h \varphi$, $V^0=\rm{Im}$$\Pi_h\varphi$. Then, for $0<\alpha<1$, we have \begin{equation}\label{e0010} \|\psi^n-\psi^n_h\|_0\leq C(\alpha,u, v,T^\alpha) (\tau^{2-\alpha}+h^{\tilde{r}+1}), \quad n=0,1,\ldots,N, \end{equation} and when $\alpha$ closes to $1$, the estimate becomes \begin{equation}\label{e0012} \|\psi^n-\psi^n_h\|_0\leq C(u,v,T)(\tau+h^{\tilde{r}+1}), \quad n=0,1,\ldots,N, \end{equation} where if $\beta\neq3/4$, $\tilde{r}=r$ and if $\beta=3/4$, $\tilde{r}=r-\epsilon$, $0<\epsilon<1/2$; $C(\alpha,u, v,T^\alpha)$ is only related to $\alpha$, $u$, $v$, $T^\alpha$, and $C(u,v,T)$ is only related to $u$, $v$, and $T$. \end{theorem} \begin{proof} We consider the case of $\lambda=0$. From $U^0=\operatorname{Re}\Pi_h \varphi$ and $V^0=\operatorname{Im}\Pi_h \varphi$, \eqref{e0010}, \eqref{e0012} are automatically fulfilled as $n=0$. The claims to $n\geq 1$ will be showed next. It follows from \eqref{s05}, \eqref{s06}, and \eqref{eee028} that \begin{gather*} \langle\mathcal{D}_{t}^\alpha u^n,\chi_1\rangle-\Lambda(v^n,\chi_1) =-\langle R_{u,\tau},\chi_1\rangle,\quad \forall \chi_1\in X_h,\\ \langle\mathcal{D}_{t}^\alpha v^n,\chi_2\rangle+\Lambda(u^n,\chi_2) =-\langle R_{v,\tau},\chi_2\rangle,\quad \forall \chi_2\in X_h, \end{gather*} where $R_{u,\tau}$, $R_{v,\tau}$ are the truncated errors in time. Let $\widetilde{U}^n=\Pi_hu^n$, $\widetilde{V}^n=\Pi_hv^n$, $\varepsilon_u^n=u^n-\widetilde{U}^n$, and $\varepsilon_v^n=v^n-\widetilde{V}^n$. Using \eqref{eq0068}, we obtain \begin{gather*} \langle\mathcal{D}_{t}^\alpha \widetilde{U}^n,\chi_1\rangle -\Lambda(\widetilde{V}^n,\chi_1) =-\langle \mathcal{D}_{t}^\alpha\varepsilon_u^n,\chi_1\rangle -\langle R_{u,\tau},\chi_1\rangle, \\ \langle\mathcal{D}_{t}^\alpha \widetilde{V}^n,\chi_2\rangle +\Lambda(\widetilde{U}^n,\chi_2) =-\langle \mathcal{D}_{t}^\alpha\varepsilon_v^n,\chi_2\rangle -\langle R_{v,\tau},\chi_2\rangle, \end{gather*} i.e., there exist \begin{gather} \begin{aligned} \langle\widetilde{U}^n,\chi_1\rangle-G_\alpha\Lambda(\widetilde{V}^n,\chi_1) &=(1-b_1) \langle\widetilde{U}^{n-1},\chi_1\rangle +\sum^{n-1}_{j=2}(b_{j-1}-b_{j})\langle\widetilde{U}^{n-j}, \chi_1\rangle\\ &\quad + b_{n-1}\langle\widetilde{U}^{0},\chi_1\rangle -G_\alpha\langle\gamma^n_u,\chi_1\rangle, \end{aligned} \label{s19}\\ \begin{aligned} \langle\widetilde{V}^n,\chi_2\rangle+G_\alpha\Lambda(\widetilde{U}^n,\chi_2) &=(1-b_1) \langle\widetilde{V}^{n-1},\chi_2\rangle +\sum^{n-1}_{j=2}(b_{j-1}-b_{j})\langle\widetilde{V}^{n-j},\chi_2\rangle\\ &\quad +b_{n-1}\langle\widetilde{V}^{0},\chi_2\rangle -G_\alpha\langle\gamma^n_v,\chi_2\rangle, \end{aligned} \label{s20} \end{gather} with terms $\gamma^n_u=\mathcal{D}_{t}^\alpha\varepsilon_u^n+R_{u,\tau}$, $\gamma^n_v=\mathcal{D}_{t}^\alpha \varepsilon_v^n+R_{v,\tau}$. Let $e^n_u=U^n-\widetilde{U}^n$, $e^n_v=V^n-\widetilde{V}^n$. Then, subtracting \eqref{s19}-\eqref{s20} from \eqref{eh003}-\eqref{s10} declares the residual equations \begin{gather} \begin{aligned} \langle e^n_u,\chi_1\rangle-G_\alpha\Lambda(e^n_v,\chi_1) &=(1-b_1)\langle e^{n-1}_u,\chi_1\rangle +\sum^{n-1}_{j=2}(b_{j-1}-b_{j})\langle e^{n-j}_u,\chi_1\rangle \\ &\quad + b_{n-1}\langle e^0_u,\chi_1\rangle + G_\alpha\langle \gamma^n_u,\chi_1\rangle, \end{aligned} \label{s17}\\ \begin{aligned} \langle e^n_v,\chi_2\rangle+G_\alpha\Lambda(e^n_u,\chi_2) &=(1-b_1)\langle e^{n-1}_v,\chi_2\rangle +\sum^{n-1}_{j=2}(b_{j-1}-b_{j})\langle e^{n-j}_v,\chi_2\rangle \\ &\quad + b_{n-1}\langle e^0_v,\chi_2\rangle + G_\alpha\langle\gamma^n_v,\chi_2\rangle, \end{aligned} \label{s18} \end{gather} which yield, by setting $\chi_1=e^n_u$, $\chi_2=e^n_v$, and adding \eqref{s17}, \eqref{s18}, that \begin{align*} \|e^n_u\|^2_0+\|e^n_v\|^2_0 &=\sum^{n-1}_{j=1}(b_{j-1}-b_{j})(\langle e^{n-j}_u,e^n_u \rangle+\langle e^{n-j}_v,e^n_v\rangle)\\ &\quad +b_{n-1}(\langle e^0_u,e^n_u\rangle + \langle e^0_v,e^n_v\rangle) +G_\alpha(\langle\gamma^n_u,e^n_u\rangle+\langle\gamma^n_v,e^n_v\rangle), \end{align*} where Theorem \ref{Th3} is also applied. Handled by H\"older's inequality, i.e., for all $a_1,b_1,a_2,b_2 >0$: \[ a_1b_1+a_2b_2\leq(a^2_1+a^2_2)^{1/2}(b^2_1+b^2_2)^{1/2}, \] it leads to \begin{align*} \|\boldsymbol{e}^n_h\|^2_0 &\leq\sum^{n-1}_{j=1}( b_{j-1}-b_{j})(\|e^{n-j}_u\|_0\|e^n_u\|_0 +\|e^{n-j}_v\|_0\|e^n_v\|_0) \\ &\quad+b_{n-1}(\|e^0_u\|_0\|e^n_u\|_0+\|e^0_v\|_0\|e^n_v\|_0) +G_\alpha(\|\gamma^n_u\|_0\|e^n_u\|_0+\|\gamma^n_v\|_0\|e^n_v\|_0)\\ &\leq\sum^{n-1}_{j=1}( b_{j-1}-b_{j})\|\boldsymbol{e}^{n-j}_h \|_0\|\boldsymbol{e}^{n}_h\|_0+b_{n-1}\|\boldsymbol{e}^0_h\|_0 \|\boldsymbol{e}^n_h\|_0+G_\alpha\|\boldsymbol{\gamma}^n_\tau\|_0 \|\boldsymbol{e}^n_h\|_0, \end{align*} with norms $\|\boldsymbol{e}^{n-j}_h\|^2_0=\|e^{n-j}_u\|^2_0+\|e^{n-j}_v\|^2_0$, $\|\boldsymbol{\gamma}^{n}_\tau\|^2_0=\|\gamma^{n}_u\|^2_0 +\|\gamma^{n}_v\|^2_0$, $j=0,1,\ldots,n$. This further implies the following result \begin{align*} \|\boldsymbol{e}^n_h\|_0\leq (1-b_1)\|\boldsymbol{e}^{n-1}_h\|_0 +\sum^{n-1}_{j=2}(b_{j-1}-b_{j}) \|\boldsymbol{e}^{n-j}_h\|_0+ b_{n-1}\|\boldsymbol{e}^0_h\|_0 +G_\alpha\|\boldsymbol{\gamma}^n_\tau\|_0. \end{align*} Also, it can be deduced that \begin{align*} \|\mathcal{D}_{t}^\alpha\varepsilon^n_u\|_0 &\leq \|\partial^\alpha_t\varepsilon_u^n\|_0 +C\tau^{2-\alpha}\max _{0\leq t\leq t_n}\|\partial^2_t\varepsilon_u\|_0 \\ &\leq C h^{\tilde{r}+1}\|\partial^\alpha_tu^n\|_{r+1} +C\tau^{2-\alpha}h^{\tilde{r}+1} \max _{0\leq t\leq t_n}\|\partial^2_tu\|_{r+1}, \end{align*} which gives \begin{align*} \|\gamma^n_u\|_0 &\leq C h^{\tilde{r}+1}\|\partial^\alpha_tu^n\|_{r+1} + C\tau^{2-\alpha}\max _{0\leq t\leq t_n}\|\partial^2_tu\|_{0}\\ &\quad + C\tau^{2-\alpha}h^{\tilde{r}+1} \max _{0\leq t\leq t_n}\|\partial^2_tu\|_{r+1}. \end{align*} Omitting $C\tau^{2-\alpha}h^{\tilde{r}+1} \max _{0\leq t\leq t_n}\|\partial^2_tu\|_{r+1}$, we thus obtain \begin{align}\label{e0016} \|\gamma^n_u\|_0\leq C\tau^{2-\alpha}\max _{0\leq t\leq t_n}\|\partial^2_tu\|_{0} +C h^{\tilde{r}+1}\|\partial^\alpha_tu^n\|_{r+1}, \end{align} and $\|\gamma^n_u\|_0\sim\|\gamma^n_v\|_0\sim\|\boldsymbol{\gamma}^n_\tau\|_0$. Now, on the foregoing discussion, with Lemma \ref{le068}, \eqref{e0016}, and $\|\boldsymbol{e}^0_h\|_0=0$, we sum up the results as \cite{Reh06}, i.e., for $n=1,2,\ldots,N$, and $\|\partial^\theta_t\psi\|^2_s=\|\partial^\theta_tu\|^2_s +\|\partial^\theta_tv\|^2_s$, $\theta=2$ or $\alpha$, $s=0$ or $r+1$, when $0<\alpha<1$, one has \begin{align*} \|\boldsymbol{e}^n_h\|_0 &\leq b^{-1}_{n-1}G_\alpha\max _{0\leq j\leq n}\|\boldsymbol{\gamma}^j_\tau\|_{0}\\ &\leq \frac{ n^{\alpha}\tau^\alpha\Gamma(2-\alpha)}{1-\alpha}\max _{0\leq j\leq n}\|\boldsymbol{\gamma}^j_\tau\|_{0}\\ &\leq\frac{T^\alpha\Gamma(2-\alpha)}{1-\alpha}\Big(C\tau^{2-\alpha}\max _{0\leq t\leq t_n}\|\partial^2_t\psi\|_{0} +C h^{\tilde{r}+1}\max _{0\leq j\leq n}\|\partial^\alpha_t\psi^j\|_{r+1}\Big)\\ &\leq C(\alpha,u,v,T^\alpha)(\tau^{2-\alpha}+h^{\tilde{r}+1}), \end{align*} and when $\alpha$ is very close to $1$, the estimate turns into \begin{align*} \|\boldsymbol{e}^n_h\|_0 &\leq n\tau\max _{0\leq j\leq n}\|\boldsymbol{\gamma}^j_\tau\|_{0}\\ &\leq Cn\tau\big(\tau\max _{0\leq t\leq t_n}\|\partial^2_t\psi\|_{0} +h^{\tilde{r}+1}\max _{0\leq j\leq n}\|\partial_t\psi^j\|_{r+1}\big)\\ &\leq C(u,v,T)(\tau+h^{\tilde{r}+1}). \end{align*} Hence, with Lemma \ref{le008} and the triangle inequality \begin{align*} \|\psi^n-\psi^n_h\|_0 &\leq \|\psi^n_h-\widetilde{\psi}^n\|_0+\|\psi^n-\widetilde{\psi}^n\|_0\\ &\leq\|\boldsymbol{e}^n_h\|_0+\|u^n-\widetilde{U}^n\|_0+\|v^n-\widetilde{V}^n\|_0, \end{align*} \eqref{e0010}, \eqref{e0012} can be established. The proof is complete. \end{proof} \section{Extension to two-dimensional fractional NLS}\label{NE} In this section, we extend the derived method to the two-dimensional time-space-fractional NLS, which has the following form \begin{align}\label{kkh1} \mathrm{i}\frac{\partial^{\alpha} \psi}{\partial t^{\alpha}} +\frac{\partial^{2\beta} \psi} {\partial|x|^{2\beta}}+\frac{\partial^{2\beta} \psi}{\partial|y|^{2\beta}} +\lambda f(|\psi|^2)\psi &=0, \quad (x,y;t) \in\Omega\times(0,T], \end{align} with the real factors in \eqref{eq001}-\eqref{e009} and the initial and boundary conditions \begin{gather} \psi(x,y,0)=\varphi(x,y), \quad (x,y)\in \Omega\cup\partial\Omega, \label{kk008}\\ \psi(x,y,t)=0, \quad (x,y;t) \in\partial\Omega\times(0,T], \label{kk009} \end{gather} where $\Omega=(a,b)\times(c,d)$, $\frac{\partial^{\alpha} \psi}{\partial t^{\alpha}}$ is defined as \eqref{yy02}, and Riesz derivative is defined by \begin{gather*} \frac{\partial^{2\beta} \psi(x,y,t)}{\partial|x|^{2\beta}} =\frac{-1}{2\cos(\beta \pi)} \big({_X}D^{2\beta}_L \psi(x,y,t)+{_X}D^{2\beta}_R \psi(x,y,t)\big),\\ {_X}D^{2\beta}_L \psi(x,y,t) = \frac{1}{\Gamma(2-{2\beta})} \frac{\partial^2}{\partial x^2} \int^x_a(x-\xi)^{1-{2\beta}}\psi(\xi,y,t)d\xi, \\ {_X}D^{2\beta}_R \psi(x,y,t) = \frac{1}{\Gamma(2-{2\beta})} \frac{\partial^2}{\partial x^2} \int^b_x(\xi-x)^{1-{2\beta}}\psi(\xi,y,t)d\xi. \end{gather*} $\frac{\partial^{2\beta}\psi}{\partial|y|^{2\beta}}$ is similar to $\frac{\partial^{2\beta}\psi}{\partial|x|^{2\beta}}$. Decompose $\psi(x,y,t)$ into its real and imaginary parts by $\psi(x,y,t)=u(x,y,t)+iv(x,y,t)$, whence, \eqref{kkh1} can be rewritten, i.e., \begin{gather*} \frac{\partial^\alpha u}{\partial t^\alpha} +\frac{\partial^{2\beta} v}{\partial|x|^{2\beta}} +\frac{\partial^{2\beta} v}{\partial|y|^{2\beta}}+\lambda f(u^2+v^2)v=0,\\ \frac{\partial^\alpha v}{\partial t^\alpha} -\frac{\partial^{2\beta}u}{\partial|x|^{2\beta}} -\frac{\partial^{2\beta}u}{\partial|y|^{2\beta}}-\lambda f(u^2+v^2)u=0, \end{gather*} so are the initial and boundary values. Let $\Sigma_h$ be a quasi-uniform family of subdivisions of $\Omega$ and the finite element subspace $X_h$ belong to $\mathcal{H}^\beta_0(\Omega)$, denoted as \begin{equation*} X_h=\{v_h\in \mathcal{H}^\beta_0(\Omega)\cap C^0(\bar{\Omega}): v_h|_D\in P_{r}(D), \quad \forall D\in \Sigma_h \}, \end{equation*} where $\mathcal{H}^\beta_0(\Omega)$ is the fractional Sobolev space. Then, the schemes for \eqref{kkh1}-\eqref{kk009} are constructed, which read: seek $U^n$, $V^n\in X_h$, for any $\chi_1$, $\chi_2\in X_h$, to satisfy \begin{gather} \begin{aligned} \langle U^n, \chi_1\rangle-G_\alpha\mathbf{\Lambda}(V^n,\chi_1) &= \sum^{n-1}_{j=1} (b_{j-1}-b_{j})\langle U^{n-j},\chi_1\rangle + b_{n-1}\langle U^0,\chi_1\rangle \\ &\quad -\lambda G_\alpha \langle f((U^n)^2+(V^n)^2)V^{n},\chi_1\rangle, \end{aligned} \label{s21} \\ \begin{aligned} \langle V^n, \chi_2\rangle+G_\alpha\mathbf{\Lambda}(U^n,\chi_2) &=\sum^{n-1}_{j=1} (b_{j-1}-b_{j})\langle V^{n-j},\chi_2\rangle + b_{n-1}\langle V^0,\chi_2\rangle \\ &\quad +\lambda G_\alpha\langle f((U^n)^2+(V^n)^2)U^{n},\chi_2\rangle, \end{aligned} \label{s22} \end{gather} subjected to the initial conditions \begin{equation}\label{s26} U^0=\operatorname{Re}\mathbf{\Pi}_h\varphi(x,y),\quad V^0=\operatorname{Im}\mathbf{\Pi}_h\varphi(x,y),\quad (x,y)\in\Omega\cup\partial\Omega, \end{equation} and zero boundary conditions \begin{equation}\label{s27} U^j = 0, \quad V^j=0, \quad 1\leq j\leq N, \quad \text{on }\partial\Omega, \end{equation} where $\mathbf{\Pi}_h$ is an operator like $\Pi_h$, and $\mathbf{\Lambda}(\cdot,\cdot)$ in \eqref{s21}-\eqref{s22} takes the form \begin{align*} \mathbf{\Lambda}(u,v) &=\frac{1}{2\cos(\beta\pi)}\Big\{\langle{_X}D^{\beta}_L u,{_X}D^{\beta}_R v\rangle +\langle{_X}D^{\beta}_R u,{_X}D^{\beta}_L v\rangle\Big\} \\ &\quad +\frac{1}{2\cos(\beta\pi)}\Big\{\langle{_Y}D^{\beta}_L u,{_Y} D^{\beta}_R v\rangle +\langle{_Y}D^{\beta}_R u,{_Y}D^{\beta}_L v\rangle\Big\}. \end{align*} \begin{theorem}[\cite{Ref0017}] \label{Th7} The bilinear form $\mathbf{\Lambda}(\cdot,\cdot)$ preserves the symmetric, continuous, and coercive properties, but with the energy norm \begin{align*} \|u\|_{E}=\big(\|u\|^2_{0}+|({_X}D^{\beta}_L u, {_X}D^{\beta}_R u)| +|({_Y}D^{\beta}_L u, {_Y}D^{\beta}_R u)|\big)^{1/2}. \end{align*} \end{theorem} \begin{theorem} The method \eqref{s21}-\eqref{s27} is stable in sense that, for $n=1,2,\ldots,N$, there exists $\|U^{n}\|^2_{0}+\|V^{n}\|^2_{0}\leq \|U^{0}\|^2_{0}+\|V^{0}\|^2_{0}$. \end{theorem} The above theorem follows from Theorem \ref{Th7} and the proof of Theorem \ref{Th5}. \begin{remark}\label{re003} \rm Solving \eqref{kkh1}-\eqref{kk009} by \eqref{s21}-\eqref{s27}, we confront severe challenge and computing burden; to cut the costs as much as possible, a feasible linearized strategy is needed to treat the nonlinear part, where, for $n\geq2$, we employ \begin{equation*} f(|\psi^n|^2)\psi^n=f(|\widehat{\psi}^n|^2)\widehat{\psi}^n +O(\tau^2),\quad\widehat{\psi}^n=2\psi^{n-1}-\psi^{n-2}. \end{equation*} Decompose $\psi^n$ and inserting it into \eqref{s21}-\eqref{s27}, we obtain explicit-implicit schemes \begin{gather} \langle\mathcal{D}_{t}^\alpha U^n,\chi_1\rangle -\mathbf{\Lambda}(V^n,\chi_1)+\lambda\langle f((\widehat{U}^n)^2 +(\widehat{V}^n)^2)\widehat{V}^n,\chi_1\rangle =0, \label{s24}\\ \langle\mathcal{D}_{t}^\alpha V^n,\chi_2\rangle+\mathbf{\Lambda}(U^n,\chi_2) -\lambda\langle f((\widehat{U}^n)^2 +(\widehat{V}^n)^2)\widehat{U}^n,\chi_2\rangle=0. \label{s25} \end{gather} The first step ought to be remained as itself, because available $U^1,V^1$ are required to start \eqref{s24}-\eqref{s25}. Moreover, if the solution $\psi(x,y,t)$ to \eqref{kkh1}-\eqref{kk009} is sufficiently regular, \eqref{s21}-\eqref{s27} and \eqref{s24}-\eqref{s25} maintain the optimal accuracy $O(\tau^{2-\alpha}+h^{r+1})$ in $L^2$-sense. The convergence will be confirmed in the subsequent experiments. \end{remark} \section{Numerical experiments}\label{NY} In this part, several numerical examples are performed to gauge the practical performance of \eqref{eh003}-\eqref{e007} and \eqref{s21}-\eqref{s27}, which also suffice to exhibit the accuracy of those methods. We use the algorithm as in \cite{Ref0099} with structured triangular meshes to assemble the stiffness matrix for the second method. In all the tests, $X_h$ is chosen to be piecewise linear; the nonlinear system of equations is solved by iteration with tolerant error 1.0e-012 and the convergent orders are computed as follows \begin{align*} \text{C. Order}=\begin{cases} \frac{\log\{e(\tau_1)/e(\tau_{2})\}}{\log\{\tau_1/\tau_{2}\}} &\text{in time},\\[4pt] \frac{\log\{e(h_1)/e(h_{2})\}}{\log\{h_1/h_{2}\}} &\text{in space}, \end{cases} \end{align*} where $e(\tau_1)$, $e(\tau_2)$, $e(h_1)$, $e(h_2)$ are the global $L^2$-norm errors (abbreviated as ``L. Error'') at stepsizes $\tau_1$, $\tau_2$, $h_1$, $h_2$, and $\tau_1\neq\tau_2$, $h_1\neq h_2$. To obtain more insights, we discuss the convergence separately by the real and imaginary parts. In accord with Theorem \ref{Th6} and Remark \ref{re003}, $O(\tau^{2-\alpha}+h^2)$ are anticipated in the sequel. \begin{example} \label{examp6.1} \rm Consider the one-dimensional Schr\"odinger type equation \begin{equation*} \mathrm{i}\frac{\partial^{\alpha}\psi}{\partial t^{\alpha}} +\frac{\partial^{2\beta}\psi}{\partial|x|^{2\beta}}+|\psi|^2\psi=g(x,t), \end{equation*} on the interval $(0,1)$, $0