A Mass- and Energy-Conserving Numerical Model for a Fractional Gross–Pitaevskii System in Multiple Dimensions

Abstract

This manuscript studies a double fractional extended p-dimensional coupled Gross–Pitaevskii-type system. This system consists of two parabolic partial differential equations with equal interaction constants, coupling terms, and spatial derivatives of the Riesz type. Associated with the mathematical model, there are energy and non-negative mass functions which are conserved throughout time. Motivated by this fact, we propose a finite-difference discretization of the double fractional Gross–Pitaevskii system which inherits the energy and mass conservation properties. As the continuous model, the mass is a non-negative constant and the solutions are bounded under suitable numerical parameter assumptions. We prove rigorously the existence of solutions for any set of initial conditions. As in the continuous system, the discretization has a discrete Hamiltonian associated. The method is implicit, multi-consistent, stable and quadratically convergent. Finally, we implemented the scheme computationally to confirm the validity of the mass and energy conservation properties, obtaining satisfactory results.

1. Introduction

Fractional calculus emerged as a generalization of conventional calculus. The classical differential and integral operators proposed by Leibniz more than three centuries ago have been extended to the fractional scenario in various different forms. In the way, various applications have been proposed to within mathematics, such as Cauchy problems with Caputo Hadamard fractional derivatives [1], the synchronization for a class of fractional-order hyperchaotic systems [2], the inequality estimates for the boundedness of multilinear singular and fractional integral operators [3], the analysis of unstructured-mesh Galerkin finite element method for the two-dimensional multi-term time–space fractional Bloch–Torrey equations on irregular convex domains [4], and the design of numerically efficient and conservative model for a Riesz space-fractional Klein–Gordon–Zakharov system [5], among various other interesting problems [6].

It is worth pointing out that each fractional operator is fully characterized by its own special kernel, and they can be used in a range of particular problems. In this respect, the analysis on the uniqueness of solutions for fractional-order differential equations may be carried out through the use of integral inequalities of fractional order, which obviously depend on the type of operator. The literature provides accounts of many applications potential applications, such as the periodic orbit analysis for the delayed Filippov systems [7], the bifurcation of limit cycles at infinity in piecewise polynomial systems [8], and the number and stability of limit cycles for planar piecewise linear systems of node–saddle type [9]. Indeed, these inequalities are often used in the theory of applied mathematics and differential equations. More well-known applications of integral inequalities are found in applied sciences, such as statistical problems, transform theory, numerical quadrature, and probability [10,11].

In the last few years, many researchers have established various types of integral inequalities by employing different approaches. For example, there are reports on some new inequalities for $(k,s)$-fractional integrals [12], certain weighted integral inequalities involving the fractional hypergeometric operators [13], and some generalized integral inequalities for Hadamard fractional integrals [14]. Moreover, those results provide direct applications to other areas of mathematics, such as conformable fractional integral inequalities [15], including those of the Hermite–Hadamard type [16]. Furthermore, there are reports in the literature on integral inequalities for a fractional operator of a function with respect to another function with nonsingular kernel [17] and Hermite–Hadamard-type inequalities for Riemann–Liouville fractional integrals [18].

Let $p\in \mathbb{N}$ and $T\in {\mathbb{R}}^{+}$. Define the set $I_n=\{1,\dots ,n\}$, for each $n\in \mathbb{N}$, and let ${\overline{I}}_n=I_n\cup \left\{0\right\}$. Suppose that $a_i,b_i\in \mathbb{R}$ satisfy $a_i<b_i$, for each $i\in I_p$. Throughout this work, we let $\mathsf{\Omega }={\mathsf{\Pi }}_{i=1}^p(a_i,b_i)\subseteq {\mathbb{R}}^p$ and ${\mathsf{\Omega }}_T=\mathsf{\Omega }\times (0,T)\subseteq {\mathbb{R}}^{p+1}$, and we assume that the sets $\overline{\mathsf{\Omega }}$ and ${\overline{\mathsf{\Omega }}}_T$ represent, respectively, the minimum closed sets containing $\mathsf{\Omega }$ and ${\mathsf{\Omega }}_T$ under the standard topology of ${\mathbb{R}}^{p+1}$, and let $\partial \mathsf{\Omega }$ be the boundary of $\mathsf{\Omega }$. Moreover, we assume that the functions ${\psi }_1:{\overline{\mathsf{\Omega }}}_T\to \mathbb{C}$ and ${\psi }_2:{\overline{\mathsf{\Omega }}}_T\to \mathbb{C}$ are sufficiently smooth. In general, if $\psi :{\overline{\mathsf{\Omega }}}_T\to \mathbb{C}$ is any function, then we extend its domain to all of ${\mathbb{R}}^p\times [0,T]$ by letting $\psi (x,t)=0$, for each $(x,t)\in ({\mathbb{R}}^p\backslash \overline{\mathsf{\Omega }})\times [0,T]$.

Definition 1

(Podlubny [19]). Suppose that $\psi :{\overline{\mathsf{\Omega }}}_T\to \mathbb{C}$, let $\alpha \in \mathbb{R}$ satisfy $\alpha \in (-1,\infty )$, and let $n\in \mathbb{Z}$ be such that the following hold: $n-1<\alpha \le n$. The partial derivative in the sense of Riesz of the function ψ with respect to $x_i$ of order α at $\left(x,t\right)\in {\mathsf{\Omega }}_T$ is

(1)$\begin{array}{cc}\hfill \frac{{\partial }^{\alpha }\psi }{\partial {\left|x_i\right|}^{\alpha }}(x,t)& =\frac{-1}{2cos\left(\frac{\pi \alpha }{2}\right)\mathsf{\Gamma }\left(n-\alpha \right)}\frac{{\partial }^n}{\partial x_i^n}{\int }_{-\infty }^{\infty }\frac{\psi (x_1,\dots ,x_{i-1},\xi ,x_{i+1},\dots ,x_p,t)}{|x_i{-\xi |}^{\alpha -1}}d\xi .\hfill \end{array}$

Here, $x=(x_1,x_2,...,x_p)\in \mathsf{\Omega }$. If $\alpha =n\in \mathbb{N}\cup \left\{0\right\}$, then the Riesz partial derivative of ψ of order α with respect to $x_i$ coincides with the classical partial derivative of ψ of order n. The Riesz fractional Laplacian and gradient operators of order α are defined, respectively, by

(2) $\begin{array}{cc}\hfill {△}^{\alpha }\psi (x,t)& =\sum _{i=1}^p\frac{{\partial }^{\alpha }\psi (x,t)}{\partial {\left|x_i\right|}^{\alpha }},\forall (x,t)\in {\mathsf{\Omega }}_T,\hfill \end{array}$

(3) $\begin{array}{cc}\hfill {▽}^{\alpha }\psi (x,t)& =\left(\frac{{\partial }^{\alpha }\psi (x,t)}{\partial |x_1{|}^{\alpha }},\frac{{\partial }^{\alpha }\psi (x,t)}{\partial |x_2{|}^{\alpha }},\dots ,\frac{{\partial }^{\alpha }\psi (x,t)}{\partial |x_p{|}^{\alpha }}\right),\forall (x,t)\in {\mathsf{\Omega }}_T.\hfill \end{array}$

For the remainder, $D,{\beta }_{11},{\beta }_{12},{\beta }_{22}\in {\mathbb{R}}^{+}\cup \left\{0\right\}$, $\lambda \in \mathbb{R}$ and $V:\overline{\mathsf{\Omega }}\to \mathbb{R}$. Suppose that ${\alpha }_1,{\alpha }_2\in (1,2]$, and let ${\varphi }_1:\overline{\mathsf{\Omega }}\to \mathbb{C}$ and ${\varphi }_2:\overline{\mathsf{\Omega }}\to \mathbb{C}$ be two functions which physically describe initial profiles for ${\psi }_1$ and ${\psi }_2$, respectively. The purpose of this work is to investigate a Riesz space-fractional two-coupled Gross–Pitaevskii system in dimensionless form, which is given by the p-dimensional initial-boundary-value problem

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

This system is particularly interesting due to the fact that it has quantities which are conserved through time, one of them being the energy function. More precisely, the total energy at time $t\in [0,T]$ is provided as

(5)$\begin{array}{cc}\hfill \mathcal{E}(t)& ={\int }_{\overline{\mathsf{\Omega }}}\left[\frac{1}{2}\sum _{i=1}^p\left(\frac{{\partial }^{{\alpha }_1}{\psi }_1}{\partial |x_i{|}^{{\alpha }_1}}{\overline{\psi }}_1+\frac{{\partial }^{{\alpha }_2}{\psi }_2}{\partial |x_i{|}^{{\alpha }_2}}{\overline{\psi }}_2\right)+V(x)\left({\left|{\psi }_1\right|}^2+{\left|{\psi }_2\right|}^2\right)+D{\left|{\psi }_1\right|}^2\right.\hfill \\ \hfill & \left.+\frac{1}{2}{\beta }_{11}{\left|{\psi }_1\right|}^4+\frac{1}{2}{\beta }_{22}{\left|{\psi }_2\right|}^4+{\beta }_{12}{\left|{\psi }_1\right|}^2{\left|{\psi }_2\right|}^2+2\lambda Re\left({\psi }_1{\overline{\psi }}_2\right)\right]dx.\hfill \end{array}$

Definition 2.

For each $z\in \mathbb{C}$, $\overline{z}$ represents the complex conjugate of z. Let $L_{x,2}(\overline{\mathsf{\Omega }})$ denote the set of all functions $f:\overline{\mathsf{\Omega }}\to \mathbb{C}$ with the property that, for each $t\in [0,T]$, the function $f(·,t):\overline{\mathsf{\Omega }}\to \mathbb{C}$ belongs to $L_2(\overline{\mathsf{\Omega }})$. For each pair $f,g\in L_{x,2}(\overline{\mathsf{\Omega }})$, the inner product of f and g is the function of t defined by

(6) ${⟨f,g⟩}_x={\int }_{\overline{\mathsf{\Omega }}}f(x,t)\overline{g(x,t)}dx,\forall t\in [0,T].$

In addition, if $u=(u_1,u_2,\dots ,u_p)$ and $v=(v_1,v_2,\dots ,v_p)$ are vector functions such that $u_i,v_i\in L_{x,2}(\overline{\mathsf{\Omega }})$ for each $i\in I_p$, then we set ${⟨u,v⟩}_x={⟨u_1,v_1⟩}_x+{⟨u_2,v_2⟩}_x+\dots +{⟨u_p,v_p⟩}_x$. Let ${\parallel f\parallel }_{x,2}=\sqrt{⟨f,f⟩}$, for each $f\in L_{x,2}(\overline{\mathsf{\Omega }})$. In similar fashion, we can define the space $L_{x,p}(\overline{\mathsf{\Omega }})$. If $f\in L_{x,p}(\overline{\mathsf{\Omega }})$, then let

(7) ${\parallel f\parallel }_{x,p}={\left({\int }_{\overline{\mathsf{\Omega }}}\left|f(x,t)\right|dx\right)}^{1/p},\forall t\in [0,T].$

The next result summarizes the property of conservation of energy satisfied by (4).

Theorem 1

(Energy conservation). The energy function (5) can be rewritten alternatively as

(8)$\begin{array}{cc}\hfill \mathcal{E}(t)& =\frac{1}{2}\left({∥{▽}^{{\alpha }_1/2}{\psi }_1∥}_{x,2}^2+{∥{▽}^{{\alpha }_2/2}{\psi }_2∥}_{x,2}^2\right)+{⟨V(x),{\left|{\psi }_1\right|}^2+{\left|{\psi }_2\right|}^2⟩}_x+D{\parallel {\psi }_1\parallel }_{x,2}^2\hfill \\ \hfill & +\frac{1}{2}{\beta }_{11}\parallel {\psi }_1{\parallel }_{x,4}^4+\frac{1}{2}{\beta }_{22}\parallel {\psi }_2{\parallel }_{x,4}^4+{\beta }_{12}{\parallel {\psi }_1{\psi }_2\parallel }_{x,2}^2+2\lambda Re{⟨{\psi }_1,{\psi }_2⟩}_x,\hfill \end{array}$

for each $t\in (0,T)$. Moreover, if ${\psi }_1$ and ${\psi }_2$ satisfy (4), then the function $\mathcal{E}(t)$ is a constant.

Proof. 

The proof of this result was given by Serna-Reyes and Macías-Díaz [20]. □

The energy density of the system (4) is defined as the function in the integrand of the expression of $\mathcal{E}(t)$. More precisely, the energy density of the physical model (4) is given by the function $\mathcal{H}({\psi }_1,{\psi }_2)=\mathcal{H}(x,t)$, where, for each $(x,t)\in {\mathsf{\Omega }}_T$,

(9)$\begin{array}{cc}\hfill \mathcal{H}(x,t)& {\displaystyle =\frac{1}{2}\sum _{i=1}^p\left({\left|\frac{{\partial }^{{\alpha }_1/2}{\psi }_1}{\partial |x_i{|}^{{\alpha }_1/2}}\right|}^2+{\left|\frac{{\partial }^{{\alpha }_2/2}{\psi }_2}{\partial |x_i{|}^{{\alpha }_2/2}}\right|}^2\right)+V(x)\left({\left|{\psi }_1\right|}^2+{\left|{\psi }_2\right|}^2\right)+D{\left|{\psi }_1\right|}^2}\hfill \\ \hfill & {\displaystyle +\frac{1}{2}{\beta }_{11}{\left|{\psi }_1\right|}^4+\frac{1}{2}{\beta }_{22}{\left|{\psi }_2\right|}^4+{\beta }_{12}{\left|{\psi }_1\right|}^2{\left|{\psi }_2\right|}^2+2\lambda Re\left({\psi }_1{\overline{\psi }}_2\right).}\hfill \end{array}$

On the other hand, the total mass of the continuous system (4) at the time t is defined as $\mathcal{M}(t)={\mathcal{M}}_1(t)+{\mathcal{M}}_2(t)$, for each $t\in (0,T)$. Here, we let ${\mathcal{M}}_1(t)={\parallel {\psi }_1\parallel }_{x,2}^2$ and ${\mathcal{M}}_2(t)={\parallel {\psi }_2\parallel }_{x,2}^2$, for each $t\in (0,T)$.

Theorem 2

(Mass conservation). If ${\psi }_1$ and ${\psi }_2$ are solutions of (4), then the function $\mathcal{M}(t)$ is non-negative and constant. Moreover, if $\lambda =0$, then the functions ${\mathcal{M}}_1(t)$ and ${\mathcal{M}}_2(t)$ are also non-negative and constant.

Proof. 

Beforehand, recall that the square-root operator of $-{△}^{\alpha }$ is the operator ${▽}^{\alpha /2}$, for each $\alpha \in (1,2]$. This means that ${⟨-{△}^{\alpha }u,v⟩}_x={⟨{▽}^{\alpha /2}u,{▽}^{\alpha /2}v⟩}_x$, for any functions $u,v\in H^2({\overline{\mathsf{\Omega }}}_T)$. Notice now that the following identities hold:

(10)$\begin{array}{cc}\hfill \mathrm{i}{⟨\frac{\partial {\psi }_k}{\partial t},{\psi }_k⟩}_x& =\mathrm{i}\frac{d}{dt}{\parallel {\psi }_k\parallel }_{x,2}^2,\forall t\in (0,T),\forall k\in I_2,\hfill \end{array}$

(11)$\begin{array}{cc}\hfill {⟨-{△}^{{\alpha }_k}{\psi }_k,{\psi }_k⟩}_x& ={∥{▽}^{{\alpha }_k/2}{\psi }_k∥}_{x,2}^2,\forall t\in (0,T),\forall k\in I_2,\hfill \end{array}$

(12)$\begin{array}{cc}\hfill {⟨(V+D){\psi }_1,{\psi }_1⟩}_x& ={⟨V+D,|{\psi }_1{|}^2⟩}_x,\forall t\in (0,T),\hfill \end{array}$

(13)$\begin{array}{cc}\hfill {⟨V{\psi }_2,{\psi }_2⟩}_x& ={⟨V,|{\psi }_2{|}^2⟩}_x,\forall t\in (0,T),\hfill \end{array}$

(14)$\begin{array}{cc}\hfill {⟨({\beta }_{11}|{\psi }_1{|}^2+{\beta }_{12}|{\psi }_2{|}^2){\psi }_1,{\psi }_1⟩}_x& ={⟨{\beta }_{11}|{\psi }_1{|}^2+{\beta }_{12}|{\psi }_2{|}^2,{\left|{\psi }_1\right|}^2⟩}_x,\forall t\in (0,T),\hfill \end{array}$

(15)$\begin{array}{cc}\hfill {⟨({\beta }_{12}|{\psi }_1{|}^2+{\beta }_{22}|{\psi }_2{|}^2){\psi }_2,{\psi }_2⟩}_x& ={⟨{\beta }_{12}|{\psi }_1{|}^2+{\beta }_{22}|{\psi }_2{|}^2,{\left|{\psi }_2\right|}^2⟩}_x,\forall t\in (0,T).\hfill \end{array}$

These inner products are real numbers, except the first one which is purely imaginary. On both sides of the first equation of (4), obtain the inner product with ${\psi }_1$ and take the imaginary parts. In addition, calculate the inner product on both sides of the second differential equation of (4) with ${\psi }_2$, and take also the imaginary parts on both sides. Then,

(16)$\begin{array}{cc}\hfill \frac{d{\mathcal{M}}_1(t)}{dt}& =\lambda Im\overline{{⟨{\psi }_1,{\psi }_2⟩}_x},\forall t\in (0,T),\hfill \end{array}$

(17)$\begin{array}{cc}\hfill \frac{d{\mathcal{M}}_2(t)}{dt}& =\lambda Im{⟨{\psi }_1,{\psi }_2⟩}_x,\forall t\in (0,T).\hfill \end{array}$

Notice that the individual masses ${\mathcal{M}}_1(t)$ and ${\mathcal{M}}_2(t)$ are constant when the internal atomic Josephson junction is equal to zero, as mentioned in the conclusion of this result. Moreover, if we add these two identities, it readily follows that

(18)$\frac{d\mathcal{M}(t)}{dt}=\lambda Im\left({⟨{\psi }_1,{\psi }_2⟩}_x+\overline{{⟨{\psi }_1,{\psi }_2⟩}_x}\right)=2\lambda Im\left(Re{⟨{\psi }_1,{\psi }_2⟩}_x\right)=0,\forall t\in (0,T).$

Conclude that the total mass of the system is conserved throughout time, as desired. □

We introduce the function of mass density of the continuous fractional model (4) as $\mathcal{R}({\psi }_1,{\psi }_2)=\mathcal{R}(x,t)={\mathcal{R}}_1(x,t)+{\mathcal{R}}_2(x,t)$, for each $(x,t)\in {\mathsf{\Omega }}_T$. In these expressions, ${\mathcal{R}}_1(x,t)={\left|{\psi }_2(x,t)\right|}^2$ and ${\mathcal{R}}_2(x,t)={\left|{\psi }_2(x,t)\right|}^2$, for each $(x,t)\in {\mathsf{\Omega }}_T$. It is obvious that

(19)$\begin{array}{cc}\hfill \mathcal{E}(t)& ={\int }_{\mathsf{\Omega }}\mathcal{H}(x,t)dx,\forall t\in [0,T],\hfill \end{array}$

(20)$\begin{array}{cc}\hfill \mathcal{M}(t)& ={\int }_{\mathsf{\Omega }}\mathcal{R}(x,t)dx,\forall t\in [0,T].\hfill \end{array}$

The following is an easy consequence of Theorem 2.

Corollary 1

(Boundedness). If ${\psi }_1$ and ${\psi }_2$ are solutions of the (4), then there exists a constant $M_0\ge 0$ such that $max\left\{\parallel {\psi }_1{\parallel }_2^2,{\parallel {\psi }_2\parallel }_2^2\right\}\le M_0$, for each $t\in (0,T)$. □

In the present work, we design a finite-difference method to solve the system (4), in such way that the main structural properties of the solutions of that continuous system are also satisfied in the discrete domain. More precisely, we design a numerical method which is able of preserving the discrete energy and mass of the system, and such that the solutions satisfy similar boundedness properties. The scheme is thoroughly analyzed for consistency, stability, and consistency, and we illustrate these features through some computer simulations. As it turns out, the proposed numerical model has an easy implementation, as shown by the code in Appendix A.

2. Numerical Model

In the following, we propose a linearly-implicit numerical method to approximate the solutions of the fractional system (4), along with the energy and mass functions. To that end, the concept of fractional centered differences is of utmost importance.

Definition 3

(Ortigueira [21]). Let $f:\mathbb{R}\to \mathbb{R}$ be any function, $h>0$, and $\alpha >-1$. The fractional-order centered difference of f of order α at x is defined as

(21)${\Delta }_h^{(\alpha )}f(x)=\sum _{k=-\infty }^{\infty }{g}_k^{(\alpha )}f(x-kh),\forall x\in \mathbb{R},$

where

(22)$g_k^{(\alpha )}=\frac{{(-1)}^k\mathsf{\Gamma }(\alpha +1)}{\mathsf{\Gamma }(\frac{\alpha }{2}-k+1)\mathsf{\Gamma }(\frac{\alpha }{2}+k+1)},\forall k\in \mathbb{Z}.$

Lemma 1

(Wang et al. [22]). Let $0<\alpha \le 2$ and $\alpha \ne 1$. Then,

(23)$\begin{array}{cc}\hfill g_0^{(\alpha )}& >0,\hfill \end{array}$

(24)$\begin{array}{cc}\hfill g_k^{(\alpha )}& =g_{-k}^{(\alpha )}<0,\forall k\ne 0,\hfill \end{array}$

(25)$\begin{array}{cc}\hfill \sum _{k=-\infty }^{\infty }{g}_k^{(\alpha )}& =0.\hfill \end{array}$

Lemma 2

(Wang et al. [22]). Let $0<\alpha \le 2$ and $\alpha \ne 1$. If $f:\mathbb{R}\to \mathbb{R}$ has integrable derivatives up to order 5, then, for almost all x,

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

Next, we present the discrete notation used to approximate the solutions of the fractional problem (4). We let $h_i,\tau \in \mathbb{R}$, for each $i\in I_p$. Moreover, we assume that $N=T/\tau$ and $M_i=(b_i-a_i)/h_i$ are natural numbers, for each $i\in I_p$, and fix uniform partitions of $[a_i,b_i]$ and $[0,T]$ given, respectively, by

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

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

Let $J={\prod }_{i=1}^p{I}_{M_i-1}$ and $\overline{J}={\prod }_{i=1}^p{\overline{I}}_{M_i}$, and let $\partial J$ represent the boundary of the mesh $\overline{J}$. Define $x_j=(x_{1,j_1},\dots ,x_{p,j_p})$, for each multi-index $j=(j_1,\dots ,j_p)\in \overline{J}$. The symbol $(u_j^n,v_j^n)$ denotes a numerical estimate to $({\psi }_1(x_j,t_n),{\psi }_2(x_j,t_n))$, for each $(j,n)\in \overline{J}\times {\overline{I}}_N$.

Definition 4.

Define the following, for $w=u,v$ and $\alpha \in (0,1)\cup (1,2]$, $(j,n)\in J\times I_{N-1}$:

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

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

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

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

and

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

Moreover, for the sake of convenience, we agree that

(34) $\begin{array}{cc}\hfill {△}_h^{(\alpha )}{w}_j^n& ={\delta }_{x_1}^{(\alpha )}{w}_j^n+{\delta }_{x_2}^{(\alpha )}{w}_j^n+\dots +{\delta }_{x_p}^{(\alpha )}{w}_j^n,\forall (j,n)\in J\times I_{N-1},\hfill \end{array}$

(35) $\begin{array}{cc}\hfill {▽}_h^{(\alpha )}{w}_j^n& =({\delta }_{x_1}^{(\alpha )}{w}_j^n,{\delta }_{x_2}^{(\alpha )}{w}_j^n,\dots ,{\delta }_{x_p}^{(\alpha )}{w}_j^n),\forall (j,n)\in J\times I_{N-1}.\hfill \end{array}$

The finite-difference method to approximate the solution of (4) on ${\mathsf{\Omega }}_T$ is the discrete system with initial-boundary conditions

(36)$\begin{array}{c}\begin{array}{cc}\hfill \mathrm{i}{\delta }_t^{(1)}{u}_j^n& =\left(-\frac{1}{2}{△}_h^{({\alpha }_1)}+V_j+D+{\beta }_{11}{\left|u_j^n\right|}^2+{\beta }_{12}{\left|v_j^n\right|}^2\right){\mu }_t^{(1)}{u}_j^n+\lambda v_j^n,\hfill \\ \hfill \mathrm{i}{\delta }_t^{(1)}{v}_j^n& =\left(-\frac{1}{2}{△}_h^{({\alpha }_2)}+V_j+{\beta }_{22}{\left|v_j^n\right|}^2+{\beta }_{12}{\left|u_j^n\right|}^2\right){\mu }_t^{(1)}{v}_j^n+\lambda u_j^n,\hfill \end{array}\\ \mathrm{such}\mathrm{that}\left\{\begin{array}{cc}{u}_j^0={\mu }_t^{(1)}{u}_j^0={\varphi }_1(x_j),\hfill & \forall j\in \overline{J},\hfill \\ v_j^0={\mu }_t^{(1)}{v}_j^0={\varphi }_2(x_j),\hfill & \forall j\in \overline{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}$

where $(j,n)\in J\times I_{N-1}$. Notice that this scheme is an implicit three-step method. Moreover, the approximations at the time $t_0$ are provided by the initial data $u_j^0={\varphi }_1(x_j)$ and $v_j^0={\varphi }_2(x_j)$, for each $j\in \overline{J}$. On the other hand, the discrete model requires knowledge of fictitious approximations at the time $t_{-1}$. To avoid this shortcoming, notice that the initial conditions ${\mu }_t^{(1)}{u}_j^0={\varphi }_1(x_j)$ and ${\mu }_t^{(1)}{v}_j^0={\varphi }_2(x_j)$ show that $u_j^{-1}=2{\varphi }_1(x_j)-u_j^1$ and $v_j^{-1}=2{\varphi }_2(x_j)-v_j^1$, respectively, for each $j\in J$. Substituting these identities into the recursive difference equations of (36) when $n=0$, we obtain the explicit formulas

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

and

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

For the sake of convenience, the stencil of the scheme (36) is provided in Figure 1.

We study the existence and uniqueness of solutions of the finite-difference method (36). To that end, we require some additional nomenclature and technical results. To start with, we let $h=(h_1,\dots ,h_p)$ and $h_{\ast }={\prod }_{i=1}^p{h}_i$, and fix the spatial mesh ${\left\{x_j\right\}}_{j\in J}\subseteq {\mathbb{R}}^p$. Let ${\mathcal{V}}_h$ be the complex vector space of all grid functions which vanish at the boundary of the mesh. For any $w\in {\mathcal{V}}_h$ and $j\in J$, convey that $w_j=w(x_j)$.

Definition 5.

Define the product $⟨·,·⟩:{\mathcal{V}}_h\times {\mathcal{V}}_h\to \mathbb{C}$ and the norm ${\parallel ·\parallel }_1:{\mathcal{V}}_h\to \mathbb{R}$ by

(39) $\begin{array}{cc}\hfill ⟨u,v⟩& =h_{\ast }\sum _{j\in I}{u}_j\overline{v_j},\forall u,v\in {\mathcal{V}}_h,\hfill \end{array}$

(40) $\begin{array}{cc}\hfill {\parallel u\parallel }_1& =h_{\ast }\sum _{j\in I}\left|u_j\right|,\forall u,v\in {\mathcal{V}}_h.\hfill \end{array}$

The Euclidean norm induced by the inner product $⟨·,·⟩$ is denoted by ${\parallel ·\parallel }_2$.

In the following, we represent the solutions of the finite-difference method (36) by ${(u^n,v^n)}_{n=0}^N$, where we agree that $u^n={(u_j^n)}_{j\in J}$ and $v^n={(v_j^n)}_{j\in J}$, for each $n\in {\overline{I}}_N$. We use the following lemma to prove that solutions of (36) exist.

Lemma 3

(Browder fixed-point theorem [23]). Let $(H,⟨·,·⟩)$ be a finite-dimensional inner-product space, let $\parallel ·\parallel :H\to H$ be the norm induced by the inner product, and suppose that $g:H\to H$ is continuous. Assume that there exists $\beta >0$ such that $Re⟨g(z),z⟩>0$, for all $z\in H$ with $\parallel z\parallel =\beta$. Then, there is $z^{\ast }\in H$, such that $g(z^{\ast })=0$ and $\parallel z^{\ast }\parallel \le \beta$.

Theorem 3

(Solubility). For any set of initial conditions, the discrete system (36) is solvable.

Proof. 

The proof requires mathematical induction. Notice that the approximations at the times $t_0$ and $t_1$ are provided by the initial profiles and the explicit systems (37) and (38), respectively. Thus, let us suppose that the numerical solutions of the discrete model (36) at the times $n-1$ and n are known, for some $n\in {\overline{I}}_{N-1}$. Rearranging the recursive equations of the method, we can readily check that the following equations hold:

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

and

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

for each $j\in J$. Define the function $G:{\mathcal{V}}_h\times {\mathcal{V}}_h\to {\mathcal{V}}_h\times {\mathcal{V}}_h$ by $G(\eta ,\nu )=(G_1(\eta ),G_2(\nu ))$, for each $(\eta ,\nu )\in {\mathcal{V}}_h\times {\mathcal{V}}_h$, where the functions $G_1,G_2:{\mathcal{V}}_h\to {\mathcal{V}}_h$ are defined by

(43)$\begin{array}{cc}\hfill {\left[G_1(\eta )\right]}_j& ={\eta }_j^n-u_j^{n-1}+\mathrm{i}\tau \left[-\frac{1}{2}{△}_h^{({\alpha }_1)}+V_j+D+{\beta }_{11}|u_j^n{|}^2+{\beta }_{12}{\left|v_j^n\right|}^2\right]{\eta }_j^n{+\mathrm{i}\lambda \tau v|}_j^n,\hfill \end{array}$

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

for each $\eta \in {\mathcal{V}}_h$, $\nu \in {\mathcal{V}}_h$, $j\in J$. After some calculations, we readily reach that

(45)$\begin{array}{cc}\hfill Re⟨G_1(\eta ),\eta ⟩& ={\parallel \eta \parallel }^2-Re⟨u^{n-1},\eta ⟩+Re\left(\mathrm{i}\lambda \tau ⟨v^n,\eta ⟩\right),\forall \eta \in {\mathcal{V}}_h,\hfill \end{array}$

(46)$\begin{array}{cc}\hfill Re⟨G_2(\nu ),\nu ⟩& ={\parallel \nu \parallel }^2-Re⟨v^{n-1},\nu ⟩+Re\left(\mathrm{i}\lambda \tau ⟨u^n,\nu ⟩\right),\forall \nu \in {\mathcal{V}}_h.\hfill \end{array}$

Using now the Cauchy–Schwarz inequality and bounding from below, it is easy to check that the following hold:

(47)$\begin{array}{cc}\hfill Re⟨G(\eta ,\nu ),(\eta ,\nu )⟩& {\displaystyle \ge {\parallel (\eta ,\nu )\parallel }^2-|⟨(u^{n-1},v^{n-1}),(\eta ,\nu )⟩|+|\lambda \tau |·|⟨(v^n,u^n),(\eta ,\nu )⟩|}\hfill \\ \hfill & {\displaystyle \ge \parallel (\eta ,\nu )\parallel \left(\parallel (\eta ,\nu )\parallel -\parallel (u^{n-1},v^{n-1})\parallel -\lambda \tau \parallel (v^n,u^n)\parallel \right).}\hfill \end{array}$

If we let $\beta =\parallel (u^{n-1},v^{n-1})\parallel +\lambda \tau \parallel (v^n,u^n)\parallel +1$ in the previous lemma, it follows that there exists $(u^{n+1},v^{n+1})\in {\mathcal{V}}_h\times {\mathcal{V}}_h$ which satisfies the discrete system (36). The conclusion of this result readily follows now by induction. □

3. Structural Properties

The present section is devoted to showing that the discrete model (36) satisfies physical properties which are similar to those satisfied by the continuous system (4). More precisely, we propose numerical energy and mass functional associated to the scheme (36) which are conserved conditionally. Various technical lemmas are required to that end.

Lemma 4

(Macías-Díaz [24]). For each $i\in I_p$ and $k\in I_2$, there exists a unique positive operator ${\delta }_{x_i}^{{\alpha }_k/2}:{\mathcal{V}}_h\to {\mathcal{V}}_h$ such that $⟨-{\delta }_{x_i}^{{\alpha }_k}u,v⟩=⟨{\delta }_{x_i}^{{\alpha }_k/2}u,{\delta }_{x_i}^{{\alpha }_k/2}v⟩$, for each $u,v\in {\mathcal{V}}_h$.

Lemma 5.

The following equations hold, for each $n\in I_{N-1}$:

(48) $\begin{array}{cc}\hfill Re⟨\mathrm{i}{\delta }_t^{(1)}{u}^n,2{\delta }_t^{(1)}{u}^n⟩& =0,\hfill \end{array}$

(49) $\begin{array}{c}\hfill Re⟨V{\mu }_t^{(1)}{u}^n,2{\delta }_t^{(1)}{u}^n⟩={\delta }_t⟨V,{\mu }_t|u^{n-1}{|}^2⟩,\end{array}$

(50) $\begin{array}{cc}\hfill Re⟨D{\mu }_t^{(1)}{u}^n,2{\delta }_t^{(1)}{u}^n⟩& =D{\delta }_t{\mu }_t{\parallel u^{n-1}\parallel }_2^2,\hfill \end{array}$

(51) $\begin{array}{cc}\hfill Re⟨{\beta }_{11}|u^n{|}^2{\mu }_t^{(1)}{u}^n,2{\delta }_t^{(1)}{u}^n⟩& =\frac{1}{2}{\beta }_{11}{\delta }_t⟨|u^n{|}^2,|u^{n-1}{|}^2⟩,\hfill \end{array}$

(52) $\begin{array}{cc}\hfill Re⟨{\beta }_{12}|v^n{|}^2{\mu }_t^{(1)}{u}^n,2{\delta }_t^{(1)}{u}^n⟩& =\frac{1}{2}{\beta }_{12}⟨|v^n{|}^2,{\delta }_t^{(1)}|u^n{|}^2⟩.\hfill \end{array}$

Proof. 

To establish (48), notice that the following identities trivially hold:

(53)$Re⟨\mathrm{i}{\delta }_t^{(1)}{u}^n,2{\delta }_t^{(1)}{u}^n⟩=Re\left(2\mathrm{i}\parallel {\delta }_t^{(1)}{u}^n{\parallel }_2^2\right)=0.$

To prove the identity (49), observe that

(54)$\begin{array}{cc}\hfill & Re⟨V{\mu }_t^{(1)}{u}^n,2{\delta }_t^{(1)}{u}^n⟩\hfill \\ \hfill & =\frac{1}{2\tau }Re\left(⟨V,{\left|u^{n+1}\right|}^2⟩+2\mathrm{i}Im⟨V{u}^{n-1},u^{n+1}⟩-⟨V,{\left|u^{n-1}\right|}^2⟩\right)\hfill \\ \hfill & =\frac{1}{\tau }\left[⟨V,\frac{1}{2}\left({\left|u^{n+1}\right|}^2+{\left|u^n\right|}^2\right)⟩-⟨V,\frac{1}{2}\left({\left|u^n\right|}^2+{\left|u^{n-1}\right|}^2\right)⟩\right]\hfill \\ \hfill & ={\delta }_t⟨V,{\mu }_t{\left|u^{n-1}\right|}^2⟩,\hfill \end{array}$

which establishes the second identity. We prove now the formula in part (51). To that end, notice that

(55)$\begin{array}{cc}\hfill & Re⟨{\beta }_{11}{\left|u^n\right|}^2{\mu }_t^{(1)}{u}^n,2{\delta }_t^{(1)}{u}^n⟩\hfill \\ \hfill & =\frac{{\beta }_{11}}{2\tau }Re\left(⟨{\left|u^n\right|}^2,{\left|u^{n+1}\right|}^2⟩+2\mathrm{i}Im⟨{\left|u^n\right|}^2{u}^{n-1},u^{n+1}⟩-⟨{\left|u^n\right|}^2,{\left|u^{n-1}\right|}^2⟩\right)\hfill \\ \hfill & =\frac{{\beta }_{11}}{2}\left[\frac{1}{\tau }\left(⟨{\left|u^{n+1}\right|}^2,{\left|u^n\right|}^2⟩-⟨{\left|u^n\right|}^2,{\left|u^{n-1}\right|}^2⟩\right)\right]=\frac{{\beta }_{11}}{2}{\delta }_t⟨{\left|u^n\right|}^2,{\left|u^{n-1}\right|}^2⟩.\hfill \end{array}$

Finally, we must point out that the proofs for the identities (50) and (52) are similar to the proofs of (49) and (51), respectively, and, for that reason, we omit them here. □

Lemma 6.

The following equations hold, for each $n\in I_{N-2}$ and $k\in I_2$:

(56) $\begin{array}{cc}\hfill ⟨|v^n{|}^2,{\delta }_t^{(1)}|u^n{|}^2⟩+⟨|u^n{|}^2,{\delta }_t^{(1)}|v^n{|}^2⟩& ={\delta }_t(⟨|u^{n-1}{|}^2,|v^n{|}^2⟩+⟨|v^{n-1}{|}^2,|u^n{|}^2⟩),\hfill \end{array}$

(57) $\begin{array}{cc}\hfill Re⟨-\frac{1}{2}{\delta }_{x_i}^{({\alpha }_k)}{\mu }_t^{(1)}{u}^n,2{\delta }_t^{(1)}{u}^n⟩& =\frac{1}{2}{\delta }_t{\mu }_t{\parallel {\delta }_{x_i}^{({\alpha }_k/2)}{u}^{n-1}\parallel }_2^2,\hfill \end{array}$

(58) $\begin{array}{cc}\hfill Re(⟨v^n,2{\delta }_t^{(1)}{u}^n⟩+⟨u^n,2{\delta }_t^{(1)}{v}^n⟩)& ={\delta }_{t}Re(⟨u^{n-1},v^n⟩+⟨v^{n-1},u^n⟩).\hfill \end{array}$

Proof. 

Applying distributivity of the inner product and rearranging terms,

(59)$\begin{array}{cc}\hfill & ⟨{\left|v^n\right|}^2,{\delta }_t^{(1)}{\left|u^n\right|}^2⟩+⟨{\left|u^n\right|}^2,{\delta }_t^{(1)}{\left|v^n\right|}^2⟩\hfill \\ \hfill & =\frac{1}{\tau }\left(⟨{\left|v^n\right|}^2,{\left|u^{n+1}\right|}^2⟩-⟨{\left|v^n\right|}^2,{\left|u^{n-1}\right|}^2⟩+⟨{\left|u^n\right|}^2,{\left|v^{n+1}\right|}^2⟩\right.\hfill \\ \hfill & \left.-⟨{\left|u^n\right|}^2,{\left|v^{n-1}\right|}^2⟩\right)\hfill \\ \hfill & ={\delta }_t\left(⟨{\left|u^{n-1}\right|}^2,{\left|v^n\right|}^2⟩+⟨{\left|v^{n-1}\right|}^2,{\left|u^n\right|}^2⟩\right),\hfill \end{array}$

whence the first identity readily follows. To prove the second identity, we use once again the distributivity of the inner product and Lemma 4 to see that

(60)$\begin{array}{cc}\hfill & Re⟨-\frac{1}{2}{\delta }_{x_i}^{({\alpha }_k)}{\mu }_t^{(1)}{u}^n,2{\delta }_t^{(1)}{u}^n⟩\hfill \\ \hfill & =\frac{1}{4\tau }Re⟨{\delta }_{x_i}^{({\alpha }_k/2)}{u}^{n+1}+{\delta }_{x_i}^{({\alpha }_k/2)}{u}^{n-1},{\delta }_{x_i}^{({\alpha }_k/2)}{u}^{n+1}-{\delta }_{x_i}^{({\alpha }_k/2)}{u}^{n-1}⟩\hfill \\ \hfill & =\frac{1}{4\tau }Re\left({∥{\delta }_{x_i}^{({\alpha }_k/2)}{u}^{n+1}∥}_2^2+2\mathrm{i}Im⟨{\delta }_{x_i}^{({\alpha }_k/2)}{u}^{n+1},{\delta }_{x_i}^{({\alpha }_k/2)}{u}^{n-1}⟩\right.\hfill \\ \hfill & \left.-{∥{\delta }_{x_i}^{({\alpha }_k/2)}{u}^{n-1}∥}_2^2\right),\hfill \end{array}$

which readily establishes the identity. Finally, using again distributivity and rearranging terms conveniently,

(61)$\begin{array}{cc}\hfill & Re\left(⟨v^n,2{\delta }_t^{(1)}{u}^n⟩+⟨u^n,2{\delta }_t^{(1)}{v}^n⟩\right)\hfill \\ \hfill & =\frac{1}{\tau }Re\left(⟨u^n,v^{n+1}⟩-⟨v^n,u^{n-1}⟩+⟨v^n,u^{n+1}⟩-⟨u^n,v^{n-1}⟩\right)\hfill \\ \hfill & ={\delta }_{t}Re\left(⟨u^{n-1},v^n⟩+⟨v^{n-1},u^n⟩\right).\hfill \end{array}$

This last chain of identities readily establishes the last property of the lemma. □

The next theorem proves the existence of energy-like invariants for (36).

Theorem 4

(Energy conservation). Let ${(u^n,v^n)}_{n=0}^N$ be a solution of (36), and define

(62)$\begin{array}{cc}\hfill E^n& =\frac{1}{2}{\mu }_t\left({∥{▽}_h^{({\alpha }_1/2)}{u}^n∥}_2^2+{∥{▽}_h^{({\alpha }_2/2)}{v}^n∥}_2^2\right)+⟨V,{\mu }_t\left({\left|u^n\right|}^2+{\left|v^n\right|}^2\right)⟩\hfill \\ \hfill & +D{\mu }_t{\parallel u^n\parallel }_2^2+\frac{{\beta }_{11}}{2}⟨{\left|u^n\right|}^2,{\left|u^{n+1}\right|}^2⟩+\frac{{\beta }_{22}}{2}⟨{\left|v^n\right|}^2,{\left|v^{n+1}\right|}^2⟩\hfill \\ \hfill & +\frac{{\beta }_{12}}{2}\left(⟨{\left|u^n\right|}^2,{\left|v^{n+1}\right|}^2⟩+⟨{\left|v^n\right|}^2,{\left|u^{n+1}\right|}^2⟩\right)\hfill \\ \hfill & +\lambda Re\left(⟨u^n,v^{n+1}⟩+⟨v^n,u^{n+1}⟩\right),\hfill \end{array}$

for each $n\in I_{N-1}$. Then, ${\delta }_t{E}^n=0$, for each $n\in I_{N-1}$.

Proof. 

Denote the left-hand sides of the two difference equations in (36) by ${\mathsf{\Theta }}_j^n$ and ${\mathsf{\Phi }}_j^n$, respectively, for each $j\in J$ and $n\in I_{N-1}$, and define ${\mathsf{\Theta }}^n={({\mathsf{\Theta }}_j^n)}_{j\in J}$, ${\mathsf{\Phi }}^n={({\mathsf{\Phi }}_j^n)}_{j\in J}$. Using that ${(u^n,v^n)}_{n=0}^N$ is a solution of (36), calculating the inner product of ${\mathsf{\Theta }}^n$ with $2{\delta }_t^{(1)}{u}^n$ and of ${\mathsf{\Phi }}^n$ with $2{\delta }_t^{(1)}{v}^n$, using the identities above and collecting terms, we reach

(63)$\begin{array}{cc}\hfill 0& =Re⟨\mathrm{i}{\delta }_t^{(1)}{u}^n,2{\delta }_t^{(1)}{u}^n⟩=Re⟨{\mathsf{\Theta }}^n,2{\delta }_t^{(1)}{u}^n⟩\hfill \\ \hfill & ={\delta }_t\left[\frac{1}{2}{\mu }_t{∥{▽}_h^{({\alpha }_1/2)}{u}^{n-1}∥}_2^2+⟨V,{\mu }_t{\left|u^{n-1}\right|}^2⟩+D{\mu }_t{\parallel u^{n-1}\parallel }_2^2\right.\hfill \\ \hfill & \left.+\frac{{\beta }_{11}}{2}⟨{\left|u^n\right|}^2,{\left|u^{n-1}\right|}^2⟩+\frac{{\beta }_{12}}{2}⟨{\left|v^n\right|}^2,{\left|u^{n-1}\right|}^2⟩\right]\hfill \\ \hfill & +2\lambda Re⟨v^n,{\delta }_t^{(1)}{u}^n⟩,\forall n\in I_{N-1}.\hfill \end{array}$

In similar fashion, notice that the following hold, for each $n\in I_{N-1}$:

(64)$\begin{array}{cc}\hfill 0& =Re⟨\mathrm{i}{\delta }_t^{(1)}{v}^n,2{\delta }_t^{(1)}{v}^n⟩=Re⟨{\mathsf{\Phi }}^n,2{\delta }_t^{(1)}{v}^n⟩\hfill \\ \hfill & ={\delta }_t\left[\frac{1}{2}{\mu }_t{∥{▽}_h^{({\alpha }_2/2)}{v}^{n-1}∥}_2^2+⟨V,{\mu }_t{\left|v^{n-1}\right|}^2⟩+\frac{{\beta }_{22}}{2}⟨{\left|v^n\right|}^2,{\left|v^{n-1}\right|}^2⟩\right.\hfill \\ \hfill & \left.+\frac{{\beta }_{12}}{2}⟨{\left|u^n\right|}^2,{\left|v^{n-1}\right|}^2⟩\right]+2\lambda Re⟨u^n,{\delta }_t^{(1)}{v}^n⟩.\hfill \end{array}$

Adding these two last identities, we readily obtain that ${\delta }_t{E}^{n-1}=0$, for each $n\in {\overline{I}}_{N-1}$. The conclusion of this result is readily reached now using mathematical induction. □

Motivated by this result, the quantities $E^n$ are called the discrete total energy of the system (36) at the time $t_n$, for each $n\in {\overline{I}}_{N-1}$. In addition, we define the discrete energy density of the system (36) as

(65)$\begin{array}{cc}\hfill H_j^n& =H(u_j^n,v_j^n)\hfill \\ \hfill & =\frac{1}{2}{\mu }_t\sum _{i=1}^p\left({\left|{\delta }_{x_i}^{({\alpha }_1/2)}{u}_j^n\right|}^2+{\left|{\delta }_{x_i}^{({\alpha }_2/2)}{v}_j^n\right|}^2\right)+V_j{\mu }_t\left({\left|u_j^n\right|}^2+{\left|v_j^n\right|}^2\right)+D{\mu }_t{\left|u_j^n\right|}^2\hfill \\ \hfill & +\frac{{\beta }_{11}}{2}{\left|u_j^n\right|}^2{\left|u_j^{n+1}\right|}^2+\frac{{\beta }_{22}}{2}{\left|v_j^n\right|}^2{\left|v_j^{n+1}\right|}^2\hfill \\ \hfill & +\frac{{\beta }_{12}}{2}\left({\left|u_j^n\right|}^2{\left|v_j^{n+1}\right|}^2+{\left|v_j^n\right|}^2{\left|u_j^{n+1}\right|}^2\right)+\lambda Re\left(u_j^n{\overline{v}}_j^{n+1}+v_j^n{\overline{u}}_j^{n+1}\right),\hfill \end{array}$

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

Lemma 7.

The following equations hold, for each $n\in I_{N-2}$ and $k\in I_2$:

(66) $\begin{array}{cc}\hfill Im⟨-\frac{1}{2}{\delta }_{x_i}^{({\alpha }_k)}{\mu }_t^{(1)}{u}^n,2{\mu }_t^{(1)}{u}^n⟩& =0,\hfill \end{array}$

(67) $\begin{array}{cc}\hfill Im⟨V{\mu }_t^{(1)}{u}^n,2{\mu }_t^{(1)}{u}^n⟩& =0,\hfill \end{array}$

(68) $\begin{array}{cc}\hfill Im⟨D{\mu }_t^{(1)}{u}^n,2{\mu }_t^{(1)}{u}^n⟩& =0,\hfill \end{array}$

(69) $\begin{array}{cc}\hfill Im⟨{\beta }_{11}|u^n{|}^2{\mu }_t^{(1)}{u}^n,2{\mu }_t^{(1)}{u}^n⟩& =0,\hfill \end{array}$

(70) $\begin{array}{cc}\hfill Im⟨{\beta }_{12}|v^n{|}^2{\mu }_t^{(1)}{u}^n,2{\mu }_t^{(1)}{u}^n⟩& =0,\hfill \end{array}$

(71) $\begin{array}{cc}\hfill Im⟨\mathrm{i}{\delta }_t^{(1)}{u}^n,2{\mu }_t^{(1)}{u}^n⟩& ={\delta }_t{\mu }_t{\parallel u^{n-1}\parallel }_2^2,\hfill \end{array}$

(72) $\begin{array}{cc}\hfill Im(⟨v^n,2{\mu }_t^{(1)}{u}^n⟩+⟨u^n,2{\mu }_t^{(1)}{v}^n⟩)& =\tau {\delta }_{t}Im(⟨u^{n-1},v^n⟩+⟨v^{n-1},u^n⟩).\hfill \end{array}$

Proof. 

Identity (a) readily follows from the fact that

(73)$Im⟨-\frac{1}{2}{\delta }_{x_i}^{({\alpha }_k)}{\mu }_t^{(1)}{u}^n,2{\mu }_t^{(1)}{u}^n⟩=-Im{\parallel {\delta }_{x_i}^{({\alpha }_k/2)}{\mu }_t^{(1)}{u}^n\parallel }_2^2,$

which obviously yields zero. On the other hand, applying distributivity of the inner product on both factors, we obtain

(74)$\begin{array}{cc}\hfill Im⟨\mathrm{i}{\delta }_t^{(1)}{u}^n,2{\mu }_t^{(1)}{u}^n⟩& =\frac{1}{2\tau }Im\left[\mathrm{i}\left({∥u^{n+1}∥}_2^2+2\mathrm{i}Im⟨u^{n-1},u^{n+1}⟩-{∥u^{n-1}∥}_2^2\right)\right]\hfill \\ \hfill & =\frac{1}{\tau }\left[\frac{1}{2}\left({∥u^{n+1}∥}_2^2+{∥u^n∥}_2^2\right)-\frac{1}{2}\left({∥u^n∥}_2^2+{∥u^{n-1}∥}_2^2\right)\right]\hfill \\ \hfill & ={\delta }_t{\mu }_t{\parallel u^{n-1}\parallel }_2^2,\hfill \end{array}$

which establishes Property (f). Finally, notice that distributivity also yields

(75)$\begin{array}{cc}\hfill Im\left(⟨v^n,2{\mu }_t^{(1)}{u}^n⟩+⟨u^n,2{\mu }_t^{(1)}{v}^n⟩\right)& =\tau {\delta }_{t}Im\left(⟨u^{n-1},v^n⟩+⟨v^{n-1},u^n⟩\right),\hfill \end{array}$

which shows the validity of (g). The proofs of Properties (b)–(e) are similar to that of (a). For that reason, we omit the proofs of the remaining identities. □

Theorem 5

(Mass conservation). If ${(u^n,v^n)}_{n=0}^N$ is a solution of (36), then ${\delta }_t{M}^n=0$ for each $n\in I_{N-2}$, where

(76)$M^n={\mu }_t\left(\parallel u^n{\parallel }_2^2+{\parallel v^n\parallel }_2^2\right)-\lambda \tau Im\left(⟨u^n,v^{n+1}⟩+⟨v^n,u^{n+1}⟩\right),\forall n\in I_{N-1},$

is the discrete total mass of the system at the time $t_n$. The discrete individual discrete masses as the quantities

(77)$\begin{array}{cc}\hfill M_1^n& ={\mu }_t{\parallel u^n\parallel }_2^2-\lambda \tau Im⟨u^n,v^{n+1}⟩,\forall n\in I_{N-1},\hfill \end{array}$

(78)$\begin{array}{cc}\hfill M_2^n& ={\mu }_t{\parallel v^n\parallel }_2^2-\lambda \tau Im⟨v^n,u^{n+1}⟩,\forall n\in I_{N-1},\hfill \end{array}$

and they are conserved if $\lambda =0$.

Proof. 

Let ${\mathsf{\Theta }}^n$ and ${\mathsf{\Phi }}^n$ be as in the proof of Theorem 4, for each $n\in I_{N-2}$. Calculating the inner product of ${\mathsf{\Theta }}^n$ and ${\mathsf{\Phi }}^n$ with $2{\mu }_t^{(1)}{u}^n$ and $2{\mu }_t^{(1)}{v}^n$, respectively, using the identities in the previous lemma. Collecting terms, we note that

(79)$\begin{array}{cc}\hfill & {\delta }_t{\mu }_t{\parallel u^{n-1}\parallel }_2^2=Im⟨\mathrm{i}{\delta }_t^{(1)}{u}^n,2{\mu }_t^{(1)}{u}^n⟩=Im⟨{\mathsf{\Theta }}^n,2{\mu }_t^{(1)}{u}^n⟩\hfill \\ \hfill & {\displaystyle =Im⟨\left(-\frac{1}{2}{△}_h^{({\alpha }_1)}+V+D+{\beta }_{11}{\left|u^n\right|}^2+{\beta }_{12}{\left|v^n\right|}^2\right){\mu }_t^{(1)}{u}^n+\lambda v^n,2{\mu }_t^{(1)}{u}^n⟩}\hfill \\ \hfill & {\displaystyle =2\lambda Im⟨v^n,{\mu }_t^{(1)}{u}^n⟩,}\hfill \end{array}$

and

(80)$\begin{array}{cc}\hfill & {\delta }_t{\mu }_t{\parallel v^{n-1}\parallel }_2^2=Im⟨\mathrm{i}{\delta }_t^{(1)}{v}^n,2{\mu }_t^{(1)}{v}^n⟩=Im⟨{\mathsf{\Phi }}^n,2{\mu }_t^{(1)}{v}^n⟩\hfill \\ \hfill & =Im⟨\left(-\frac{1}{2}{△}_h^{({\alpha }_2)}+V+{\beta }_{22}{\left|v^n\right|}^2+{\beta }_{12}{\left|u^n\right|}^2\right){\mu }_t^{(1)}{v}^n+\lambda u^n,2{\mu }_t^{(1)}{v}^n⟩\hfill \\ \hfill & =2\lambda Im⟨u^n,{\mu }_t^{(1)}{v}^n⟩,\hfill \end{array}$

for each $n\in I_{N-2}$. Observe that, indeed, the individual discrete masses are conserved if $\lambda =0$. On the other hand, adding the last two identities, we readily check that

(81)$\begin{array}{cc}\hfill {\delta }_t{M}^{n-1}& ={\delta }_t\left[{\mu }_t\left(\parallel u^{n-1}{\parallel }_2^2+{\parallel v^{n-1}\parallel }_2^2\right)-\lambda \tau Im\left(⟨u^{n-1},v^n⟩+⟨v^{n-1},u^n⟩\right)\right]\hfill \\ \hfill & ={\delta }_t{\mu }_t\left(\parallel u^{n-1}{\parallel }_2^2+{\parallel v^{n-1}\parallel }_2^2\right)-2\lambda Im⟨v^n,{\mu }_t^{(1)}{u}^n⟩-2\lambda Im⟨u^n,{\mu }_t^{(1)}{v}^n⟩\hfill \\ \hfill & =0,\hfill \end{array}$

for each $n\in I_{N-2}$, which is what we wanted to prove. □

After this last result, we refer to the quantities $M^n$ as the total mass of the discrete model (36) at the time $t_n$, for each $n\in {\overline{I}}_{N-1}$. Moreover, we define the density of mass at the point $(x_j,t_n)$ as the number

(82)$R_j^n=R(u_j^n,v_j^n)={\mu }_t(|u_j^n{|}^2+|v_j^n{|}^2)-\lambda \tau Im\left(u_j^n{\overline{v}}_j^{n+1}+v_j^n{\overline{u}}_j^{n+1}\right),$

for each $(j,n)\in J\times I_{N-1}$. It is worth noticing that the quantities $E^n$ and $M^n$ in the last theorems are clearly the discrete counterparts of the total energy and the total mass of associated to the continuous system (36), respectively.

Corollary 2

(Mass non-negativity). Let ${(u^n,v^n)}_{n=0}^N$ be a solution of (36), and suppose that $\left|\lambda \right|\tau <1$. Then, the quantities $M^n$ are non-negative and constant.

Proof. 

It only remains to prove that the quantities $M^n$ are non-negative. Using Young’s inequality, we obtain that

(83)$\begin{array}{cc}\hfill M^n& \ge {\mu }_t\left(\parallel u^n{\parallel }_2^2+{\parallel v^n\parallel }_2^2\right)-\frac{1}{2}\lambda \tau \left(\parallel u^n{\parallel }_2^2+\parallel v^{n+1}{\parallel }_2^2+\parallel v^n{\parallel }_2^2+{\parallel u^{n+1}\parallel }_2^2\right)\hfill \\ \hfill & =(1-|\lambda |\tau ){\mu }_t\left(\parallel u^n{\parallel }_2^2+{\parallel v^n\parallel }_2^2\right)\ge 0,\hfill \end{array}$

for each $n\in {\overline{I}}_{N-1}$. This obviously completes the proof. □

The proof of the following result is similar to that of Corollary 2.

Corollary 3

(Boundedness). Suppose that ${(u^n,v^n)}_{n=0}^N$ is a solution of (36), and that ${\eta }_0\in {\mathbb{R}}^{+}$ satisfies $\left|\lambda \right|\tau <{\eta }_0<1$. Then, $max\left\{\parallel u^n{\parallel }_2^2,{\parallel v^n\parallel }_2^2\right\}\le 2{(1-{\eta }_0)}^{-1}{M}^0$, for each $n\in {\overline{I}}_N$, □

4. Numerical Properties

To start with, we introduce the continuous functional

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

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

for each $(x,t)\in {\mathsf{\Omega }}_T$, and the next discrete functions, for each $(j,n)\in J\times I_{N-1}$:

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

and

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

Using this nomenclature, we agree that

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

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

for each $(x,t)\in {\mathsf{\Omega }}_T$ and $(j,n)\in J\times I_{N-1}$.

Theorem 6

(Consistency). If ${\psi }_1,{\psi }_2\in {\mathcal{C}}_{x,t}^{5,4}(\overline{{\mathsf{\Omega }}_T})$, then there exist constants $C,C^{\prime },C^{″}>0$ which are independent of h and τ, such that, for each $(j,n)\in J\times I_{N-1}$, the following hold:

(90)$\begin{array}{cc}\hfill \left|L({\psi }_1(x_j,t_n),{\psi }_2(x_j,t_n))-\mathcal{L}({\psi }_1(x_j,t_n),{\psi }_2(x_j,t_n))\right|& \le C({\tau }^2+\parallel h{\parallel }_2^2),\hfill \end{array}$

(91)$\begin{array}{cc}\hfill \left|H({\psi }_1(x_j,t_n),{\psi }_2(x_j,t_n))-\mathcal{H}({\psi }_1(x_j,t_{n+\frac{1}{2}}),{\psi }_2(x_j,t_{n+\frac{1}{2}}))\right|& \le C^{\prime }({\tau }^2+{\parallel h\parallel }_2^2),\hfill \end{array}$

(92)$\begin{array}{cc}\hfill \left|R({\psi }_1(x_j,t_n),{\psi }_2(x_j,t_n))-\mathcal{R}({\psi }_1(x_j,t_{n+\frac{1}{2}}),{\psi }_2(x_j,t_{n+\frac{1}{2}}))\right|& \le C^{″}\tau .\hfill \end{array}$

Proof. 

Use Taylor’s theorem, Lemma 2 and the conditions of regularity, there are constants $C_{1,k},C_{2,k,i},C_{3,k},C_{4,l}\in \mathbb{R}$ which are independent of h and $\tau$, for $i\in I_p$, $k\in I_2$ and $l\in I_3$, such that the next hold, for each $(j,n)\in J\times I_{N-1}$:

(93)$\begin{array}{cc}\hfill \left|{\delta }_t^{(1)}{\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}$

(94)$\begin{array}{cc}\hfill \left|{\mu }_t^{(1)}{\delta }_{x_i}^{({\alpha }_1)}{\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}$

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

(96)$\begin{array}{cc}\hfill \left|{\left|{\psi }_1(x_j,t_n)\right|}^2{\mu }_t^{(1)}{\psi }_1(x_j,t_n)-{\left|{\psi }_1(x_j,t_n)\right|}^2{\psi }_1(x_j,t_n)\right|& \le C_{4,1}{\tau }^2,\hfill \end{array}$

(97)$\begin{array}{cc}\hfill \left|{\left|{\psi }_1(x_j,t_n)\right|}^2{\mu }_t^{(1)}{\psi }_2(x_j,t_n)-{\left|{\psi }_1(x_j,t_n)\right|}^2{\psi }_2(x_j,t_n)\right|& \le C_{4,2}{\tau }^2,\hfill \end{array}$

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

The first inequality of the conclusion follows from the triangle inequality. The proof of the second inequality is similar. Finally, to prove the last inequality, observe firstly that the regularity conditions on ${\psi }_1$ and ${\psi }_2$ assure that there exist constants $C_{5,k},C_6\in \mathbb{R}$ which are independent of h and $\tau$, for $k\in I_2$, such that

(99)$\begin{array}{cc}\hfill \left|{\mu }_t|{\psi }_k(x_j,t_n){|}^2-{|{\psi }_k(x_j,t_{n+\frac{1}{2}})|}^2\right|& \le C_{5,k}{\tau }^2,\hfill \end{array}$

(100)$\begin{array}{cc}\hfill \left|\lambda \tau Im\left({\psi }_1(x_j,t_n)\overline{{\psi }_2(x_j,t_{n+1})}-{\psi }_2(x_j,t_n)\overline{{\psi }_1(x_j,t_{n+1})}\right)\right|& \le C_6\tau .\hfill \end{array}$

The third inequality of this result readily follows now assuming that $\tau <1$, using the triangle inequality, and letting $C^{″}$ be a suitable constant in terms of $C_{5,1}$, $C_{5,2}$ and $C_6$. □

Lemma 8

(Ferreira [25]). Let ${({\omega }^n)}_{n=0}^N$ and ${({\rho }^n)}_{n=0}^N$ be arrays of non-negative numbers for which there exists $C_0\ge 0$ with the property that the following hold:

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

If $n\in {\overline{I}}_N$, then ${\omega }^n\le {\rho }^n{e}^{C_{0}n\tau }$. □

To prove stability, take two solutions of our scheme (36) associated to two sets of initial data. In notation, we assume that $(u,v)$ represents a solution for (36), and that $(\tilde{u},\tilde{v})$ solves the discrete problem but with initial profiles described by functions ${\tilde{\varphi }}_1,{\tilde{\varphi }}_2:\overline{\mathsf{\Omega }}\to \mathbb{C}$. Concretely, the following discrete initial-boundary-value problem holds:

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

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

Lemma 9.

Let $ϵ={({ϵ}^n)}_{n=0}^N$ and $\zeta ={({\zeta }^n)}_{n=0}^N$ be sequences in ${\mathcal{V}}_h$. If $m\in I_{N-1}$ then

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

Proof. 

After substituting the definition of the operator ${\delta }_t$ at the left-hand side of the identity (103) and performing algebraic operations, we obtain that

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

Observe that most of the terms inside the sums cancel out. The conclusion readily follows now. □

Lemma 10.

Let $(u,v)$ and $(\tilde{u},\tilde{v})$ be solutions of (36) and (102), respectively, and let us define ${ϵ}^n=u^n-{\tilde{u}}^n$ and ${\zeta }^n=v^n-{\tilde{v}}^n$, for each $n\in {\overline{I}}_N$. For each $(j,n)\in J\times I_{N-1}$, define

(105) $\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){\mu }_t^{(1)}{u}_j^n-\left({\beta }_{11}{\left|{\tilde{u}}_j^n\right|}^2+{\beta }_{12}{\left|{\tilde{v}}_j^n\right|}^2\right){\mu }_t^{(1)}{\tilde{u}}_j^n,\hfill \end{array}$

(106) $\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){\mu }_t^{(1)}{v}_j^n-\left({\beta }_{22}{\left|{\tilde{v}}_j^n\right|}^2+{\beta }_{12}{\left|{\tilde{u}}_j^n\right|}^2\right){\mu }_t^{(1)}{\tilde{v}}_j^n.\hfill \end{array}$

There is a constant $C\ge 0$ which is independent of h and τ, such that, for each $(j,n)\in J\times I_{N-1}$,

(107) $max\left\{\right|P_j^n|,|Q_j^n\left|\right\}\le C\left(|{ϵ}_j^{n-1}|+|{ϵ}_j^n|+|{ϵ}_j^{n+1}|+|{\zeta }_j^{n-1}|+|{\zeta }_j^n|+|{\zeta }_j^{n+1}|\right).$

Proof. 

It is straightforward. □

Theorem 7

(Stability). Assume that $(u,v)$ and $(\tilde{u},\tilde{v})$ are solutions of the problems (36) and (102), respectively, and define ${ϵ}^n=u^n-{\tilde{u}}^n$ and ${\zeta }^n=v^n-{\tilde{v}}^n$, for each $n\in {\overline{I}}_N$. Assume that the condition

(108)$\left(2\lambda +4C^{\prime }+1\right)\tau <1$

is satisfied. Then, there exists a constant $C^{\prime }\ge 0$ which is independent of τ and h, such that

(109)${\mu }_t\left(\parallel {ϵ}^n{\parallel }_2^2+{\parallel {\zeta }^n\parallel }_2^2\right)\le {\mu }_t(\parallel {ϵ}^0{\parallel }_2^2+\parallel {\zeta }^0{\parallel }_2^2)e^{\left(\lambda +\frac{4C^{\prime }+1}{2}\right)T},\forall n\in {\overline{I}}_N.$

Proof. 

Obviously, the next problem is satisfied, for each $(j,n)\in J\times I_{N-1}$:

(110)$\begin{array}{c}\begin{array}{cc}\hfill {\displaystyle \mathrm{i}{\delta }_t^{(1)}{ϵ}_j^n}& =-\frac{1}{2}{△}_h^{({\alpha }_1)}{\mu }_t^{(1)}{ϵ}_j^n+V_j{\mu }_t^{(1)}{ϵ}_j^n+D{\mu }_t^{(1)}{ϵ}_j^n+P_j^n+\lambda {\zeta }_j^n,\hfill \\ \hfill {\displaystyle \mathrm{i}{\delta }_t^{(1)}{\zeta }_j^n}& =-\frac{1}{2}{△}_h^{({\alpha }_2)}{\mu }_t^{(1)}{\zeta }_j^n+V_j{\mu }_t^{(1)}{\zeta }_j^n+Q_j^n+\lambda {ϵ}_j^n,\hfill \end{array}\\ \mathrm{such}\mathrm{that}\left\{\begin{array}{cc}{ϵ}_j^0={\varphi }_1(x_j)-{\tilde{\varphi }}_1(x_j),\hfill & \forall j\in \overline{J},\hfill \\ {\zeta }_j^0={\varphi }_2(x_j)-{\tilde{\varphi }}_2(x_j),\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 $P_j^n$ and $Q_j^n$ are as in Lemma 10. Then, there is $C^{\prime }\ge 0$ independent of h and $\tau$, with the property that, for each $n\in I_{N-1}$,

(111)$\begin{array}{cc}\hfill & max\{\parallel P^n{\parallel }_2^2,\parallel Q^n{\parallel }_2^2\}\hfill \\ \hfill & \le C^{\prime }\left(\parallel {ϵ}^{n-1}{\parallel }_2^2+\parallel {ϵ}^n{\parallel }_2^2+\parallel {ϵ}^{n+1}{\parallel }_2^2+\parallel {\zeta }^{n-1}{\parallel }_2^2+\parallel {\zeta }^n{\parallel }_2^2+{\parallel {\zeta }^{n+1}\parallel }_2^2\right).\hfill \end{array}$

Obtain the inner product of the first vector Equation of (110) with $2{\mu }_t^{(1)}{ϵ}^n$ and calculate the imaginary parts. Meanwhile, in the second equation, obtain the inner product with $2{\mu }_t^{(1)}{\zeta }^n$ and calculate imaginary parts. By Theorem 5,

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

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

It is readily checked then that that there exists $C^{″}\ge 0$ independent of $\tau$ and h, such that

(114)$\begin{array}{cc}\hfill & 2{\mu }_t\left(\parallel {ϵ}^{m+1}{\parallel }_2^2+{\parallel {\zeta }^{m+1}\parallel }_2^2\right)\hfill \\ \hfill & =2{\mu }_t\left(\parallel {ϵ}^0{\parallel }_2^2+{\parallel {\zeta }^0\parallel }_2^2\right)+4\tau \sum _{n=1}^{m}Im⟨P^n,{\mu }_t^{(1)}{ϵ}^n⟩\hfill \\ \hfill & +4\tau \sum _{n=1}^{m}Im⟨Q^n,{\mu }_t^{(1)}{\zeta }^n⟩+2\tau \lambda \sum _{n=1}^m{\delta }_{t}Im\left[⟨{ϵ}^{n-1},{\zeta }^n⟩+⟨{\zeta }^{n-1},{ϵ}^n⟩\right]\hfill \\ \hfill & \le 2{\mu }_t\left(\parallel {ϵ}^0{\parallel }_2^2+{\parallel {\zeta }^0\parallel }_2^2\right)+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 & +\tau \sum _{n=1}^m\left(2\parallel P^n{\parallel }_2^2+\parallel {ϵ}^{n-1}{\parallel }_2^2+\parallel {ϵ}^n{\parallel }_2^2+2\parallel Q^{n+1}{\parallel }_2^2+\parallel {\zeta }^{n-1}{\parallel }_2^2+{\parallel {\zeta }^{n+1}\parallel }_2^2\right)\hfill \\ \hfill & \le 2\left(1+\lambda \tau \right){\mu }_t\left(\parallel {ϵ}^0{\parallel }_2^2+{\parallel {\zeta }^0\parallel }_2^2\right)\hfill \\ \hfill & +\left(4C^{\prime }+1\right)\tau \sum _{n=1}^m\left(\parallel {ϵ}^{n-1}{\parallel }_2^2+\parallel {ϵ}^n{\parallel }_2^2+\parallel {ϵ}^{n+1}{\parallel }_2^2+\parallel {\zeta }^{n-1}{\parallel }_2^2+\parallel {\zeta }^n{\parallel }_2^2+{\parallel {\zeta }^{n+1}\parallel }_2^2\right)\hfill \\ \hfill & \le 2\left(1+\left(\lambda +\frac{4C^{\prime }+1}{2}\right)\tau \right){\mu }_t\left(\parallel {ϵ}^0{\parallel }_2^2+{\parallel {\zeta }^0\parallel }_2^2\right)\hfill \\ \hfill & +4\left(4C^{\prime }+1\right)\tau \sum _{n=1}^{m-1}{\mu }_t\left(\parallel {ϵ}^n{\parallel }_2^2+{\parallel {\zeta }^n\parallel }_2^2\right)+\left(4C^{\prime }+1\right)\tau {\mu }_t\left(\parallel {ϵ}^{m+1}{\parallel }_2^2+{\parallel {\zeta }^{m+1}\parallel }_2^2\right).\hfill \end{array}$

Subtract $\left(4C^{\prime }+1\right)\tau {\mu }_t\left(\parallel {ϵ}^{m+1}{\parallel }_2^2+{\parallel {\zeta }^{m+1}\parallel }_2^2\right)$. Then, (101) is satisfied, for each $m\in I_{N-1}$, using $C_0=2\left(4C^{\prime }+1\right)$,

(115)$\begin{array}{cc}\hfill {\omega }^m& ={\mu }_t\left(\parallel {ϵ}^m{\parallel }_2^2+{\parallel {\zeta }^m\parallel }_2^2\right),\forall m\in {\overline{I}}_N,\hfill \end{array}$

(116)$\begin{array}{cc}\hfill {\rho }^m& =\left(1+\left(\lambda +\frac{4C^{\prime }+1}{2}\right)\tau \right){\mu }_t\left(\parallel {ϵ}^0{\parallel }_2^2+{\parallel {\zeta }^0\parallel }_2^2\right),\forall m\in {\overline{I}}_N.\hfill \end{array}$

The conclusion of this result is reached now using Lemma 8. □

The following result is now straightforward.

Corollary 4

(Uniqueness). For any initial data, the problem (36) has a unique solution. □

We establish next the convergence properties of our finite-difference scheme. To that end, we assume that $({\psi }_1,{\psi }_2)$ represents a pair of sufficiently smooth solutions of the initial-value problem (4). As a consequence,

(117)$\left\{\begin{array}{cc}{L}_1({\psi }_1(x_j,t_n),{\psi }_2(x_j,t_n))={\zeta }_j^n,\hfill & \forall (j,n)\in J\times {\overline{I}}_N,\hfill \\ L_2({\psi }_1(x_j,t_n),{\psi }_2(x_j,t_n))={\xi }_j^n,\hfill & \forall (j,n)\in J\times {\overline{I}}_N.\hfill \end{array}\right.$

Obviously, ${\zeta }_j^n$ and ${\xi }_j^n$ represent here the local truncation errors.

Theorem 8

(Convergence). Let ${\psi }_1,{\psi }_2\in {\mathcal{C}}_{x,t}^{5,4}({\overline{\mathsf{\Omega }}}_T)$, and let $({\psi }_1,{\psi }_2)$ solve (4). There is $C_5\ge 0$ independent of τ and h, such that the numerical solution converges to that of the continuous problem (4), whenever $2C_5\tau <1$.

Proof. 

Let $(u,v)$ be the solution of problem (36). For the sake of convenience, let us define ${\eta }_j^n={\psi }_1(x_j,t_n)-u_j^n$ and ${\theta }_j^n={\psi }_2(x_j,t_n)-v_j^n$, for each $n\in {\overline{I}}_N$. Notice that $(\eta ,\theta )$ satisfies the problem

(118)$\begin{array}{c}\begin{array}{cc}\hfill {\displaystyle \mathrm{i}{\delta }_t^{(1)}{\eta }_j^n}& =-\frac{1}{2}{△}_h^{({\alpha }_1)}{\mu }_t^{(1)}{\eta }_j^n+V_j{\mu }_t^{(1)}{\eta }_j^n+D{\mu }_t^{(1)}{\eta }_j^n+{\mathcal{P}}_j^n+\lambda {\theta }_j^n+{\zeta }_j^n,,\hfill \\ \hfill {\displaystyle \mathrm{i}{\delta }_t^{(1)}{\theta }_j^n}& =-\frac{1}{2}{△}_h^{({\alpha }_2)}{\mu }_t^{(1)}{\theta }_j^n+V_j{\mu }_t^{(1)}{\theta }_j^n+{\mathcal{Q}}_j^n+\lambda {\eta }_j^n+{\xi }_j^n,\hfill \end{array}\\ \mathrm{such}\mathrm{that}\left\{\begin{array}{cc}{\eta }_j^0={\theta }_j^0=0,\hfill & \forall j\in \overline{J},\hfill \\ {\eta }_j^n={\theta }_j^n=0,\hfill & \forall (j,n)\in \partial J\times {\overline{I}}_N,\hfill \end{array}\right.\end{array}$

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

(119)$\begin{array}{cc}\hfill {\mathcal{P}}_j^n& =\left({\beta }_{11}|{\psi }_1(x_j,t_n){|}^2+{\beta }_{12}{\left|{\psi }_2(x_j,t_n)\right|}^2\right){\mu }_t^{(1)}{\psi }_1(x_j,t_n)\hfill \\ \hfill & -\left({\beta }_{11}{\left|u_j^n\right|}^2+{\beta }_{12}{\left|v_j^n\right|}^2\right){\mu }_t^{(1)}{u}_j^n,\hfill \end{array}$

and

(120)$\begin{array}{cc}\hfill {\mathcal{Q}}_j^n& =\left({\beta }_{22}|{\psi }_2(x_j,t_n){|}^2+{\beta }_{12}{\left|{\psi }_1(x_j,t_n)\right|}^2\right){\mu }_t^{(1)}{\psi }_2(x_j,t_n)\hfill \\ \hfill & -\left({\beta }_{22}{\left|v_j^n\right|}^2+{\beta }_{12}{\left|u_j^n\right|}^2\right){\mu }_t^{(1)}{v}_j^n.\hfill \end{array}$

Using the regularity of ${\psi }_1$ and ${\psi }_2$ and Theorem 6, there exists a constant $C\ge 0$ independent of $\tau$ and h, such that $max\left\{\right|{\zeta }_j^n|,|{\xi }_j^n\left|\right\}\le C({\tau }^2+{\parallel h\parallel }_2^2)$, for each $(j,n)\in J\times {\overline{I}}_N$. As in the previous theorem, there is $C_1\ge 0$ independent of h and $\tau$, such that

(121)$max\left\{\right|{\mathcal{P}}_j^n|,|{\mathcal{Q}}_j^n\left|\right\}\le C_1\left(|{\eta }_j^{n-1}|+|{\eta }_j^n|+|{\eta }_j^{n+1}|+|{\theta }_j^{n-1}|+|{\theta }_j^n|+|{\theta }_j^{n+1}|\right),$

for each $(j,n)\in J\times I_{N-1}$. This means that there is $C_2\ge 0$, such that, for each $n\in I_{N-1}$,

(122)$\begin{array}{cc}\hfill & max\{\parallel {\mathcal{P}}^n{\parallel }_2^2,\parallel {\mathcal{Q}}^n{\parallel }_2^2\}\hfill \\ \hfill & \le C_2\left(\parallel {\eta }^{n-1}{\parallel }_2^2+\parallel {\eta }^n{\parallel }_2^2+\parallel {\eta }^{n+1}{\parallel }_2^2+\parallel {\theta }^{n-1}{\parallel }_2^2+\parallel {\theta }^n{\parallel }_2^2+{\parallel {\theta }^{n+1}\parallel }_2^2\right).\hfill \end{array}$

From here, it is easy to reach the formulas

(123)$\begin{array}{cc}\hfill {\delta }_t^{(1)}{\parallel {\eta }^n\parallel }_2^2& =2Im⟨{\mathcal{P}}^n,{\mu }_t^{(1)}{\eta }^n⟩+2\lambda Im⟨{\theta }^n,{\mu }_t^{(1)}{\eta }^n⟩+2Im⟨{\zeta }^n,{\mu }_t^{(1)}{\eta }^n⟩,\hfill \end{array}$

(124)$\begin{array}{cc}\hfill {\delta }_t^{(1)}{\parallel {\theta }^n\parallel }_2^2& =2Im⟨{\mathcal{Q}}^n,{\mu }_t^{(1)}{\theta }^n⟩+2\lambda Im⟨{\eta }^n,{\mu }_t^{(1)}{\theta }^n⟩+2Im⟨{\xi }^n,{\mu }_t^{(1)}{\theta }^n⟩,\hfill \end{array}$

for each $n\in I_{N-1}$. Moreover, there exist $C_3,C_4\ge 0$ independent of $\tau$ and h, such that

(125)$\begin{array}{cc}\hfill & \left(1-\left(\lambda +\frac{4C_2+1}{2}\right)\tau \right){\mu }_t\left(\parallel {\eta }^m{\parallel }_2^2+{\parallel {\theta }^m\parallel }_2^2\right)\hfill \\ \hfill & \le 2\left(4C_2+1\right)\tau \sum _{n=1}^{m-1}{\mu }_t\left(\parallel {\eta }^n{\parallel }_2^2+{\parallel {\theta }^n\parallel }_2^2\right)\hfill \\ \hfill & +\frac{\tau }{4}\sum _{n=1}^m\left(2\parallel {\zeta }^n{\parallel }_2^2+2\parallel {\xi }^n{\parallel }_2^2+\parallel {\eta }^{n-1}{\parallel }_2^2+\parallel {\eta }^{n+1}{\parallel }_2^2+\parallel {\theta }^{n-1}{\parallel }_2^2+{\parallel {\theta }^{n-1}\parallel }_2^2\right)\hfill \\ \hfill & \le C_3({\tau }^2+{\parallel h\parallel }_2^2{)}^2+2\left(4C_2+\frac{5}{4}\right)\tau \sum _{n=1}^{m-1}{\mu }_t\left(\parallel {\eta }^n{\parallel }_2^2+{\parallel {\theta }^n\parallel }_2^2\right)\hfill \\ \hfill & +\frac{\tau }{2}{\mu }_t\left(\parallel {\eta }^m{\parallel }_2^2+{\parallel {\theta }^m\parallel }_2^2\right).\hfill \end{array}$

Subtract $\frac{\tau }{2}{\mu }_t\left(\parallel {\eta }^m{\parallel }_2^2+{\parallel {\theta }^m\parallel }_2^2\right)$ from both ends of this inequalities, we notice now that (101) is satisfied, for each $m\in I_{N-1}$, letting $C_0=2\left(4C_2+\frac{5}{4}\right)$, ${\omega }^m={\mu }_t\left(\parallel {\eta }^n{\parallel }_2^2+{\parallel {\theta }^n\parallel }_2^2\right)$ and ${\rho }^m=C_3({\tau }^2+{\parallel h\parallel }_2^2{)}^2$, for each $m\in {\overline{I}}_N$. Lemma 8 guarantees now that

(126)${\mu }_t\left(\parallel {\eta }^n{\parallel }_2^2+{\parallel {\theta }^n\parallel }_2^2\right)\le C_5({\tau }^2+{\parallel h\parallel }_2^2{)}^2,\forall (j,n)\in J\times {\overline{I}}_N,$

where $C_5=C_3{e}^{C_{0}T}$. From this, it follows that $\parallel {\eta }^n{\parallel }_2,\parallel {\theta }^n{\parallel }_2\le \sqrt{C_5}({\tau }^2+{\parallel h\parallel }_2^2)$, for each $(j,n)\in J\times {\overline{I}}_N$, which implies that the solution of (36) converges to that of (4). □

5. Numerical Simulations

The present section is devoted to describe our computational implementation of the finite-difference scheme (36), and to provide numerical simulations that illustrate the capability of the scheme to preserve the energy and the mass of the system. Our description of our implementation and the simulations are carried out in the one-dimensional case, that is, when $p=1$. Moreover, for simplification purposes, let $a=a_1$, $b=b_1$, $h=h_1$, and $M=M_1$. Throughout this section, we represent the approximations in vector form as

(127)$w^n={(w_0,w_1,w_2,\dots ,w_M)}^{\top }\in {\mathbb{C}}^{M+1},\forall n\in {\overline{I}}_N,\forall w=u,v.$

For computational purposes, we define next various matrices, all of them with sizes equal to $(J+1)\times (J+1)$. To start with, we let ${\mathbf{A}}_u^n=\mathbf{I}+{\mathbf{D}}_u^{-}+\frac{1}{2}\tau {\mathbf{G}}^{({\alpha }_1)}-\tau {\mathbf{E}}_u^n$ and ${\mathbf{B}}_u^n={\mathbf{D}}_u^{+}-\frac{1}{2}\tau {\mathbf{G}}^{({\alpha }_1)}+\tau {\mathbf{E}}_u^n$, for each $n\in {\overline{I}}_N$, where $\mathbf{I}$, ${\mathbf{D}}_u^{\pm }$ and ${\mathbf{E}}_u^n$ are diagonal matrices defined through their diagonal entries by the formulas

(128)$\begin{array}{cc}\hfill \mathbf{I}& =diag\left(1,0,0,0,\dots ,0,1\right),\hfill \end{array}$

(129)$\begin{array}{cc}\hfill {\mathbf{D}}_u^{\pm }& =diag\left(0,\mathrm{i}\pm \tau [V(x_1)+D],\mathrm{i}\pm \tau [V(x_2)+D],\dots ,\mathrm{i}\pm \tau [V(x_{J-1})+D],0\right),\hfill \end{array}$

and

(130)$\begin{array}{cc}\hfill {\mathbf{E}}_u^n& =diag(0,{\beta }_{11}|u_1^n{|}^2+{\beta }_{12}|v_1^n{|}^2,{\beta }_{11}|u_2^n{|}^2+{\beta }_{12}{\left|v_2^n\right|}^2,\dots ,\hfill \\ \hfill & \dots ,{\beta }_{11}|u_{J-1}^n{|}^2+{\beta }_{12}{\left|v_{J-1}^n\right|}^2,0),\forall n\in {\overline{I}}_N,\hfill \end{array}$

and for each $\alpha \in (1,2]$, we define

(131)${\mathbf{G}}^{(\alpha )}=-\frac{1}{h^{\alpha }}\left(\begin{array}{ccccccc}0& 0& 0& 0& \cdots & 0& 0\\ g_{-1}^{(\alpha )}& g_0^{(\alpha )}& g_1^{(\alpha )}& g_2^{(\alpha )}& \cdots & g_{J-2}^{(\alpha )}& g_{J-1}^{(\alpha )}\\ g_{-2}^{(\alpha )}& g_{-1}^{(\alpha )}& g_0^{(\alpha )}& g_1^{(\alpha )}& \cdots & g_{J-3}^{(\alpha )}& g_{J-2}^{(\alpha )}\\ g_{-3}^{(\alpha )}& g_{-2}^{(\alpha )}& g_{-1}^{(\alpha )}& g_0^{(\alpha )}& \cdots & g_{J-4}^{(\alpha )}& g_{J-3}^{(\alpha )}\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ g_{1-J}^{(\alpha )}& g_{2-J}^{(\alpha )}& g_{3-J}^{(\alpha )}& g_{4-J}^{(\alpha )}& \cdots & g_0^{(\alpha )}& g_1^{(\alpha )}\\ 0& 0& 0& 0& \cdots & 0& 0\end{array}\right).$

For each $n\in {\overline{I}}_N$, we also introduce the matrices ${\mathbf{A}}_v^n=\mathbf{I}+{\mathbf{D}}_v^{-}+\frac{1}{2}\tau {\mathbf{G}}^{({\alpha }_2)}-\tau {\mathbf{E}}_v^n$ and ${\mathbf{B}}_v^n={\mathbf{D}}_n^{+}-\frac{1}{2}\tau {\mathbf{G}}^{({\alpha }_2)}+\tau {\mathbf{E}}_v^n$, for each $n\in {\overline{I}}_N$, where the new matrices ${\mathbf{D}}_v^{\pm }$ and ${\mathbf{E}}_v^n$ are diagonal, and they are defined as

(132)$\begin{array}{cc}\hfill {\mathbf{D}}_v^{\pm }& =diag\left(0,\mathrm{i}\pm \tau V(x_1),\mathrm{i}\pm \tau V(x_2),\mathrm{i}\pm \tau V(x_3),\dots ,\mathrm{i}\pm \tau V(x_{J-1}),0\right),\hfill \end{array}$

and

(133)$\begin{array}{cc}\hfill {\mathbf{E}}_v^n& =diag(0,{\beta }_{22}|v_1^n{|}^2+{\beta }_{12}|u_1^n{|}^2,{\beta }_{22}|v_2^n{|}^2+{\beta }_{12}{\left|u_2^n\right|}^2,\dots ,\hfill \\ \hfill & \dots ,{\beta }_{22}|v_{J-1}^n{|}^2+{\beta }_{12}{\left|u_{J-1}^n\right|}^2,0),\forall n\in {\overline{I}}_N,\hfill \end{array}$

Observe that the approximations at the times $t_0$ of the discrete model are provided exactly through the initial conditions in (36). Moreover, the approximations at the time $t_1$ are provided explicitly through Formulas (37) and (38). Meanwhile, the recursive equations of the discrete model can be alternatively expressed in terms of the nomenclature introduced in this work as

(134)$\left\{\begin{array}{cc}\hfill {\mathbf{A}}_u^n{u}^{n+1}& ={\mathbf{B}}_u^n{u}^{n-1}+2\tau \lambda v_{\ast }^n,\forall n\in {\overline{I}}_{N-1},\hfill \\ \hfill {\mathbf{A}}_v^n{v}^{n+1}& ={\mathbf{B}}_v^n{v}^{n-1}+2\tau \lambda u_{\ast }^n,\forall n\in {\overline{I}}_{N-1},\hfill \end{array}\right.$

where the $(J+1)$-dimensional vectors $u_{\ast }^n$ and $v_{\ast }^n$ are given by

(135)$\begin{array}{cc}\hfill u_{\ast }^n& ={\left(0,u_1^n,u_2^n,u_3^n,\dots ,u_{J-1}^n,0\right)}^{\top },\forall n\in {\overline{I}}_N,\hfill \end{array}$

(136)$\begin{array}{cc}\hfill v_{\ast }^n& ={\left(0,v_1^n,v_2^n,v_3^n,\dots ,b_{J-1}^n,0\right)}^{\top },\forall n\in {\overline{I}}_N.\hfill \end{array}$

It is worth pointing out that our computer implementation in Appendix A is based on the recursive system (134).

Example 1.

In the following, we use the constants in Table 1. Meanwhile, the parameters ${\alpha }_1$ and ${\alpha }_2$ change values between simulations, and $V(x)=\frac{1}{2}{x}^2$, for each $x\in \mathsf{\Omega }$. On the other hand, ${\varphi }_1(x)={\varphi }_2(x)=\frac{1}{\sqrt{\pi }}exp(-x^2)$, for each $x\in \mathsf{\Omega }$. Throughout, we fix the computer parameters $\tau =0.05$ and $h_1=0.1$. In our simulations, we consider three cases for the differentiation orders ${\alpha }_1$ and ${\alpha }_2$, namely,

• Case 1.

  ${\alpha }_1=1.1$ and ${\alpha }_2=1.9$.

• Case 2.

  ${\alpha }_1=1.5$ and ${\alpha }_2=1.5$.

• Case 3.

  ${\alpha }_1=1.9$ and ${\alpha }_2=1.1$.

Under these circumstances, Figure 2, Figure 3 and Figure 4 provide the approximate solutions for the problem (4), using the values of Cases 1–3, respectively. In each figure, the graphs provide the behavior of the following, versus $(x,t)$: (a) $Re{\psi }_1$; (b) $Re{\psi }_2$; (c) $Im{\psi }_1$; (d) $Im{\psi }_2$; (e) $|{\psi }_1|$; (f) $|{\psi }_2|$. Figure 5 provides approximations for the energy functionals for these problems. More precisely, the left column of that figure provides approximate graphs of the discrete energy density $\mathcal{H}$ versus $(x,t)$, while the graphs of the total energy $\mathcal{E}$ versus t are shown on the right column. We used Case 1 (top row), Case 2 (middle row) and Case 3 (bottom row). Obviously, the total energy of the system is approximately constant, as expected from the theoretical results. Finally, Figure 6 provides a similar analysis for the mass functionals of the mathematical model (4). The results show that the mass is approximately constant in all the cases, in agreement with the theoretical results proved in this manuscript.

Example 2.

Next, we illustrate computationally the rate of convergence of (36). To that end, we use the values in Table 1 and the other parameters in the previous examples. Take a temporal period equal to $0.5$ and set $\alpha ={\alpha }_1=1.5$. The exact solution of the problem is approximated letting $\tau =1\times {10}^{-5}$ and $h=h_1=0.005$. Define ${ϵ}_{\tau ,h}(x_j)={\psi }_1(x_j,t_N)-u_j^N$ and ${\zeta }_{\tau ,h}(x_j)={\psi }_2(x_j,t_N)-v_j^N$, for each $j\in J$. For convenience, define ${ϵ}_{\tau ,h}={({ϵ}_{\tau ,h}(x_j))}_{j\in J}$ and ${\zeta }_{\tau ,h}={({\zeta }_{\tau ,h}(x_j))}_{j\in J}$. Let ${\eta }_{\tau ,h}=({ϵ}_{\tau ,h},{\zeta }_{\tau ,h})$. Notice that

(137) $\parallel {\eta }_{\tau ,h}{\parallel }_2^2=\parallel {ϵ}_{\tau ,h}{\parallel }_2^2+{\parallel {\zeta }_{\tau ,h}\parallel }_2^2.$

Introduce the convergence rates

(138) ${\rho }_{\tau ,h}={log}_2\left(\frac{\parallel {\eta }_{2\tau ,h}{\parallel }_2}{\parallel {\eta }_{\tau ,h}{\parallel }_2}\right)$

and

(139) ${\sigma }_{\tau ,h}={log}_2\left(\frac{\parallel {\eta }_{\tau ,2h}{\parallel }_2}{\parallel {\eta }_{\tau ,h}{\parallel }_2}\right).$

Table 2 summarizes our computational investigation on the convergence of the scheme. The simulations show that (36) is quadratically convergent, as predicted by Theorem 8.

Before closing this section, it is important to point out that the extension of the Gross–Pitaevskii system to the fractional scenario (4) was proposed as a natural generalization of the classical Bose–Einstein condensates. Within that context and to the best of our knowledge, the theory of fractional operators has not found a solid justification via experimentation. However, the mathematical and numerical investigation of fractional extensions of this system is a difficult problem whose interest circumscribed to mathematical and numerical areas. Nevertheless, the authors hope that some potential applications of the present work will surface in the near future using some well-known applications of fractional derivatives (see Appendix B), and that the mathematical and numerical ideas developed in this work will lead to propose and analyze new computational methods for other physically relevant systems.

6. Conclusions

In this work, we study a finite-difference method for the coupled system of Gross–Pitaevskii equations with space-fractional derivatives in the Riesz sense. We prove the existence of solutions of our discretization for any set of initial conditions using a complex fixed-point theorem. Inspired by the existence of invariants in the continuous system, we show that some discrete forms of energy and mass functions are constant in time. As its continuous counterpart, we also prove the non-negativity of the mass function and the uniform boundedness of the approximations by imposing suitable numerical constraints. We present a rigorous numerical analysis for the method, obtaining a second order of consistency in space and time for the solutions. Moreover, we show that the discrete energy and mass densities are consistent approximations of the continuous functionals. We also provide conditions to guarantee the stability and the uniqueness of the numerical solutions, and we prove the convergence of the system. The theoretical results show that the method is quadratically convergent. It is worth pointing out that analyses of stability and convergence were carried out using a discrete version of Gronwall’s inequality. A Matlab implementation of the numerical model is provided in the Appendix of this work, and we used to produce approximations to the solutions of our continuous model. The results show that the energy and mass are approximately constant, in agreement with the theoretical results derived in this work, and improving computationally efforts already reported [20,26].

It is worth pointing out that there are already some reports by these same authors on the design of numerical methods to solve double Gross–Pitaevskii-type systems with fractional diffusion, most notably the papers by Serna-Reyes et al. [26] and Serna-Reyes and Macías-Díaz [20]. In both works, the same system of fractional partial differential equations is investigated using a finite-difference methodology, the fractional derivatives are of the Riesz type, and they are discretized using fractional-order centered differences. To start with, the proofs for the existence of solutions of the numerical model are more complicated in the case of the published papers [20,26]. Indeed, the arguments in those cases require using some fixed-point theorem in view of the nonlinear nature of those works. Meanwhile, in the present work, the existence of solutions is readily established using matrix arguments. This is due to the linear nature of the numerical solved. Moreover, the uniqueness in the present case is derived at the same time as the existence. In terms of conservation of energy, all the numerical models are capable of preserving this quantity except [26]. All the numerical models have a consistency of the second order in both space and time, and they are all stable and quadratically convergent in both space and time, for sufficiently small values of the computer parameters. Finally, in terms of implementation, the present numerical methodology is easy to implement computationally, in view that the code only requires solving a couple of linear systems at each time steps. On the other hand, the finite-difference schemes reported in [20,26] are harder to code. This is again due to the nonlinear nature of the numerical models. Since the solution functions are complex-valued, implementations of Newton’s method are useless in this point. To overcome this shortcoming, those methods were implemented using a recursive approach.

After the completion of this work, various avenues of research remain open in the investigation of fractional systems with associated energy functions. It is worth recalling that there are already various reports available in the literature on energy-conserving methods for a wide family of mathematical models and their fractional extensions [27,28]. Those fractional extensions have been mainly established for the case of Riesz space-fractional mathematical models [29]. Unfortunately, the property of energy conservation has not been studied in depth for Caputo time-fractional hyperbolic systems. There are very few reports in the literature which address this topic [30]. In fact, only some few reports for the parabolic case have shown some success in preserving the dissipation of free energy in those systems [31]. To this day, this problem remains open, and its solution is one of the avenues of research that should be tackled after the successful completion of so many on energy-conserving methods for space-fractional partial differential equations. In particular, the authors of this work are interested in developing dissipation-preserving methods for time-space-fractional extensions of the model investigated in the present manuscript, using Caputo fractional operators in time and Riesz fractional derivatives in space.

Author Contributions

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

Funding

The first author (A.J.S.-R.) would like to acknowledge the financial support of the National Council for Science and Technology of Mexico (CONACYT). The corresponding author (J.E.M.-D.) acknowledges the financial support from CONACYT through grant A1-S-45928.

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 present manuscript is a submission to the Special Issue of Mathematics MDPI on “Fractional Calculus—Theory and Applications”. The authors would like to thank the anonymous reviewers and the associated editor in charge of handling this work 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.

Appendix A. Matlab Code

The following is a Matlab implementation of the finite-difference method (36), using the vector representation (134) for the recursive equations. It is important to point out here that we chose to provide a Matlab implementation of the numerical model in view that scientists are more familiarized with this computer language rather than with C/C++ or Fortran. However, the translation of this code to those programming languages is relatively straightforward for programmers with minimal knowledge of C/C++ or Fortran.

function [X,Y,phi1,phi2,H,R,t,Energy,Mass]=gp

 

 function H=fracmatrix(M,alpha,h)

 g=zeros(1,M);

 g(1)=gamma(alpha+1)/gamma(0.5∗alpha+1)^2/h^alpha;

 for l=1:M-1

 g(l+1)=(1-(alpha+1)/(0.5∗alpha+l))∗g(l);

 end

 H=zeros(M,M);

 for l=1:M

 for m=1:M

 H(l,m)=g(abs(l-m)+1);

 end

 end

 end

 

 a1=-5;

 b1=5;

 T=10;

 

 alpha1=1.9;

 alpha2=1.1;

 beta11=1.5;

 beta12=0.5;

 beta22=1.5;

 lambda=-0.5;

 D=2;

 

 tau=0.05;

 h=0.1;

 

 x=a1:h:b1;

 t=0:tau:T;

 M=length(x);

 N=length(t);

 

 Ha1=fracmatrix(M,alpha1,h);

 Ha2=fracmatrix(M,alpha2,h);

 Ha3=fracmatrix(M,0.5∗alpha1,h);

 Ha4=fracmatrix(M,0.5∗alpha2,h);

 

 V=0.5.∗x.^2;

 

 u1=exp(-x.^2)./sqrt(pi);

 v1=exp(-x.^2)./sqrt(pi);

 

 u2=u1+1i.∗tau.∗(-0.5.∗u1∗Ha1-V.∗u1-D.∗u1-lambda.∗v1..

 -(beta11.∗abs(u1).^2+beta12.∗abs(v1).^2).∗u1);

 v2=v1+1i.∗tau.∗(-0.5.∗v1∗Ha2-V.∗v1-lambda.∗u1..

 -(beta12.∗abs(u1).^2+beta22.∗abs(v1).^2).∗v1);

 

 A1=diag(1i-tau.∗V-tau.∗D)-0.5.∗tau.∗Ha1;

 A2=diag(1i+tau.∗V+tau.∗D)+0.5.∗tau.∗Ha1;

 B1=diag(1i-tau.∗V)-0.5.∗tau.∗Ha2;

 B2=diag(1i+tau.∗V)+0.5.∗tau.∗Ha2;

 

 Energy=zeros(size(t));

 Mass=zeros(size(t));

 

 [X,Y]=meshgrid(x,t);

 phi1=zeros(size(X));

 phi2=zeros(size(X));

 H=zeros(size(X));

 R=zeros(size(X));

 phi1(1,:)=u1;

 phi1(2,:)=u2;

 phi2(1,:)=v1;

 phi2(2,:)=v2;

 

 H(1,:)=0.25.∗(abs(u1∗Ha3).^2+abs(u2∗Ha3).^2+abs(v1∗Ha4).^2+abs(v2∗Ha4).^2) ...

 +0.5.∗V.∗(abs(u1).^2+abs(u2).^2+abs(v1).^2+abs(v2).^2) ...

 +0.5.∗D.∗(abs(u1).^2+abs(u2).^2)+0.5.∗beta11.∗abs(u1).^2.∗abs(u2).^2 ...

 +0.5.∗beta22.∗abs(v1).^2.∗abs(v2).^2 ...

 +0.5.∗beta12.∗(abs(u1).^2.∗abs(v2).^2+abs(v1).⌵2.∗abs(u2).^2) ...

 +lambda.∗real(u1.∗conj(v2)+v1.∗conj(u2));

 R(1,:)=abs(u1).^2+abs(u2).^2+abs(v1).^2+abs(v2).^2 ...

 -lambda.∗tau.∗imag(u1.∗conj(v2)+v1.∗conj(u2));

 Energy(1)=h∗sum(H(1,:));

 Mass(1)=h∗sum(R(1,:));

 

 for n=3:N

 a=tau.∗(beta11.∗abs(u2).^2+beta12.∗abs(v2).^2);

 b=tau.∗(beta12.∗abs(u2).^2+beta22.∗abs(v2).^2);

 u3=linsolve(A1′-diag(a), (A2′+diag(a))∗u1′ + 2.∗lambda.∗tau.∗v2′ )′;

 v3=linsolve(B1′-diag(b), (B2′+diag(b))∗v1′ + 2.∗lambda.∗tau.∗u2′ )′;

 

 u1=u2;

 u2=u3;

 v1=v2;

 v2=v3;

 

 phi1(n,:)=u3;

 phi2(n,:)=v3;

 

 H(n-1,:)=0.25.∗(abs(u1∗Ha3).^2+abs(u2∗Ha3).^2+abs(v1∗Ha4).^2+abs(v2∗Ha4).^2) ...

 +0.5.∗V.∗(abs(u1).^2+abs(u2).^2+abs(v1).^2+abs(v2).^2) ...

 +0.5.∗D.∗(abs(u1).^2+abs(u2).^2)+0.5.∗beta11.∗abs(u1).^2.∗abs(u2).^2 ...

 +0.5.∗beta22.∗abs(v1).^2.∗abs(v2).^2 ...

 +0.5.∗beta12.∗(abs(u1).^2.∗abs(v2).^2+abs(v1).^2.∗abs(u2).^2) ...

 +lambda.∗real(u1.∗conj(v2)+v1.∗conj(u2));

 R(n-1,:)=abs(u1).^2+abs(u2).^2+abs(v1).^2+abs(v2).^2 ...

 -lambda.∗tau.∗imag(u1.∗conj(v2)+v1.∗conj(u2));

 Energy(n-1)=h∗sum(H(n-1,:));

 Mass(n-1)=h∗sum(R(1,:));

 end

 Energy(N)=Energy(N-1);

 Mass(N)=Mass(N-1);

 H(N,:)=H(N-1,:);

 R(N,:)=R(N-1,:);

end

Appendix B. Fractional Calculus

It is important to recall that Riesz space-fractional operators have been derived mathematically from physical systems with long-range interactions. Indeed, it has been established thoroughly that (parabolic or hyperbolic) partial differential equations with Riesz diffusion are obtained as the continuum limit of some particle systems with long-range interactions. This continuum-limit process involves the application of several operators, like Fourier series transform, the limit when the distance between consecutive particles goes to zero along with the inverse Fourier transform [32,33].

More precisely, suppose that s is a real number in the set $\{1,2\}$ and allow $\alpha >0$. We restrict our attention to the one-dimensional case and consider an array of particles on a straight line whose motion equations are described by the coupled system of ordinary differential equations

(A1)

$\frac{d^s{u}_n}{d{t}^s}(t)=I_n(u(t))+F(u_n(t)),\forall t\in {\mathbb{R}}^{+},\forall n\in \mathbb{Z}.$

In this context, $u_n$ is a function that stands for the displacement of the nth particle with respect to the equilibrium point and F denotes the nonlinear interaction of the oscillators with some external force: Meanwhile, we use h to denote the distance existing between two consecutive oscillators. Moreover, a general form for the interaction function $I_n$ is fixed as

(A2)

$I_n(u(t))=\sum _{\begin{array}{c}m=-\infty \\ m\ne n\end{array}}^{\infty }J(n,m)\left[u_n(t)-u_m(t)\right],\forall t\in {\mathbb{R}}^{+},\forall n\in \mathbb{Z},$

where $J\in {\mathcal{L}}^2(\mathbb{Z})$ is such that $J(n,m)=J(n-m)=J(m-n)$ for each pair $m,n\in \mathbb{Z}$. For the sake of convenience, we introduce the function

(A3)

$J_{\alpha }(k)=\sum _{\begin{array}{c}n=-\infty \\ n\ne 0\end{array}}^{\infty }{e}^{-ikn}J(n),\forall k\in \mathbb{R},$

and assume that the following is satisfied:

(A4)

$A_{\alpha }=\underset{k\to 0}{lim}\frac{J_{\alpha }(k)-J_{\alpha }(0)}{{\left|k\right|}^{\alpha }}\in \mathbb{R}\backslash \left\{0\right\}.$

On both ends of (A1), apply the Fourier series transform. Next, take the limit $h\to 0$ and obtain finally the inverse Fourier transform. In such way, the following fractional partial differential equation is reached:

(A5)

$\frac{{\partial }^{s}u}{\partial t^s}(x,t)-h^{\alpha }{A}_{\alpha }\frac{{\partial }^{\alpha }u}{{\partial \left|x\right|}^{\alpha }}(x,t)-F(u(x,t))=0,\forall (x,t)\in \mathbb{R}\times {\mathbb{R}}^{+}.$

Here, the fractional derivative of order $\alpha$ is understood in the sense of Riesz [32]. It is worth pointing out here that more details on these derivations are provided in [34].

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

Figures and Tables

Figure 1. Stencil of the numerical model (36) at the time tn in one spatial dimension. The black circles denote the known solutions at the times tn−1 and tn. Meanwhile, the crosses represent the unknown approximations at the time tn+1. For the sake of convenience, we agree that xji=x1,ji.

Figure 2. Approximate solutions for the problem (4) using the parameters in Table 1, and the computational time-steps h=0.1 and τ=0.05 for the finite-difference method (36). We used V(x)=12x2, and initial data ϕ1(x)=ϕ2(x)=1πexp(−x2). The graphs provide the behavior of (a) Reψ1, (b) Reψ2, (c) Imψ1, (d) Imψ2, (e) |ψ1|, and (f) |ψ2|, versus (x,t). In these simulations, we used α1=1.1 and α2=1.9.

Figure 3. Approximate solutions for the problem (4) using the parameters in Table 1, and the computational time-steps h=0.1 and τ=0.05 for the finite-difference method (36). We used V(x)=12x2, and initial data ϕ1(x)=ϕ2(x)=1πexp(−x2). The graphs provide the behavior of (a) Reψ1, (b) Reψ2, (c) Imψ1, (d) Imψ2, (e) |ψ1|, and (f) |ψ2|, versus (x,t). In these simulations, we used α1=α2=1.5.

Figure 4. Approximate solutions for the problem (4) using the parameters in Table 1, and the computational time-steps h=0.1 and τ=0.05 for the finite-difference method (36). We used V(x)=12x2, and initial data ϕ1(x)=ϕ2(x)=1πexp(−x2). The graphs provide the behavior of (a) Reψ1, (b) Reψ2, (c) Imψ1, (d) Imψ2, (e) |ψ1|, and (f) |ψ2|, versus (x,t). In these simulations, we used α1=1.9 and α2=1.1.

Figure 5. Approximate solutions for the problem (4) using the parameters in Table 1, and the computational time-steps h=0.1 and τ=0.05 for the finite-difference method (36). We used V(x)=12x2, and initial data ϕ1(x)=ϕ2(x)=1πexp(−x2). (left column) Approximate graphs of the discrete energy density H versus (x,t). (right column) Graphs of the total energy E versus t. We used α1=1.1 and α2=1.9 (a,b), α1=1.5 and α2=1.5 (c,d), and α1=1.9 and α2=1.1 (e,f).

Figure 6. Approximate solutions for the problem (4) using the parameters in Table 1, and the computational time-steps h=0.1 and τ=0.05 for the finite-difference method (36). We used V(x)=12x2, and initial data ϕ1(x)=ϕ2(x)=1πexp(−x2). (left column) Approximate graphs of the discrete mass density R versus (x,t). (right column) Graphs of the total mass M versus t. We used α1=1.1 and α2=1.9 (a,b), α1=1.5 and α2=1.5 (c,d), and α1=1.9 and α2=1.1 (e,f).

Table 1

Table of the parameter values employed in the numerical simulations of this work.

p	${\mathit{a}}_1$	${\mathit{b}}_1$	T	${\mathit{\beta }}_{11}$	${\mathit{\beta }}_{12}$	${\mathit{\beta }}_{22}$	$\mathit{\lambda }$	D
1	$-5$	5	10	$1.5$	$0.5$	$1.5$	$0.5$	2

Table 2

Table of absolute errors in the Euclidean norm and (a) temporal and (b) spatial rates of convergence for various values of the parameters $\tau$ and h. The experiment corresponds to that described in Example 2.

(a) Temporal study of convergence
	$0.04$		$h=0.02$		$h=0.01$
$\tau$	$\parallel {ϵ}_{t,h}{\parallel }_2$	${\rho }_{\tau ,h}$	$\parallel {ϵ}_{t,h}{\parallel }_2$	${\rho }_{\tau ,h}$	$\parallel {ϵ}_{t,h}{\parallel }_2$	${\rho }_{\tau ,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) Spatial study of convergence
	$\tau =0.02$		$\tau =0.01$		$\tau =0.005$
h	$\parallel {ϵ}_{t,h}{\parallel }_2$	${\sigma }_{\tau ,h}$	$\parallel {ϵ}_{t,h}{\parallel }_2$	${\sigma }_{\tau ,h}$	$\parallel {ϵ}_{t,h}{\parallel }_2$	${\sigma }_{\tau ,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$