An Efficient Discrete Model to Approximate the Solutions of a Nonlinear Double-Fractional Two-Component Gross–Pitaevskii-Type System

Abstract

In this work, we introduce and theoretically analyze a relatively simple numerical algorithm to solve a double-fractional condensate model. The mathematical system is a generalization of the famous Gross–Pitaevskii equation, which is a model consisting of two nonlinear complex-valued diffusive differential equations. The continuous model studied in this manuscript is a multidimensional system that includes Riesz-type spatial fractional derivatives. We prove here the relevant features of the numerical algorithm, and illustrative simulations will be shown to verify the quadratic order of convergence in both the space and time variables.

Dataset License: CC-BY-NC.

1. Introduction

There have been dramatic developments in the area of fractional calculus in recent decades [1], and many areas in applied and theoretical mathematics have benefited from these developments [2,3]. In particular, there have been substantial developments in the theory and application of numerical methods for fractional partial differential equations. For example, from a theoretical point of view, theoretical analyses of conservative finite-difference schemes to solve the Riesz space-fractional Gross–Pitaevskii system have been proposed in the literature [4], along with convergent three-step numerical methods to solve double-fractional condensates, explicit dissipation-preserving methods for Riesz space-fractional nonlinear wave equations in multiple dimensions [5], energy conservative difference schemes for nonlinear fractional Schrödinger equations [6], conservative difference schemes for the Riesz space-fractional sine-Gordon equation [7], high-order central difference schemes for Caputo fractional derivatives [8], among other examples.

It is important to point out that most of the methods mentioned above refer to discretizations for partial differential equations with fractional derivatives in space. In particular, those models consider fractional partial derivatives of the Riesz type. It is worth noting that the Riesz derivatives are linear combinations of the left and right Riemann–Liouville operators. other approaches strive to provide discretizations of systems of partial differential equations with Caputo-type derivatives in time. As mentioned above, there are reports on high-order central difference schemes for Caputo fractional derivatives [8], modified integral discretization schemes for two-point boundary value problems with Caputo fractional derivatives [9], predictor–corrector schemes for solving nonlinear delay differential equations of fractional order [10], numerical integrators for Caputo fractional differential equations with infinity memory effect at initial conditions [11], second-order schemes for the fast evaluation of the Caputo fractional derivatives [12], and homotopy perturbation methods for solving the Caputo-type fractional-order Volterra–Fredholm integro-differential equations [13].

Various systems mentioned above are capable of preserving some physical quantities, such as mass and energy; the development of numerical schemes that preserve these features is an important area of research. Historically, there have been many reports in this area, including systems of integer-order partial differential equations. For example, there are reports on the numerical solutions of conservative nonlinear Klein-Gordon [14] and sine-Gordon [15,16] equations, symplectic methods for the Schrödinger equation [17], fast and structure-preserving schemes for partial differential equations based on the discrete variational derivative method [18], structure-preserving numerical methods for partial differential equations [19], dissipative or conservative Galerkin methods using discrete partial derivatives for nonlinear evolution equations [20], among other reports [21]. These approaches have been extended to the fractional-case scenario, and have been helpful in designing numerical models that are able to preserve the mass and the energy of nonlinear systems [22,23,24].

Motivated by these developments, the present paper presents an efficient discrete model to approximate the solutions of a nonlinear double-fractional two-component Gross–Pitaevskii system. The system consists of two coupled complex-valued functions, whose dynamics are described by parabolic partial differential equations with nonlinear reactions. We consider here spatial derivatives of fractional order in the Riesz sense. The mathematical model is a complicated system, and the need to approximate its solutions is an interesting task. We propose here an easy-to-implement scheme to approximate the solutions and fully analyze them theoretically. We establish the properties of stability and convergence using a discrete form of Grönwall’s inequality. Some simulations will be provided to illustrate the performance of the scheme, and a numerical test of the convergence is provided. Here, it is worth noting that the main advantage of the present discretization is that the computer implementation is easy. Moreover, the scheme is quadratically convergent, and the experiments show that such is the case. Other discretizations [4,25] are more difficult to implement computationally, in view that a fixed point technique should be coded along with the computational algorithm. Moreover, those schemes provide more complicated conditions in order to guarantee the convergence. Finally, one of the limitations of the present methodology is that, to the best of our knowledge, it is not capable of preserving the energy of the system.

Let $p\in \mathbb{N}$ and $T\in {\mathbb{R}}^{+}$. Define $I_n=\{k\in \mathbb{N}:k\le n\}$, denote $\left\{0\right\}\cup I_n$ by ${\overline{I}}_n$, for each $n\in \mathbb{N}$, and let $a_i,b_i\in \mathbb{R}$ be such that $a_i<b_i$, for all $i\in I_p$. Define $\mathrm{\Omega }={\mathrm{\Pi }}_{i=1}^p(a_i,b_i)\subseteq {\mathbb{R}}^p$ and ${\mathrm{\Omega }}_T=\mathrm{\Omega }\times (0,T)$, and let $\mathrm{i}=\sqrt{-1}$. Agree that ${\psi }_1:{\overline{\mathrm{\Omega }}}_T\to \mathbb{C}$ and ${\psi }_2:{\overline{\mathrm{\Omega }}}_T\to \mathbb{C}$, and define $x=(x_1,\cdots ,x_p)\in \mathrm{\Omega }$. Functions defined on $\mathrm{\Omega }$ will be extended to all of ${\mathbb{R}}^p$ by letting them be equal to zero outside $\mathrm{\Omega }$.

Definition 1

(Podlubny [1]). Let Γ be the Gamma function, let $f:\mathbb{R}\to \mathbb{R}$ be a function, and assume that $n\in \mathbb{N}\cup \left\{0\right\}$ and $\alpha \in \mathbb{R}$ are such that $n-1<\alpha <n$. We define

(1)$\frac{d^{\alpha }f\left(x\right)}{{d\left|x\right|}^{\alpha }}=\frac{-1}{2cos\left(\frac{\pi \alpha }{2}\right)\mathrm{\Gamma }(n-\alpha )}\frac{d^n}{d{x}^n}{\int }_{-\infty }^{\infty }\frac{f\left(\xi \right)d\xi }{{|x-\xi |}^{\alpha +1-n}}.$

Definition 2.

Let $\psi :{\overline{\mathrm{\Omega }}}_T\to \mathbb{C}$, and fix $i\in I_p$. Let $\alpha >-1$ and n be as in the definition above. For each $\left(x,t\right)\in {\mathrm{\Omega }}_T$, we define

(2) $\frac{{\partial }^{\alpha }\psi (x,t)}{\partial |x_i{|}^{\alpha }}=\frac{-1}{2cos\left(\frac{\pi \alpha }{2}\right)\mathrm{\Gamma }(n-\alpha )}\frac{{\partial }^n}{\partial x_i^n}{\int }_{-\infty }^{\infty }\frac{\psi (x_1,\dots ,x_{i-1},\eta ,x_{i+1},\dots ,x_p,t)}{|x_i{-\eta |}^{\alpha -1}}d\eta .$

Agree also that

(3) ${△}^{\alpha }\psi (x,t)=\sum _{i=1}^p\frac{{\partial }^{\alpha }\psi }{\partial {\left|x_i\right|}^{\alpha }}(x,t).$

For the remainder, D, ${\beta }_{11}$, ${\beta }_{12}$, ${\beta }_{22}$ and $\lambda$ are real parameters, and $V:\overline{\mathrm{\Omega }}\to \mathbb{R}$. Let ${\alpha }_1$ and ${\alpha }_2$ satisfy $1<{\alpha }_1\le 2$ and $1<{\alpha }_2\le 2$, and assume that ${\varphi }_1:\overline{\mathrm{\Omega }}\to \mathbb{C}$ and ${\varphi }_2:\overline{\mathrm{\Omega }}\to \mathbb{C}$. Here, we study the following problem, for each $(x,t)\in {\mathrm{\Omega }}_T$ (see [26,27])

(4)$\begin{array}{c}\begin{array}{cc}\hfill \mathrm{i}\frac{\partial {\psi }_1}{\partial t}& =\lambda {\psi }_2+\left[-\frac{1}{2}{△}^{{\alpha }_1}+V\left(x\right)+D+{\beta }_{11}{\left|{\psi }_1\right|}^2+{\beta }_{12}{\left|{\psi }_2\right|}^2\right]{\psi }_1,\hfill \\ \hfill \mathrm{i}\frac{\partial {\psi }_2}{\partial t}& =\lambda {\psi }_1+\left[-\frac{1}{2}{△}^{{\alpha }_2}+V\left(x\right)+{\beta }_{12}{\left|{\psi }_1\right|}^2+{\beta }_{22}{\left|{\psi }_2\right|}^2\right]{\psi }_2,\hfill \end{array}\\ \mathrm{subjected}\mathrm{to}\left\{\begin{array}{cc}{\psi }_i(x,0)={\varphi }_i\left(x\right),\hfill & \forall i=1,2,\forall x\in \overline{\mathrm{\Omega }},\hfill \\ {\psi }_i(x,t)=0,\hfill & \forall i=1,2,\forall (x,t)\in ({\mathbb{R}}^p\setminus \mathrm{\Omega })\times (0,T).\hfill \end{array}\right.\end{array}$

Definition 3

(Ortigueira [28]). Suppose that $h,\alpha \in {\mathbb{R}}^{+}$ and assume that $f:\mathbb{R}\to \mathbb{R}$. We define the discrete operator

(5)${\mathrm{\Delta }}_h^{\left(\alpha \right)}f\left(x\right)=\sum _{k=-\infty }^{\infty }f(x-kh)g_k^{\left(\alpha \right)},\forall x\in \mathbb{R},$

where

(6) $g_k^{\left(\alpha \right)}=\frac{{(-1)}^k\mathrm{\Gamma }(\alpha +1)}{\mathrm{\Gamma }(\frac{\alpha }{2}+k+1)\mathrm{\Gamma }(\frac{\alpha }{2}-k+1)},\forall k\in \mathbb{Z}.$

It is important to note that, if f is sufficiently smooth, and $\alpha \in (0,1)\cup (1,2]$, then

(7) $\frac{{\partial }^{\alpha }f\left(x\right)}{{\partial \left|x\right|}^{\alpha }}=-\frac{{\mathrm{\Delta }}_h^{\alpha }f\left(x\right)}{h^{\alpha }}+\mathcal{O}\left(h^2\right),$

for almost all $x\in \mathbb{R}$ (see [29]).

2. Numerical Algorithm

For the remainder, we will let N and $M_i$ be natural numbers, with $i\in I_p$. Define $\tau =T/N$ and $h_i=\left(b_i-a_i\right)/M_i$. Let us introduce

(8)$\begin{array}{cc}\hfill x_{i,j_i}& =j_i{h}_i+a_i,\forall i\in I_p,\forall j_i\in {\overline{I}}_{M_i},\hfill \end{array}$

(9)$\begin{array}{cc}\hfill t_n& =n\tau ,\forall n\in {\overline{I}}_N.\hfill \end{array}$

Let $\overline{J}={\prod }_{i=1}^p{\overline{I}}_{M_i}$, where $J={\prod }_{i=1}^p{I}_{M_i-1}$. If $j=(j_1,\dots ,j_p)\in \overline{J}$, we define $x_j$ as $(x_{1,j_1},\dots ,x_{p,j_p})$. For each $(j,n)\in \overline{J}\times {\overline{I}}_N$, we use $(u_j^n,v_j^n)$ to denote a computational approximation to $(U_j^n,V_j^n)=({\psi }_1(x_j,t_n),{\psi }_2(x_j,t_n))$. Finally, let $\partial J$ represent the collection of all $j\in J$ with the property that $x_j\in \partial \mathrm{\Omega }$.

Definition 4.

If $w=u,v$ and $\alpha \in (0,1)\cup (1,2]$, we define the averages

(10) $\begin{array}{cc}\hfill {\mu }_t{w}_j^n& =\frac{w_j^{n+1}+w_j^n}{2},\hfill \end{array}$

(11) $\begin{array}{cc}\hfill {\mu }_t^{\left(1\right)}{w}_j^n& =\frac{w_j^{n+1}+w_j^{n-1}}{2},\hfill \end{array}$

and the differences

(12) $\begin{array}{cc}\hfill {\delta }_t{w}_j^n& =\frac{w_j^{n+1}-w_j^n}{\tau },\hfill \end{array}$

(13) $\begin{array}{cc}\hfill {\delta }_t^{\left(1\right)}{w}_j^n& =\frac{w_j^{n+1}-w_j^{n-1}}{2\tau },\hfill \end{array}$

(14) $\begin{array}{cc}\hfill {\delta }_{x_i}^{\left(\alpha \right)}{w}_j^n& =-h_i^{-\alpha }\sum _{k=0}^{M_i}{g}_{j_i-k}^{\left(\alpha \right)}{w}_{j_1,\dots ,j_{i-1},k,j_{i+1},\dots ,j_p}^n.\hfill \end{array}$

Moreover, we agree that

(15) $\begin{array}{cc}\hfill {△}_h^{\left(\alpha \right)}{w}_j^n& ={\delta }_{x_1}^{\left(\alpha \right)}{w}_j^n+{\delta }_{x_2}^{\left(\alpha \right)}{w}_j^n+\dots +{\delta }_{x_p}^{\left(\alpha \right)}{w}_j^n,\hfill \end{array}$

(16) $\begin{array}{cc}\hfill {▽}_h^{\left(\alpha \right)}& =({\delta }_{x_1}^{\left(\alpha \right)}{w}_j^n,{\delta }_{x_2}^{\left(\alpha \right)}{w}_j^n,\dots ,{\delta }_{x_p}^{\left(\alpha \right)}{w}_j^n).\hfill \end{array}$

With this nomenclature, we introduce the following numerical algorithm to approximate the solution of (4) on $\overline{\mathrm{\Omega }}\times [0,T]$, for each $(j,n)\in J\times {\overline{I}}_{N-1}$:

(17)$\begin{array}{c}\begin{array}{cc}\hfill \mathrm{i}{\delta }_t^{\left(1\right)}{u}_j^n& =\left[-\frac{1}{2}{△}_h^{\left({\alpha }_1\right)}+V_j+D\right]{\mu }_t^{\left(1\right)}{u}_j^n+\left[{\beta }_{11}{\left|u_j^n\right|}^2+{\beta }_{12}{\left|v_j^n\right|}^2\right]u_j^n+\lambda v_j^n,\hfill \\ \hfill \mathrm{i}{\delta }_t^{\left(1\right)}{v}_j^n& =\left[-\frac{1}{2}{△}_h^{\left({\alpha }_2\right)}+V_j\right]{\mu }_t^{\left(1\right)}{v}_j^n+\left[{\beta }_{12}{\left|u_j^n\right|}^2+{\beta }_{22}{\left|v_j^n\right|}^2\right]v_j^n+\lambda u_j^n,\hfill \end{array}\\ \mathrm{subjected}\mathrm{to}\left\{\begin{array}{cc}{u}_j^0={\mu }_t^{\left(1\right)}{u}_j^0={\varphi }_1\left(x_j\right),\hfill & \forall j\in J,\hfill \\ v_j^0={\mu }_t^{\left(1\right)}{v}_j^0={\varphi }_2\left(x_j\right),\hfill & \forall j\in J,\hfill \\ u_j^n=v_j^n=0,\hfill & \forall (j,n)\in \partial J\times {\overline{I}}_N,\hfill \end{array}\right.\end{array}$

It is clear that, if $u^n$, $u^{n-1}$, $v^n$ and $v^{n-1}$ are known for some $n\in I_{N-1}$, then the resulting iterative formulas wield expressions where the only unknowns are $u^{n+1}$ and $v^{n+1}$. On the other hand, $u^0$ and $v^0$ are prescribed by the initial data.

In order to determine the approximations when the time is equal to $t_1$, we employ the initial conditions ${\mu }_t^{\left(1\right)}{u}_j^0={\varphi }_1\left(x_j\right)$ and ${\mu }_t^{\left(1\right)}{v}_j^0={\varphi }_2\left(x_j\right)$, for all $j\in J$. The initial conditions yield then

(18)$\begin{array}{cc}\hfill {\delta }_t^{\left(1\right)}{u}_j^0& =\frac{u_j^1-{\varphi }_1\left(x_j\right)}{\tau },\forall j\in J\hfill \end{array}$

(19)$\begin{array}{cc}\hfill {\delta }_t^{\left(1\right)}{v}_j^0& =\frac{v_j^1-{\varphi }_2\left(x_j\right)}{\tau },\forall j\in J.\hfill \end{array}$

As a consequence, we obtain

(20)$\begin{array}{cc}\hfill u_j^1& ={\varphi }_1\left(x_j\right)+\mathrm{i}\tau \left[\frac{1}{2}{△}_h^{\left({\alpha }_1\right)}-V_j-D\right]{\varphi }_1\left(x_j\right)-\mathrm{i}\tau \lambda {\varphi }_2\left(x_j\right)\hfill \\ \hfill & -\mathrm{i}\tau \left[{\beta }_{11}{\left|{\varphi }_1\left(x_j\right)\right|}^2+{\beta }_{12}{\left|{\varphi }_2\left(x_j\right)\right|}^2\right]{\varphi }_1\left(x_j\right),\forall j\in J,\hfill \end{array}$

and

(21)$\begin{array}{cc}\hfill v_j^1& ={\varphi }_2\left(x_j\right)+i\tau \left[\frac{1}{2}{△}_h^{\left({\alpha }_2\right)}-V_j\right]{\varphi }_2\left(x_j\right)-i\tau \lambda {\varphi }_1\left(x_j\right)\hfill \\ \hfill & -i\tau \left[{\beta }_{12}{\left|{\varphi }_1\left(x_j\right)\right|}^2+{\beta }_{22}{\left|{\varphi }_2\left(x_j\right)\right|}^2\right]{\varphi }_2\left(x_j\right),\forall j\in J.\hfill \end{array}$

For the sake of convenience, we let $h=(h_1,\dots ,h_p)$, and denote with ${\mathcal{V}}_h$ the complex-valued functions whose domains are $\{x_j:j\in \overline{J}\}$. Moreover, set $w_j=w\left(x_j\right)$.

Theorem 1.

For all initial conditions, the system (17) has a unique solution.

Proof. 

We provide the proof for the one-dimensional case only, the higher dimensional scenario being similar. Beforehand, notice that $u^0$ and $v^0$ are defined by the initial conditions, and so are $u^1$ and $v^1$, as shown by the identities (20) and (21), respectively. Suppose that $u^n$, $u^{n-1}$, $v^n$ and $v^{n-1}$ are obtained. The first numerical model is expressed alternatively in the following form:

(22)$\frac{1}{2}\left(\frac{\mathrm{i}}{\tau }-\frac{g_0^{\left({\alpha }_1\right)}}{2h^{{\alpha }_1}}-V_j-D\right)u_j^{n+1}-\frac{1}{2h_1^{{\alpha }_1}}\sum _{\begin{array}{c}l=0\\ l\ne j\end{array}}^{M_1}{g}_{l-j}^{\left({\alpha }_1\right)}{u}_l^{n+1}=b\left(u_j^{n-1}\right)+a(u_j^n,v_j^n),$

for each $(j,n)\in J\times I_{N-1}$. Here,

(23)$\begin{array}{cc}\hfill a(u_j^n,v_j^n)& =({\beta }_{11}|u_j^n{|}^2+{\beta }_{12}|v_j^n{|}^2)u_j^n+\lambda v_j^n,\hfill \end{array}$

(24)$\begin{array}{cc}\hfill b\left(u_j^n\right)& =\frac{1}{2}\left(\frac{\mathrm{i}}{\tau }-\frac{1}{2}{\delta }_{x_1}^{\left({\alpha }_1\right)}+V_j+D\right)u_j^n.\hfill \end{array}$

As a consequence, the first equation of the system is rewritten as

(25)$-\frac{1}{2}A{u}^{n+1}=b\left(u^{n-1}\right)+a(u^n,v^n),$

where

(26)$A=\left(\begin{array}{cccc}{\displaystyle \frac{g_0^{\left({\alpha }_1\right)}}{2h_1^{{\alpha }_1}}}+V_0+D-{\displaystyle \frac{\mathrm{i}}{\tau }}& {\displaystyle \frac{g_{-1}^{\left({\alpha }_1\right)}}{2h_1^{{\alpha }_1}}}& \cdots & {\displaystyle \frac{g_{2-M_1}^{\left({\alpha }_1\right)}}{2h_1^{{\alpha }_1}}}\\ {\displaystyle \frac{g_1^{\left({\alpha }_1\right)}}{2h_1^{{\alpha }_1}}}& {\displaystyle \frac{g_0^{\left({\alpha }_1\right)}}{2h_1^{{\alpha }_1}}}+V_1+D-{\displaystyle \frac{\mathrm{i}}{\tau }}& \cdots & {\displaystyle \frac{g_{3-M_1}^{\left({\alpha }_1\right)}}{2h_1^{{\alpha }_1}}}\\ ⋮& ⋮& \ddots & ⋮\\ {\displaystyle \frac{g_{M_1-2}^{\left({\alpha }_1\right)}}{2h_1^{{\alpha }_1}}}& {\displaystyle \frac{g_{M_1-3}^{\left({\alpha }_1\right)}}{2h_1^{{\alpha }_1}}}& \cdots & {\displaystyle \frac{g_0^{\left({\alpha }_1\right)}}{2h_1^{{\alpha }_1}}}+V_{M_1}+D-{\displaystyle \frac{\mathrm{i}}{\tau }}\end{array}\right).$

For the existence of $u^{n+1}$, it suffices show that A is invertible [30]. Observe that

(27)$\left|a_{ii}\right|=\sqrt{{\left(\frac{g_0^{\left({\alpha }_1\right)}}{2h_1^{{\alpha }_1}}+V_i+D\right)}^2+{\left(\frac{1}{\tau }\right)}^2}.$

Also, the following hold:

(28)$\sum _{j\ne i}^{M_1}\left|a_{ij}\right|=\sum _{j\ne i}^{M_1}\left|\frac{g_{l-j}^{\left({\alpha }_1\right)}}{2h_1^{{\alpha }_1}}\right|=-\sum _{j\ne i}^{M_1}\frac{g_{l-j}^{{\alpha }_1}}{2h^{{\alpha }_1}}<-\sum _{\begin{array}{c}l=-\infty \\ l\ne j\end{array}}^{\infty }\frac{g_{l-j}^{{\alpha }_1}}{2h^{{\alpha }_1}}=\frac{g_0^{{\alpha }_1}}{2h^{{\alpha }_1}}\le \left|a_{ii}\right|.$

By the previous lemma, A is non-singular, which means that there exists a unique solution $u_{n+1}$ of (25). The existence of $v^{n+1}$ is proved similarly, noting that the second recursive equation of (17) is rewritten by

(29)$-\frac{1}{2}B{v}^{n+1}=d\left(u^{n-1}\right)+c(u^n,v^n),$

where

(30)$\begin{array}{cc}\hfill c(u_j^n,v_j^n)& =\lambda u_j^n+({\beta }_{12}|u_j^n{|}^2+{\beta }_{22}|v_j^n{|}^2)v_j^n,\hfill \end{array}$

(31)$\begin{array}{cc}\hfill d\left(v_j^n\right)& =\frac{1}{2}\left(\frac{\mathrm{i}}{\tau }-\frac{1}{2}{\delta }_{x_1}^{\left({\alpha }_2\right)}+V_j\right)v_j^n.\hfill \end{array}$

Moreover, B is the square complex matrix defined by

(32)$B=\left(\begin{array}{cccc}{\displaystyle \frac{g_0^{\left({\alpha }_2\right)}}{2h_1^{{\alpha }_2}}}+V_0-{\displaystyle \frac{\mathrm{i}}{\tau }}& {\displaystyle \frac{g_{-1}^{\left({\alpha }_2\right)}}{2h_1^{{\alpha }_2}}}& \cdots & {\displaystyle \frac{g_{2-M_1}^{\left({\alpha }_2\right)}}{2h_1^{{\alpha }_2}}}\\ {\displaystyle \frac{g_1^{\left({\alpha }_2\right)}}{2h_1^{{\alpha }_2}}}& {\displaystyle \frac{g_0^{\left({\alpha }_2\right)}}{2h_1^{{\alpha }_2}}}+V_1-{\displaystyle \frac{\mathrm{i}}{\tau }}& \cdots & {\displaystyle \frac{g_{3-M_1}^{\left({\alpha }_2\right)}}{2h_1^{{\alpha }_2}}}\\ ⋮& ⋮& \ddots & ⋮\\ {\displaystyle \frac{g_{M_1-2}^{\left({\alpha }_2\right)}}{2h_1^{{\alpha }_2}}}& {\displaystyle \frac{g_{M_1-3}^{\left({\alpha }_2\right)}}{2h_1^{{\alpha }_2}}}& \cdots & {\displaystyle \frac{g_0^{\left({\alpha }_2\right)}}{2h_1^{{\alpha }_2}}}+V_{M_1}-{\displaystyle \frac{\mathrm{i}}{\tau }}\end{array}\right),$

which is also strictly diagonally dominant. The conclusion follows now by induction. □

The discretization proposed in this work is similar to a linear implicit discretization of (4) in various senses. Indeed, notice that our approach hinges on approximating the nonlinear term at the time $t_n$, while the linear terms are approximated at the time $t_{n+1}$. The difference is that the linear term of the numerical model (17) is approximated by the average of the numerical solutions at the levels $n+1$ and $n-1$ through ${\mu }_t^{\left(1\right)}{u}_j^n$ and ${\mu }_t^{\left(1\right)}{v}_j^n$. In that sense, the present discretization would seem computationally more complex than the linear implicit scheme. In this point, we would like to clarify that the linear implicit scheme has the advantage of being a two-step method, but the computational implementation would require solving systems of linear equations similar to those associated to the discrete model (17). On the other hand, as we will see in the following section, our current discretization has convergence of the second order in space, while the corresponding order of the linear implicit scheme is known to be linear.

3. Computational Properties

To prove the consistency, let us define the continuous operators

(33)$\begin{array}{cc}\hfill {\mathcal{L}}_1({\psi }_1,{\psi }_2)& =\mathrm{i}\frac{\partial {\psi }_1}{\partial t}-\lambda {\psi }_2+\left[\frac{1}{2}{△}^{{\alpha }_1}-V\left(x\right)-D-{\beta }_{11}{\left|{\psi }_1\right|}^2-{\beta }_{12}{\left|{\psi }_2\right|}^2\right]{\psi }_1,\hfill \end{array}$

(34)$\begin{array}{cc}\hfill {\mathcal{L}}_2({\psi }_1,{\psi }_2)& =\mathrm{i}\frac{\partial {\psi }_2}{\partial t}-\lambda {\psi }_1+\left[\frac{1}{2}{△}^{{\alpha }_2}-V\left(x\right)-{\beta }_{12}{\left|{\psi }_1\right|}^2-{\beta }_{22}{\left|{\psi }_2\right|}^2\right]{\psi }_2,\hfill \end{array}$

for each $(x,t)\in {\mathrm{\Omega }}_T$. Define the operators

(35)$\begin{array}{cc}\hfill L_1(u_j^n,v_j^n)& =\mathrm{i}{\delta }_t^{\left(1\right)}{u}_j^n-\lambda v_j^n+\left[\frac{1}{2}{△}_h^{\left({\alpha }_1\right)}-V_j-D\right]{\mu }_t^{\left(1\right)}{u}_j^n-\left[{\beta }_{11}{\left|u_j^n\right|}^2+{\beta }_{12}{\left|v_j^n\right|}^2\right]u_j^n,\hfill \end{array}$

(36)$\begin{array}{cc}\hfill L_2(u_j^n,v_j^n)& =\mathrm{i}{\delta }_t^{\left(1\right)}{v}_j^n-\lambda u_j^n+\left[\frac{1}{2}{△}_h^{{\alpha }_2}-V_j\right]{\mu }_t^{\left(1\right)}{v}_j^n-\left[{\beta }_{12}{\left|u_j^n\right|}^2+{\beta }_{22}{\left|v_j^n\right|}^2\right]v_j^n.\hfill \end{array}$

Finally, for each $(x,t)\in {\mathrm{\Omega }}_T$ and $(j,n)\in J\times I_{N-1}$, we let

(37)$\begin{array}{cc}\hfill \mathcal{L}({\psi }_1,{\psi }_2)& =\left({\mathcal{L}}_1({\psi }_1,{\psi }_2),{\mathcal{L}}_2({\psi }_1,{\psi }_2)\right),\hfill \end{array}$

(38)$\begin{array}{cc}\hfill L({\psi }_1,{\psi }_2)& =(L_1({\psi }_1,{\psi }_2),L_2({\psi }_1,{\psi }_2)).\hfill \end{array}$

Theorem 2.

The numerical model (17) yields quadratically consistent approximations to the solutions of (4).

Proof. 

Suppose that the functions u, v, and V are sufficiently smooth. Then there exist real numbers $C_{1,k}$, $C_{2,k,i}$, and $C_{3,k}$, such that

(39)$\begin{array}{cc}\hfill \left|{\delta }_t{\psi }_k(x_j,t_n)-\frac{\partial {\psi }_k}{\partial t}(x_j,t_n)\right|& \le C_{1,k}{\tau }^2,\hfill \end{array}$

(40)$\begin{array}{cc}\hfill \left|{\mu }_t^{\left(1\right)}{\delta }_{x_i}^{\left({\alpha }_1\right)}{\psi }_k(x_j,t_n)-\frac{{\partial }^{{\alpha }_1}{\psi }_k}{{\partial \left|x\right|}^{{\alpha }_1}}(x_j,t_n)\right|& \le C_{2,k,i}({\tau }^2+h_i^2),\hfill \end{array}$

(41)$\begin{array}{cc}\hfill \left|V\left(x_j\right){\mu }_t^{\left(1\right)}{\psi }_k(x_j,t_n)-V\left(x_j\right){\psi }_k(x_j,t_n)\right|& \le C_{3,k}{\tau }^2.\hfill \end{array}$

The conclusion follows from the triangle inequality and Taylor’s theorem. □

In the following, it is worth recalling that the square-root operator of $-{△}_h^{\left(\alpha \right)}$ is the discrete fractional operator ${△}_h^{(\alpha /2)}$, for each $\alpha \in (1,2]$.

It is easy to check that, if $\alpha \in (1,2]$ and $w={\left(w^n\right)}_{n=1}^N\subseteq {\mathcal{V}}_h$, then

(42)$\begin{array}{cc}\hfill Im⟨\mathrm{i}{\delta }_t^{\left(1\right)}{w}^n,2{\mu }_t^{\left(1\right)}{w}^n⟩& ={\delta }_t^{\left(1\right)}{\parallel w^n\parallel }_2^2,\forall n\in I_{N-1},\hfill \end{array}$

(43)$\begin{array}{cc}\hfill Im⟨-{△}_h^{\left(\alpha \right)}{\mu }_t^{\left(1\right)}{w}^n,2{\mu }_t^{\left(1\right)}{w}^n⟩& =⟨-{△}_h^{\left(\alpha \right)}{\mu }_t^{\left(1\right)}{w}^n,2{\mu }_t^{\left(1\right)}{w}^n⟩=0,\forall n\in I_{N-1}.\hfill \end{array}$

These identities will be employed in the proofs of Theorems 3 and 4. Similarly, the following result will be crucial in those proofs.

Lemma 1.

If $ϵ={\left({ϵ}^n\right)}_{n=0}^N$ and $\zeta ={\left({\zeta }^n\right)}_{n=0}^N$ belong to ${\mathcal{V}}_h$, then

(44) $\begin{array}{cc}\hfill \sum _{n=1}^m\left[Im⟨{ϵ}^n,{\mu }_t^{\left(1\right)}{\zeta }^n⟩+Im⟨{\zeta }^n,{\mu }_t^{\left(1\right)}{ϵ}^n⟩\right]& =Im⟨{ϵ}^m,{\zeta }^{m+1}⟩+Im⟨{\zeta }^m,{ϵ}^{m+1}⟩\hfill \\ \hfill & -Im⟨{ϵ}^0,{\zeta }^1⟩-Im⟨{\zeta }^0,{ϵ}^1⟩.\hfill \end{array}$

Proof. 

It is easy to check that

(45)$\begin{array}{cc}\hfill & \sum _{n=1}^m\left[Im⟨{ϵ}^n,{\mu }_t^{\left(1\right)}{\zeta }^n⟩+Im⟨{\zeta }^n,{\mu }_t^{\left(1\right)}{ϵ}^n⟩\right]\hfill \\ \hfill & =\sum _{n=1}^{m}Im\left[⟨{ϵ}^n,{\zeta }^{n+1}⟩+\overline{⟨{ϵ}^{n-1},{\zeta }^n⟩}+⟨{\zeta }^n,{ϵ}^{n+1}⟩+\overline{⟨{\zeta }^{n-1},{ϵ}^n⟩}\right]\hfill \\ \hfill & =Im⟨{ϵ}^m,{\zeta }^{m+1}⟩+\sum _{n=1}^{m-1}Im\left[⟨{ϵ}^n,{\zeta }^{n+1}⟩+\overline{⟨{ϵ}^n,{\zeta }^{n+1}⟩}\right]+\overline{⟨{ϵ}^0,{\zeta }^1⟩}\hfill \\ \hfill & +Im⟨{\zeta }^m,{ϵ}^{m+1}⟩+\sum _{n=1}^{m-1}Im\left[⟨{\zeta }^n,{ϵ}^{n+1}⟩+\overline{⟨{\zeta }^n,{ϵ}^{n+1}⟩}\right]+\overline{⟨{\zeta }^0,{ϵ}^1⟩},\hfill \end{array}$

which is what we wanted to prove. □

Next, we consider initial data of the form $({\varphi }_1,{\varphi }_2)$ and $({\tilde{\varphi }}_1,{\tilde{\varphi }}_2)$. Here, ${\tilde{\varphi }}_1$ and ${\tilde{\varphi }}_2$ are both complex functions, and the numerical approximations associated to each of these pars is represented as $(u,v)$ and $(\tilde{u},\tilde{v})$, respectively.

Lemma 2.

For each $n\in {\overline{I}}_N$, let ${ϵ}^n=u^n-{\tilde{u}}^n$ and ${\zeta }^n=v^n-{\tilde{v}}^n$. Agree that

(46) $\begin{array}{cc}\hfill P_j^n& =-\left({\beta }_{11}{\left|{\tilde{u}}_j^n\right|}^2+{\beta }_{12}{\left|{\tilde{v}}_j^n\right|}^2\right){\tilde{u}}_j^n+\left({\beta }_{11}{\left|u_j^n\right|}^2+{\beta }_{12}{\left|v_j^n\right|}^2\right)u_j^n,\hfill \end{array}$

(47) $\begin{array}{cc}\hfill Q_j^n& =-\left({\beta }_{22}{\left|{\tilde{v}}_j^n\right|}^2+{\beta }_{12}{\left|{\tilde{u}}_j^n\right|}^2\right){\tilde{v}}_j^n+\left({\beta }_{22}{\left|v_j^n\right|}^2+{\beta }_{12}{\left|u_j^n\right|}^2\right)v_j^n,\hfill \end{array}$

for each $(j,n)\in J\times I_{N-1}$. Then there exists a constant $C\ge 0$, which dependents on τ and h, such that $max\left\{\right|P_j^n|,|Q_j^n\left|\right\}\le C\left(\right|{ϵ}_j^n|+|{\zeta }_j^n\left|\right)$, for each $(j,n)\in J\times I_{N-1}$.□

Theorem 3.

Let ${ϵ}^n$ and ${\zeta }^n$ be as before. For some $C^{″}\ge 0$, whenever $2C^{″}\tau <1$ holds, then

(48) ${\mu }_t(\parallel {ϵ}^n{\parallel }_2^2+\parallel {\zeta }^n{\parallel }_2^2)\le 2{\mu }_t(\parallel {ϵ}^0{\parallel }_2^2+\parallel {\zeta }^0{\parallel }_2^2)e^{2C^{″}T},\forall n\in {\overline{I}}_N.$

Proof. 

Beforehand, using the notation above, we observe that

(49)$\begin{array}{c}\begin{array}{cc}\hfill \mathrm{i}{\delta }_t^{\left(1\right)}{ϵ}_j^n& =\left[-\frac{1}{2}{△}_h^{{\alpha }_1}+V_j+D\right]{\mu }_t^{\left(1\right)}{ϵ}_j^n+P_j^n+\lambda {\zeta }_j^n,\hfill \\ \hfill \mathrm{i}{\delta }_t^{\left(1\right)}{\zeta }_j^n& =\left[-\frac{1}{2}{△}_h^{{\alpha }_2}+V_j\right]{\mu }_t^{\left(1\right)}{\zeta }_j^n+Q_j^n+\lambda {ϵ}_j^n,\hfill \end{array}\\ \mathrm{such}\mathrm{that}\left\{\begin{array}{cc}{ϵ}_j^0={\mu }_t^{\left(1\right)}{ϵ}_j^0={\varphi }_1\left(x_j\right)-{\tilde{\varphi }}_1\left(x_j\right),\hfill & \forall j\in \overline{J},\hfill \\ {\zeta }_j^0={\mu }_t^{\left(1\right)}{\zeta }_j^0={\varphi }_2\left(x_j\right)-{\tilde{\varphi }}_2\left(x_j\right),\hfill & \forall j\in \overline{J},\hfill \\ {ϵ}_j^n={\zeta }_j^n=0,\hfill & \forall (j,n)\in \partial J\times {\overline{I}}_N.\hfill \end{array}\right.\end{array}$

Thus, there is some $C^{\prime }\in {\mathbb{R}}^{+}$, such that

(50)$max\{\parallel P^n{\parallel }_2^2,\parallel Q^n{\parallel }_2^2\}\le C^{\prime }\left(\parallel {ϵ}^n{\parallel }_2^2+{\parallel {\zeta }^n\parallel }_2^2\right),\forall n\in I_{N-1}.$

Algebraic calculations establish that

(51)$\begin{array}{cc}\hfill {\delta }_t{\mu }_t{\parallel {ϵ}^{n-1}\parallel }_2^2& =2\lambda Im⟨{\zeta }^n,{\mu }_t^{\left(1\right)}{ϵ}^n⟩+2Im⟨P^n,{\mu }_t^{\left(1\right)}{ϵ}^n⟩,\forall n\in I_{N-1},\hfill \end{array}$

(52)$\begin{array}{cc}\hfill {\delta }_t{\mu }_t{\parallel {\zeta }^{n-1}\parallel }_2^2& =2\lambda Im⟨{ϵ}^n,{\mu }_t^{\left(1\right)}{\zeta }^n⟩+2Im⟨Q^n,{\mu }_t^{\left(1\right)}{\zeta }^n⟩,\forall n\in I_{N-1}.\hfill \end{array}$

These identities and bounding, showing that there exists $C^{″}\in {\mathbb{R}}^{+}$, with the property that

(53)$\begin{array}{cc}\hfill & {\mu }_t(\parallel {ϵ}^m{\parallel }_2^2+\parallel {\zeta }^m{\parallel }_2^2)\hfill \\ \hfill & ={\mu }_t(\parallel {ϵ}^0{\parallel }_2^2+\parallel {\zeta }^0{\parallel }_2^2)+2\tau \sum _{n=1}^{m}Im(⟨P^n,{\mu }_t^{\left(1\right)}{ϵ}^n⟩+⟨Q^n,{\mu }_t^{\left(1\right)}{\zeta }^n⟩)\hfill \\ \hfill & +2\lambda \tau \left(Im⟨{ϵ}^m,{\zeta }^{m+1}⟩+Im⟨{\zeta }^m,{ϵ}^{m+1}⟩-Im⟨{ϵ}^0,{\zeta }^1⟩-Im⟨{\zeta }^0,{ϵ}^1⟩\right)\hfill \\ \hfill & \le {\mu }_t(\parallel {ϵ}^0{\parallel }_2^2+\parallel {\zeta }^0{\parallel }_2^2)+2\tau \sum _{n=0}^m\left[\parallel P^n{\parallel }_2^2+{\parallel Q^n\parallel }_2^2+{\mu }_t\left(\parallel {ϵ}^n{\parallel }_2^2+{\parallel {\zeta }^n\parallel }_2^2\right)\right]\hfill \\ \hfill & +2\lambda \tau {\mu }_t\left(\parallel {ϵ}^m{\parallel }_2^2+\parallel {\zeta }^m{\parallel }_2^2+\parallel {ϵ}^0{\parallel }_2^2+{\parallel {\zeta }^0\parallel }_2^2\right)\hfill \\ \hfill & \le {\mu }_t(\parallel {ϵ}^0{\parallel }_2^2+\parallel {\zeta }^0{\parallel }_2^2)+C^{\prime \prime }\tau \sum _{n=0}^m{\mu }_t\left(\parallel {ϵ}^n{\parallel }_2^2+{\parallel {\zeta }^n\parallel }_2^2\right),\forall m\in I_{M-1}.\hfill \end{array}$

Subtract then the term $C^{″}\tau {\mu }_t(\parallel {ϵ}^m{\parallel }_2^2+\parallel {\zeta }^m{\parallel }_2^2)$ from both ends of this chain of inequalities, regroup and observe that $1-C^{″}\tau >\frac{1}{2}$. Therefore,

(54)${\mu }_t(\parallel {ϵ}^m{\parallel }_2^2+\parallel {\zeta }^m{\parallel }_2^2)\le 2{\mu }_t(\parallel {ϵ}^0{\parallel }_2^2+\parallel {\zeta }^0{\parallel }_2^2)+2C^{″}\tau \sum _{n=0}^{m-1}{\mu }_t\left(\parallel {ϵ}^n{\parallel }_2^2+{\parallel {\zeta }^n\parallel }_2^2\right).$

The theorem follows employing a suitable discrete Grönwall inequality (see [31]) with ${\rho }^n=2{\mu }_t(\parallel {ϵ}^0{\parallel }_2^2+\parallel {\zeta }^0{\parallel }_2^2)$ and ${\omega }^n={\mu }_t(\parallel {ϵ}^n{\parallel }_2^2+\parallel {\zeta }^n{\parallel }_2^2)$. □

Theorem 4.

Let $(U,V)\in {\mathcal{C}}_{x,t}\left(\overline{{\mathrm{\Omega }}_T}\right)\times {\mathcal{C}}_{x,t}\left(\overline{{\mathrm{\Omega }}_T}\right)$ be a solution of (4). Then the solution of (17) converges to $(U,V)$ with order $\mathcal{O}({\tau }^2+\parallel h{\parallel }_2^2)$ in the $L_2$-norm if τ is sufficiently small.

Proof. 

The proof is similar to that of the stability property of (17). □

Proof. 

The proof of this result is similar to that of Theorem 3 and, for that reason, we provide here a shortened proof only. Beforehand, let ${ϵ}^n=u^n-U^n$ and ${\zeta }^n=v^n-V^n$, for each $n\in {\overline{I}}_N$. It is obvious then that the pair $(ϵ,\zeta )$ satisfies the discrete problem

(55)$\begin{array}{c}\begin{array}{cc}\hfill \mathrm{i}{\delta }_t^{\left(1\right)}{ϵ}_j^n& =\left[-\frac{1}{2}{△}_h^{{\alpha }_1}+V_j+D\right]{\mu }_t^{\left(1\right)}{ϵ}_j^n+P_j^n+\lambda {\zeta }_j^n-{\rho }_j^n,\hfill \\ \hfill \mathrm{i}{\delta }_t^{\left(1\right)}{\zeta }_j^n& =\left[-\frac{1}{2}{△}_h^{{\alpha }_2}+V_j\right]{\mu }_t^{\left(1\right)}{\zeta }_j^n+Q_j^n+\lambda {ϵ}_j^n-{\sigma }_j^n,\hfill \end{array}\\ \mathrm{such}\mathrm{that}\left\{\begin{array}{cc}{ϵ}_j^0={\mu }_t^{\left(1\right)}{ϵ}_j^0=0,\hfill & \forall j\in \overline{J},\hfill \\ {\zeta }_j^0={\mu }_t^{\left(1\right)}{\zeta }_j^0=0,\hfill & \forall j\in \overline{J},\hfill \\ {ϵ}_j^n={\zeta }_j^n=0,\hfill & \forall (j,n)\in \partial J\times {\overline{I}}_N,\hfill \end{array}\right.\end{array}$

where

(56)$\begin{array}{cc}\hfill P_j^n& =\left({\beta }_{11}{\left|u_j^n\right|}^2+{\beta }_{12}{\left|v_j^n\right|}^2\right)u_j^n-\left({\beta }_{11}{\left|U_j^n\right|}^2+{\beta }_{12}{\left|V_j^n\right|}^2\right)U_j^n,\hfill \end{array}$

(57)$\begin{array}{cc}\hfill Q_j^n& =\left({\beta }_{22}{\left|v_j^n\right|}^2+{\beta }_{12}{\left|u_j^n\right|}^2\right)v_j^n-\left({\beta }_{22}{\left|V_j^n\right|}^2+{\beta }_{12}{\left|U_j^n\right|}^2\right)V_j^n,\hfill \end{array}$

for each $(j,n)\in J\times I_{N-1}$. Proceeding now as in the proof of Theorem 3, there exists a constant $C^{\prime }\ge 0$, such that (50) is satisfied. Following the same steps used to obtain (51) and (52), we may readily obtain the following identities, valid for each $n\in I_{N-1}$:

(58)$\begin{array}{cc}\hfill {\delta }_t{\mu }_t{\parallel {ϵ}^{n-1}\parallel }_2^2& =2Im⟨P^n,{\mu }_t^{\left(1\right)}{ϵ}^n⟩+2\lambda Im⟨{\zeta }^n,{\mu }_t^{\left(1\right)}{ϵ}^n⟩-2Im⟨{\rho }^n,{\mu }_t^{\left(1\right)}{ϵ}^n⟩,\hfill \end{array}$

(59)$\begin{array}{cc}\hfill {\delta }_t{\mu }_t{\parallel {\zeta }^{n-1}\parallel }_2^2& =2Im⟨Q^n,{\mu }_t^{\left(1\right)}{\zeta }^n⟩+2\lambda Im⟨{ϵ}^n,{\mu }_t^{\left(1\right)}{\zeta }^n⟩-2Im⟨{\sigma }^n,{\mu }_t^{\left(1\right)}{\zeta }^n⟩.\hfill \end{array}$

Again, as in the proof of Theorem 3, we add these last two equations and take the sum on both sides over all $n\in {\overline{I}}_m$, for some $m\in I_{N-1}$. Using telescoping sums, rearranging terms, taking then absolute values, using Young’s inequality and (50), and applying the inequality of the hypotheses, it follows that there is a constant $C^{″}\ge 0$, such that

(60)${\omega }^m\le {\rho }^m+C^{″}\tau \sum _{n=0}^{m-1}{\omega }^n,\forall m\in I_{N-1},$

for sufficiently small values of $\tau$. In this expression, ${\omega }^n$ is the constant used in the proof of Theorem 3, for each $n\in {\overline{I}}_{N-1}$. Moreover, in this case,

(61)${\rho }^m=2{\mu }_t(\parallel {ϵ}^0{\parallel }_2^2+\parallel {\zeta }^0{\parallel }_2^2)+4\tau \sum _{n=1}^m(\parallel {\rho }^n{\parallel }_2^2+\parallel {\sigma }^n{\parallel }_2^2),\forall m\in I_{N-1}.$

Letting now $C_1=8e^{C^{″}T}$, using the discrete Grönwall inequality [31] for sufficiently small values of $\tau$ and substituting the initial conditions of (55), we readily reach that

(62)$\begin{array}{cc}\hfill \parallel {ϵ}^m{\parallel }_2^2,{\parallel {\zeta }^m\parallel }_2^2& \le 2{\mu }_t(\parallel {ϵ}^m{\parallel }_2^2+\parallel {\zeta }^m{\parallel }_2^2)\le C_1\tau \sum _{n=1}^m(\parallel {\rho }^n{\parallel }_2^2+\parallel {\sigma }^n{\parallel }_2^2)\hfill \\ \hfill & \le C_0^2{C}_{1}T({\tau }^2+{\parallel h\parallel }_2^2{)}^2,\forall m\in I_{N-1}.\hfill \end{array}$

As a consequence, there exists a constant $C\ge 0$ with the property that $\parallel {ϵ}^m{\parallel }_2,\parallel {\zeta }^m{\parallel }_2\le C({\tau }^2+{\parallel h\parallel }_2^2)$, for each $m\in {\overline{I}}_N$. We conclude that the scheme (17) converges to the exact solution of (4), with the quadratic order of convergence. □

For the remainder of this work, use ${\varphi }_1\left(x\right)={\varphi }_2\left(x\right)=\frac{1}{\sqrt{\pi }}e-x^2$, $V\left(x\right)=\frac{1}{2}{x}^2$, $p=1$, $\mathrm{\Omega }=(-7,7)$, ${\beta }_{11}=1.5$, ${\beta }_{12}=0.5$, ${\beta }_{22}=1.5$, $\lambda =0.5$ and $D=2$. Set $h_1=0.1$, $\tau =0.01$ and $T=10$. For illustration purposes, Figure 1, Figure 2 and Figure 3 depicts the approximations to the solutions of (4) for various values of ${\alpha }_1$ and ${\alpha }_2$. Finally, we study numerically the convergence of the scheme (17) using $\alpha ={\alpha }_1=1.5$ and $T=0.5$. Table 1a,b summarize the computational results on the temporal and spatial convergence, and they show that (17) is quadratically convergent.

4. Conclusions

Before concluding this paper, we must note that our discretization of the double fractional Gross–Pitaevskii system follows a standard approach [4,25]. In particular, this means that the numerical discretization uses local approximations, estimating each of the terms of the mathematical model. A natural question would be whether a non-standard approach could be used to discretize our mathematical model. Here, we understand the adjective ’non-standard’ in the sense of R.E. Mickens [32]. Such topic of investigation is outside the scope of this work, but it is an interesting avenue of research for a future study. On the other hand, the constant $C\ge 0$ in Lemma 2 are not necessarily independent of h and $\tau$. This is due to the fact that the nonlinear terms in (4) are not globally Lipschitz continuous. This is one of the limitations of our study, and it remains an open problem of research. To solve this shortcoming, it would be desirable to establish the uniform boundedness of the numerical approximations obtained through (4). Following recent reports available in the literature [24], the determination of conserved positive quantities for the finite-difference method (17) would be helpful in bounding the numerical solutions [7]. The authors have attempted to employ the discrete energy method to derive such quantities [33,34]. Unfortunately, their efforts have not yielded conserved positive quantities to this day.

Author Contributions

Conceptualization, J.E.M.-D.; methodology, J.E.M.-D.; software, J.E.M.-D., N.R. and A.J.S.-R.; validation, J.E.M.-D., N.R. and A.J.S.-R.; formal analysis, J.E.M.-D. and A.J.S.-R.; investigation, J.E.M.-D.; resources, J.E.M.-D. and N.R.; data curation, J.E.M.-D. and N.R.; writing—original draft preparation, J.E.M.-D. and A.J.S.-R.; writing—review and editing, J.E.M.-D. and N.R.; visualization, J.E.M.-D.; supervision, J.E.M.-D. and N.R.; project administration, J.E.M.-D. and N.R.; funding acquisition, J.E.M.-D. and N.R. All authors have read and agreed to the published version of the manuscript.

Funding

J.E.M.-D. thanks the National Council of Science and Technology of Mexico for financially supporting him via grant A1-S-45928. Ministerio de Ciencia e Innovación and Regional Development European Funds through project PGC2018-101443-B-I00.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author, J.E.M.-D., upon reasonable request.

Acknowledgments

The authors would like to thank the anonymous reviewers and the associate editor for their time and comments. All of their suggestions were strictly followed, and the result was a substantial improvement in the final quality of this work.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Table

Figure 1. Graphs of the approximate solutions of problem (4) vs. [Forumla omitted. See PDF.], employing [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. The graphs show (a) [Forumla omitted. See PDF.], (b) [Forumla omitted. See PDF.], (c) [Forumla omitted. See PDF.], (d) [Forumla omitted. See PDF.], (e) [Forumla omitted. See PDF.] and (f) [Forumla omitted. See PDF.].

Figure 2. Graphs of the approximate solutions of problem (4) vs. [Forumla omitted. See PDF.], employing [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. The graphs show (a) [Forumla omitted. See PDF.], (b) [Forumla omitted. See PDF.], (c) [Forumla omitted. See PDF.], (d) [Forumla omitted. See PDF.], (e) [Forumla omitted. See PDF.] and (f) [Forumla omitted. See PDF.].

Figure 3. Graphs of the approximate solutions of problem (4) vs. [Forumla omitted. See PDF.], employing [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. The graphs show (a) [Forumla omitted. See PDF.], (b) [Forumla omitted. See PDF.], (c) [Forumla omitted. See PDF.], (d) [Forumla omitted. See PDF.], (e) [Forumla omitted. See PDF.] and (f) [Forumla omitted. See PDF.].

Figure 3. Graphs of the approximate solutions of problem (4) vs. [Forumla omitted. See PDF.], employing [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. The graphs show (a) [Forumla omitted. See PDF.], (b) [Forumla omitted. See PDF.], (c) [Forumla omitted. See PDF.], (d) [Forumla omitted. See PDF.], (e) [Forumla omitted. See PDF.] and (f) [Forumla omitted. See PDF.].

Table 1

Computational study of convergence.

(a) Time
	$$\mathit{h}$$
	$$\mathbf{0.04}$$		$$\mathbf{0.02}$$		$$\mathit{h}=\mathbf{0.01}$$
$$\mathit{\tau }$$	$$\parallel {\mathit{ϵ}}_{\mathit{t},\mathit{h}}{\parallel }_{\mathbf{2}}$$	$${\mathit{\rho }}_{\mathit{\tau },\mathit{h}}$$	$$\parallel {\mathit{ϵ}}_{\mathit{t},\mathit{h}}{\parallel }_{\mathbf{2}}$$	$${\mathit{\rho }}_{\mathit{\tau },\mathit{h}}$$	$$\parallel {\mathit{ϵ}}_{\mathit{t},\mathit{h}}{\parallel }_{\mathbf{2}}$$	$${\mathit{\rho }}_{\mathit{\tau },\mathit{h}}$$
$$0.02/2^0$$	$$1.2783\times {10}^{-2}$$	−	$$3.4758\times {10}^{-3}$$	−	$$8.8854\times {10}^{-4}$$	−
$$0.02/2^1$$	$$3.5439\times {10}^{-3}$$	$$1.8508$$	$$9.2893\times {10}^{-4}$$	$$1.9037$$	$$2.2993\times {10}^{-4}$$	$$1.9502$$
$$0.02/2^2$$	$$9.5353\times {10}^{-4}$$	$$1.8940$$	$$2.5436\times {10}^{-4}$$	$$1.8687$$	$$5.8602\times {10}^{-5}$$	$$1.9722$$
$$0.02/2^3$$	$$2.4988\times {10}^{-4}$$	$$1.9320$$	$$6.6350\times {10}^{-5}$$	$$1.9387$$	$$1.4386\times {10}^{-5}$$	$$2.0263$$
$$0.02/2^4$$	$$6.7002\times {10}^{-5}$$	$$1.8990$$	$$1.7215\times {10}^{-5}$$	$$1.9464$$	$$3.7194\times {10}^{-6}$$	$$1.9515$$
(b) Space
	$$\mathit{\tau }$$
	$$\mathbf{0.02}$$		$$\mathbf{0.01}$$		$$\mathbf{0.005}$$
$$\mathit{h}$$	$$\parallel {\mathit{ϵ}}_{\mathit{t},\mathit{h}}{\parallel }_{\mathbf{2}}$$	$${\mathit{\sigma }}_{\mathit{\tau },\mathit{h}}$$	$$\parallel {\mathit{ϵ}}_{\mathit{t},\mathit{h}}{\parallel }_{\mathbf{2}}$$	$${\mathit{\sigma }}_{\mathit{\tau },\mathit{h}}$$	$$\parallel {\mathit{ϵ}}_{\mathit{t},\mathit{h}}{\parallel }_{\mathbf{2}}$$	$${\mathit{\sigma }}_{\mathit{\tau },\mathit{h}}$$
$$0.08$$	$$4.6852\times {10}^{-2}$$	−	$$1.3239\times {10}^{-2}$$	−	$$3.5547\times {10}^{-3}$$	−
$$0.04$$	$$1.2783\times {10}^{-2}$$	$$1.8739$$	$$3.5439\times {10}^{-3}$$	$$1.9003$$	$$9.5353\times {10}^{-4}$$	$$1.8984$$
$$0.02$$	$$3.4758\times {10}^{-3}$$	$$1.8508$$	$$9.2893\times {10}^{-4}$$	$$1.9316$$	$$2.5436\times {10}^{-4}$$	$$1.9064$$
$$0.01$$	$$8.8854\times {10}^{-4}$$	$$1.9958$$	$$2.2993\times {10}^{-4}$$	$$2.0143$$	$$5.8602\times {10}^{-5}$$	$$2.1178$$
$$0.005$$	$$2.4345\times {10}^{-4}$$	$$1.8678$$	$$6.6720\times {10}^{-5}$$	$$1.7850$$	$$1.6445\times {10}^{-5}$$	$$1.8333$$

References

1. Podlubny, I. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications; Elsevier: Amsterdam, The Netherlands, 1998; Volume 198.

2. Saeed, A.; Bilal, M.; Gul, T.; Kumam, P.; Khan, A.; Sohail, M.M. Fractional Order Stagnation Point Flow of the Hybrid Nanofluid towards a Stretching Sheet. Sci. Rep.; 2021; 11, 20429. [DOI: https://dx.doi.org/10.1038/s41598-021-00004-3] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34650086]

3. Gul, T.; Alghamdi, W.; Khan, I.; Ali, I. New similarity variable to transform the fluid flow from PDEs into fractional-order ODEs: Numerical study. Phys. Scr.; 2021; 96, 084009. [DOI: https://dx.doi.org/10.1088/1402-4896/abf377]

4. Serna-Reyes, A.J.; Macías-Díaz, J. Theoretical analysis of a conservative finite-difference scheme to solve a Riesz space-fractional Gross–Pitaevskii system. J. Comput. Appl. Math.; 2021; 113413. [DOI: https://dx.doi.org/10.1016/j.cam.2021.113413]

5. Macías-Díaz, J. An explicit dissipation-preserving method for Riesz space-fractional nonlinear wave equations in multiple dimensions. Commun. Nonlinear Sci. Numer. Simul.; 2018; 59, pp. 67-87. [DOI: https://dx.doi.org/10.1016/j.cnsns.2017.10.019]

6. Wang, P.; Huang, C. An energy conservative difference scheme for the nonlinear fractional Schrödinger equations. J. Comput. Phys.; 2015; 293, pp. 238-251. [DOI: https://dx.doi.org/10.1016/j.jcp.2014.03.037]

7. Xing, Z.; Wen, L. A conservative difference scheme for the Riesz space-fractional sine-Gordon equation. Adv. Differ. Equ.; 2018; 2018, 238. [DOI: https://dx.doi.org/10.1186/s13662-018-1689-5]

8. Ying, Y.; Lian, Y.; Tang, S.; Liu, W.K. High-order central difference scheme for Caputo fractional derivative. Comput. Methods Appl. Mech. Eng.; 2017; 317, pp. 42-54. [DOI: https://dx.doi.org/10.1016/j.cma.2016.12.008]

9. Cen, Z.; Huang, J.; Xu, A.; Le, A. A modified integral discretization scheme for a two-point boundary value problem with a Caputo fractional derivative. J. Comput. Appl. Math.; 2020; 367, 112465. [DOI: https://dx.doi.org/10.1016/j.cam.2019.112465]

10. Bhalekar, S.; Daftardar-Gejji, V. A predictor-corrector scheme for solving nonlinear delay differential equations of fractional order. J. Fract. Calc. Appl.; 2011; 1, pp. 1-9.

11. Mendes, E.M.; Salgado, G.H.; Aguirre, L.A. Numerical solution of Caputo fractional differential equations with infinity memory effect at initial condition. Commun. Nonlinear Sci. Numer. Simul.; 2019; 69, pp. 237-247. [DOI: https://dx.doi.org/10.1016/j.cnsns.2018.09.022]

12. Yan, Y.; Sun, Z.Z.; Zhang, J. Fast evaluation of the Caputo fractional derivative and its applications to fractional diffusion equations: A second-order scheme. Commun. Comput. Phys.; 2017; 22, pp. 1028-1048. [DOI: https://dx.doi.org/10.4208/cicp.OA-2017-0019]

13. Das, P.; Rana, S.; Ramos, H. Homotopy perturbation method for solving Caputo-type fractional-order Volterra–Fredholm integro-differential equations. Comput. Math. Methods; 2019; 1, e1047. [DOI: https://dx.doi.org/10.1002/cmm4.1047]

14. Strauss, W.; Vazquez, L. Numerical solution of a nonlinear Klein-Gordon equation. J. Comput. Phys.; 1978; 28, pp. 271-278. [DOI: https://dx.doi.org/10.1016/0021-9991(78)90038-4]

15. Ben-Yu, G.; Pascual, P.J.; Rodriguez, M.J.; Vázquez, L. Numerical solution of the sine-Gordon equation. Appl. Math. Comput.; 1986; 18, pp. 1-14. [DOI: https://dx.doi.org/10.1016/0096-3003(86)90025-1]

16. Fei, Z.; Vázquez, L. Two energy conserving numerical schemes for the sine-Gordon equation. Appl. Math. Comput.; 1991; 45, pp. 17-30. [DOI: https://dx.doi.org/10.1016/0096-3003(91)90087-4]

17. Tang, Y.F.; Vázquez, L.; Zhang, F.; Pérez-García, V. Symplectic methods for the nonlinear Schrödinger equation. Comput. Math. Appl.; 1996; 32, pp. 73-83. [DOI: https://dx.doi.org/10.1016/0898-1221(96)00136-8]

18. Furihata, D. Fast and structure-preserving schemes for PDEs based on discrete variational derivative method. Geom. Numer. Integr. Its Appl.; 2017; 74, pp. 19-22.

19. Furihata, D.; Matsuo, T. Discrete Variational Derivative Method: A Structure-Preserving Numerical Method for Partial Differential Equations; CRC Press: New York, NY, USA, 2010.

20. Matsuo, T. Dissipative/conservative Galerkin method using discrete partial derivatives for nonlinear evolution equations. J. Comput. Appl. Math.; 2008; 218, pp. 506-521. [DOI: https://dx.doi.org/10.1016/j.cam.2007.08.001]

21. Li, X.G.; Li, S.C. High-order numerical methods with mass and energy conservation for spin–orbit-coupled Bose–Einstein condensates. Int. J. Comput. Math.; 2021; 98, pp. 738-757. [DOI: https://dx.doi.org/10.1080/00207160.2020.1777407]

22. Zhai, L.; Wang, J. High-order conservative scheme for the coupled space fractional nonlinear Schrödinger equations. Int. J. Comput. Math.; 2021; pp. 1-22. [DOI: https://dx.doi.org/10.1080/00207160.2021.1925889]

23. Fu, Y.; Shi, Y.; Zhao, Y. Explicit high-order structure-preserving algorithms for the two-dimensional fractional nonlinear Schrödinger equation. Int. J. Comput. Math.; 2021; pp. 1-28. [DOI: https://dx.doi.org/10.1080/00207160.2021.1940978]

24. Macías-Díaz, J.E. On the solution of a Riesz space-fractional nonlinear wave equation through an efficient and energy-invariant scheme. Int. J. Comput. Math.; 2019; 96, pp. 337-361. [DOI: https://dx.doi.org/10.1080/00207160.2018.1438605]

25. Serna-Reyes, A.J.; Macías-Díaz, J.E.; Reguera, N. A Convergent Three-Step Numerical Method to Solve a Double-Fractional Two-Component Bose–Einstein Condensate. Mathematics; 2021; 9, 1412. [DOI: https://dx.doi.org/10.3390/math9121412]

26. Ertik, H.; Şirin, H.; Demirhan, D.; Büyükkiliç, F. Fractional mathematical investigation of Bose–Einstein condensation in dilute ⁸⁷Rb, ²³Na and ⁷Li atomic gases. Int. J. Mod. Phys. B; 2012; 26, 1250096. [DOI: https://dx.doi.org/10.1142/S0217979212500968]

27. Perez-Garcia, V.M.; Michinel, H.; Cirac, J.; Lewenstein, M.; Zoller, P. Dynamics of Bose-Einstein condensates: Variational solutions of the Gross–Pitaevskii equations. Phys. Rev. A; 1997; 56, 1424. [DOI: https://dx.doi.org/10.1103/PhysRevA.56.1424]

28. Ortigueira, M.D. Riesz potential operators and inverses via fractional centred derivatives. Int. J. Math. Math. Sci.; 2006; 2006, 048391. [DOI: https://dx.doi.org/10.1155/IJMMS/2006/48391]

29. Wang, X.; Liu, F.; Chen, X. Novel second-order accurate implicit numerical methods for the Riesz space distributed-order advection-dispersion equations. Adv. Math. Phys.; 2015; 2015, 590435. [DOI: https://dx.doi.org/10.1155/2015/590435]

30. Desplanques, J. Théoreme d’algébre. J. Math. Spec.; 1887; 9, pp. 12-13.

31. Pen-Yu, K. Numerical methods for incompressible viscous flow. Sci. Sin.; 1977; 20, pp. 287-304.

32. Mickens, R.E. Nonstandard Finite Difference Schemes: Methodology And Applications; 1st ed. World Scientific: New York, NY, USA, 2020.

33. Xiao, A.; Wang, C.; Wang, J. Conservative linearly-implicit difference scheme for a class of modified Zakharov systems with high-order space fractional quantum correction. Appl. Numer. Math.; 2019; 146, pp. 379-399. [DOI: https://dx.doi.org/10.1016/j.apnum.2019.07.019]

34. Wang, Y.; Li, Q.; Mei, L. A linear, symmetric and energy-conservative scheme for the space-fractional Klein–Gordon–Schrödinger equations. Appl. Math. Lett.; 2019; 95, pp. 104-113. [DOI: https://dx.doi.org/10.1016/j.aml.2019.03.032]