High-order schemes for the fractional coupled nonlinear Schrödinger equation

1.   Introduction

In 2000, Laskin established the nonlinear Schrödinger (NLS) model ^[18] with fractional Laplacian operator by extending the Feynmann path integral to the Lévy path ^[17]. The derived equation is more accurate than the traditional NLS equation in describing the variation law of quantum states of nonlocal physical systems with time ^[16]. In recent, scholars have done a lot of research on the fractional NLS equation in theoretical analysis and numerical algorithms, and achieved fruitful results ^[11,13,34].

In this paper, we aim to present and analyze high-order schemes with conservation properties for the fractional CNLS equation with periodic boundary conditions that can be written in the form ^[28]

where $\text{i}$ is the imaginary unit root, $\Omega \subset \mathbb{R}^d (d = 1, 2)$, $\vartheta$, $\zeta$ and $\varepsilon$ are real constants, $\varphi(x, t)$ and $\phi(x, t)$ are complex valued functions, $1 < s\leq2$, and $\mathcal{L}^{\frac{s}{2}}$ denotes the fractional Laplacian operator and can be expressed by the Fourier transform, namely

where $\widehat{\varphi}(\xi) = \int_{\Omega}\varphi(x)\text{e}^{-\text{i}\xi x}d x$ represents the Fourier transform for $\varphi(x)$. The fractional CNLS equation will reduce to the classical CNLS equation when $s = 2$. Therefore, similar to the classical system, the fractional CNLS equation with periodic boundary conditions possesses the following physical invariants that do not change with time ^[30]

where $\mathcal{M}$ is the fractional mass, and we refer to $\mathcal{H}$ as the fractional Hamiltonian energy.

By setting $\varphi = \widetilde{p}+\text{i}\widetilde{q}$, $\phi = \widehat{p}+\text{i}\widehat{q}$, the original system (1.1), (1.2) can be rewritten as the following real system

According to the variational derivative formula ^[33], the original fractional CNLS system (1.6)–(1.9) can be expressed by an infinite-dimensional Hamiltonian system

where $\frac{\delta H(y)}{\delta y}$ is the vector of variational derivatives for $y$, $y = \big(\widetilde{q}, \widehat{q}, \widetilde{p}, \widehat{p}\big)^{T}$, $I_2$ is the second order unit matrix, and

In recent years, some numerical algorithms ^{[15,23,24,37]} have been given for solving the fractional NLS equation as explicit analytical solutions of fractional differential equations can not be obtained ^[10]. However, these traditional algorithms can not inherit the conservative laws of the equation and can not be implemented in long-time numerical calculations. Inspired by the idea of the structure-preserving algorithm of the classical Hamiltonian system, scholars have developed many conservative schemes for the fractional CNLS system. For example, Wang and Xiao ^[29,30,31] first derived invariants of system (1.1) and (1.2) and constructed finite difference schemes to conserve these invariants; the existence and uniqueness of solutions and convergence of schemes are discussed in these works. In ^[22,26], Li et al. developed an efficient conservative difference and finite element scheme for solving the strongly fractional CNLS equation. There are many corresponding structure algorithms for fractional NLS-type equations, and the readers can refer to the literature ^[6,32,35].

Unfortunately, these numerical schemes mentioned above only have second-order accuracy temporally or require small time step sizes to obtain satisfactory numerical solutions for NLS-type equations. Therefore, it is necessary to construct and analyze high-accuracy conservative numerical schemes for the equation. Over the past decade, some numerical methods have been proposed to construct high-accuracy conservative methods for Hamiltonian systems ^[1,21]. These methods can also be applied to develop high order conservative schemes for fractional NLS-type equations. But, these schemes can both preserve not the mass and energy. Fortunately, some energy quadratic methods exist, including the invariant energy quadratization (IEQ) ^[9,36] and the scalar auxiliary variable (SAV) approaches ^[3,25,27], originally developed for gradient flow models. These methods have been applied to construct high-accuracy modified energy-preserving schemes for conservative equations ^{[7,8,14,19,20]} by combining with the symplectic Runge-Kutta method ^[12]. The quadratic auxiliary variable (QAV) approach is newly proposed to construct high-accuracy schemes for the KdV equation ^[2], the derived schemes can conserve the original energy of the system. This method provides us with a way to construct some numerical schemes that can preserve the original energy of conservative systems. But, the QAV method is not applicable for some equations with complex energy functionals.

As is shown in Eq (1.5), the Hamiltonian energy functional is very complicated and has high degree polynomials, and cannot be handled by introducing a single auxiliary variable. Therefore, the original QAV approach can not be used to construct numerical schemes to inherit the Hamiltonian energy for the equation. Influenced by the idea of the QAV method, we developed a new method which we called the multiple quadratic auxiliary variables (MQAV) approach by introducing more quadratic auxiliary variables and construct conservative schemes for system (1.1), (1.2) based on the symplectic Runge-Kutta method. The constructed schemes can preserve the mass and Hamiltonian energy of the system and have high-accuracy in the temporal direction. In addition, the MQAV approach also can be applied to other conservative systems with complex energy functionals and high degree polynomials.

The following is the outline of this paper. Section 2 gives an equivalent system for the fractional CNLS equation by introducing two quadratic auxiliary variables. Section 3 applies the symplectic Runge-Kutta method to the equivalent system in the time direction can derive a conservative semi-discrete scheme, and discusses the scheme's conservation properties. Section 4 uses the Fourier pseudo-spectral method to approximate the semi-discrete system to obtain a fully-discrete conservative scheme and gives the proposed schemes' stability. In Section 5, numerical results are displayed to illustrate the theoretical results. The final section contains some conclusions.

2.   Multiple quadratic auxiliary variable system

We first derive an equivalent system for the fractional CNLS system based on the idea of the MQAV method.

By setting $u = |\varphi|^{2}$, $v = |\phi|^{2}$, the original Hamiltonian energy can be rewritten as

Considering the following initial conditions

and taking the variational derivative for the new energy Eq (2.1), the original fractional CNLS system (1.1), (1.2) can be rewritten as the following equivalent MQAV system

here $(\ast, \bullet): = \int_{\Omega}\ast \bar{\bullet} dx$, $\bar{\bullet}$ is the conjugate of $\bullet$, $\text{Re}(\ast)$ represents taking the real part of $\ast$.

Theorem 2.1. The newly proposed MQAV system can inherit some invariants, namely

where

and $\mathcal{M}(t)$, $\mathcal{E}(t)$ are given in Eqs (1.4) and (2.1), respectively.

Proof. First, system (2.4), (2.5) can be rewritten as

and adding the above two formulas we can derive

Then, we compute the inner products of equation (2.2), (2.3) with $\varphi$ and $\phi$, and take the imaginary part of the derivation, leading to

Further, by computing the inner products of system (2.2)–(2.5) with $\varphi_t, \phi_t, u, v$, respectively, we derive

Remark 2.1. It is worth noting that the equivalent MQAV system not only preserves the modified energy and mass but also conserves a new invariant $\mathcal{I}(t)$, which implies the modified energy $\mathcal{E}(t)$ is equivalent to the original Hamiltonian energy $\mathcal{H}(t)$ in the continuous case. For the constructed new system, the auxiliary variables relation can be regarded as weak properties, while the original energy reduces to a weak invariant ^[12].

3.   Conservative semi-discrete MQAV-RK scheme

We all know that the construction of the high-accuracy conservative schemes for the fractional CNLS equation is challenging. Fortunately, the MQAV system (2.2)–(2.5) can preserve the quadratic conserved quantities. This allows us to use the symplectic RK method that can conserve the system's quadratic invariants to construct a high-order accuracy conservative scheme for the fractional CNLS equation. First, we outline some notations of the symplectic RK method.

3.1.   Symplectic RK methods

We consider the following ordinary differential equation

We define $\tau = \frac{T}{N}$ as the time step, and set $t_n = n\tau$, $n = 0, 1, \cdots, N$. An $s$-stage RK method for system (3.1) in $[t_n, t_{n+1}]$ is given by

where $a_{ij}, b_i$, $c_i = \sum\limits_{i = 1}^{s}a_{ij}$ $(i, j = 1, 2, \cdots, s)$ are real coefficients. If these coefficients have the following relationship

the proposed RK method is called as the symplectic RK method ^[12].

3.2.   Semi-discrete high-order conservative scheme

To develop a high-order conservative scheme for the MQAV system that preserves some quadratic invariants in the time direction, a natural idea is to use the symplectic RK method to approximate system (2.2)–(2.5). The resulting semi-discrete RK scheme is given as follows

Then $(\varphi^{n+1}, \phi^{n+1}, u^{n+1}, v^{n+1})^T$ can be updated by

where $z^n$ represents the numerical approximation to $z(\cdot, t)$ and $t_n$ for any function $z$. For the semi-discrete scheme, we can derive the following theorem for the structure-preserving properties.

Theorem 3.1. The system (3.5), (3.6) can inherit the following invariants, i.e.,

where ${\mathcal{M}}^n$ is the mass function and is defined by

and

$\mathcal{E}^n$ is the energy function, which has the form

Proof. According to system (3.6) and the definition of $\mathcal{M}^n$, we have

Plugging $k_{\varphi}^i = -\text{i}\Big(\vartheta\mathcal{L}^{\frac{s}{2}} \varphi_{in}-\zeta\big(u_{in}+\varepsilon v_{in}\big) \varphi_{in}\Big)$, $k_{\phi}^i = -\text{i}\Big(\vartheta\mathcal{L}^{\frac{s}{2}} \phi_{in}-\zeta\big(v_{in}+\varepsilon u_{in}\big) \phi_{in}\Big)$ into above system, we derive

where symplectic condition (3.4) and $\sum\limits_{i, j = 1}^s a_{ij}b_ik_ik_j$ = $\sum\limits_{i, j = 1}^s a_{ji}b_jk_ik_j$ were used. Noticing that

and

This together with Eq (3.12) can derive

Then, we further obtain

where the symplectic condition (3.4) and $\sum\limits_{i, j = 1}^s a_{ij}b_ik_ik_j$ = $\sum\limits_{i, j = 1}^s a_{ji}b_jk_ik_j$ were used. Furthermore, similar discussions, we obtain

Therefore, we derive

Based on the definition of $\mathcal{E}^{n}$, we can deduce

The proof is completed.

Then, under the consistent initial condition $u^0 = |\varphi^0|$, $v^0 = |\phi^0|$, we have the following original energy conservation theorem for the semi-discrete MQAV-RK system.

Theorem 3.2. The MQAV-RK system (3.5), (3.6) can also inherit the original semi-discrete Hamiltonian energy of the fractional CNLS system, namely

where ${\mathcal{H}}^n$ is the original semi-discrete energy function and given by

Proof. Based on Theorem 3.1 and the two introduced quadratic auxiliary variable $u = |\varphi|^2$, $v = |\phi|^2$, one can obtain $u^n = |\varphi^n|^2$, $v^n = |\phi^n|^2$ for $0\leq n \leq N$, which leads to

and we can immediately obtain $\mathcal{H}^n = \mathcal{H}^0$.

4.   Fully-discrete conservative scheme

4.1.   Fourier pseudo-spectral method

The proposed conservative scheme has fourth-order accuracy in time, and the Fourier pseudo-spectral method is chosen to approximate the semi-discrete system in space so that the accuracy of time direction is not affected. We consider system $(1.1)$, $(1.2)$ in a periodic region $\Omega = [-\mathbf{L}, \mathbf{L}]$, and define $x_j = -\mathbf{L}+jh = -\mathbf{L}+j\frac{2\mathbf{L}}{\mathbf{N}}, 0\leq j\leq \mathbf{N}-1$, $\mathbf{N}$ is a positive even integer, and $h: = \frac{2\mathbf{L}}{\mathbf{N}}$. Furthermore, we denote

be a vector space of grid functions with ${v}_{j} = v(x_j)$. Then defining the interpolation space $S^{p}_N: = \text {span}\{g_{j}(x) \}$ with

where $\mu = \frac{\pi}{\mathbf{L}}$, $a_{l} = \left\{\begin{aligned} & 1, \; \; |l| < N/2, \\ & 2, \; \; |l| = N/2, \end{aligned}\right.$. Then, we set $I_{N}:L^{2}(\Omega)\rightarrow S^{p}_{N}$ as the interpolation operator

with

where $\widehat{v}_{-\mathbf{N}/2} = \widehat{v}_{\mathbf{N}/2}$. Applying the Laplacian operator $\mathcal{L}^{\frac{s}{2}}$ to $I_{\mathbf{N}}v(x)$ leads to

where ${{{{{\boldsymbol{v}}}}}}\in \mathcal{V}_h$, and

and the discrete Fourier transform and its inverse are given by

It is can be computed by using the FFT technique in long time simulations.

4.2.   Conservative Fourier pseudo-spectral scheme

This subsection aims to give a high-order fully-discrete scheme for the Eqs (1.1), (1.2), and the scheme is derived by using the Fourier pseudo-spectral method to approximate equation (3.5), (3.6) in space, namely

Then $({\boldsymbol \varphi}^{n+1}, {\boldsymbol \phi}^{n+1}, {{{\boldsymbol{u}}}}^{n+1}, {{{\boldsymbol{v}}}}^{n+1})^T$ can be updated by

For a fully-discrete scheme, we can derive the following theorem for the structure-preserving properties.

Theorem 4.1. System (4.4), (4.5) can conserve some conservation laws, i.e.,

where ${{M}}^n$ is the mass function and defined by

and

and ${E}^n$ is the energy function, which has the form

where $\langle\star, \bullet\rangle = h\sum\limits_{j = 0}^{\mathbf{N}-1} \star_j\bar{\bullet}_j$.

Proof. The proof process is the same as the Theorem 3.1.

In addition, we can also prove the developed schemes can inherit the original Hamiltonian energy.

Theorem 4.2. System (4.4), (4.5) can conserve the original Hamiltonian energy of the fractional CNLS system, namely

where

Remark 4.1. Actually, the QAV method is similar to the IEQ and SAV methods, both of which are energy quadratic methods. Therefore, they can also be used to develop high-order conservative schemes for conservative equations. However, the derived schemes based on the IEQ and SAV methods can only inherit the modified energy ^[5,8]. Theorem 4.1 and Theorem 4.2 indicate that the constructed high-order MQAV-RK scheme can inherit both the modified energy and the Hamiltonian energy.

According to the mass conservation law of Theorem 4.1, we can conclude that the proposed schemes are stable.

5.   Numerical experiments

We give some examples to confirm the conservation properties and the high-accuracy of the developed schemes. The exact solution of the system is not given, and thus we compute the numerical errors by using the formula

The accuracy of the constructed schemes can be computed by

where $\tau_{j}$, $\text {error}_{j}, ({j} = 1, 2)$ are the time step and the maximum-norm errors with $\tau_{j}$, respectively. Furthermore, to test the conservation properties, we define the relative errors of some invariants as

where $M^{n}$, $H^{n}$ and $E^{n}$ denote the invariant at $t = n\tau$. Without losing generality, the paper will consider the following $2s$-4th order symplectic RK method ^[12]

To compare the advantages of the constructed scheme with some conservative schemes, we define

● MQAV-RK4: The constructed fourth-order conservative MQAV-RK scheme in this paper;

● SAV-RK4: The proposed fourth-order conservative RK scheme based on the SAV method ^[8].

● CN2: The proposed second-order Crank-Nicolson conservative scheme for the Eqs (1.1) and (1.2) ^[30].

Example 5.1. First, consider a decoupled fractional NLS equation in $\Omega = (-\pi, \pi)$ with the form

with the exact solution

In the simulation, we take $\eta = 2$, $\gamma = 1$, $\theta = 2$ to exhibit the theoretical results.

First, we set $s = 1.6$ and $T = 1$ to test the convergence of the given three numerical schemes. The maximum-norm errors and convergence rates are shown in Table 2. The listed data in the table indicate that the MQAV-RK4 and SAV-RK4 both have fourth-order convergence rates in time, and the CN2 scheme only gets second-order accuracy. Furthermore, we investigate the maximum-norm errors of the numerical solution versus the CPU time using the three selected schemes. The results are summarized in Figure 1. For a given global error, we observe that (ⅰ) the cost of the CN2 scheme is the most expensive because of the low-order accuracy in time; (ⅱ) the cost of SAV-RK4 is much cheaper than the CN2 scheme; (ⅲ) the cost of our proposed MQAV-RK4 scheme is the cheapest. We test the conservative properties of four conservative schemes at $T = 100$ by using the denoted formulas in (5.3). Figure 2 reveals relative errors of conservation laws. As is shown above, the SAV-RK4 scheme can only conserve the modified energy and mass of the equation. But, our proposed methods can both inherit the original energy and mass.

Example 5.2. This example considers the fractional CNLS equation with initial conditions ^[30]

In this example, $x\in(-20, 20)$, and we choose $\delta_{0} = 5, \omega_{0} = 3, \vartheta = 1, \zeta = 2$, $\varepsilon = 1$.

First, Eqs (5.1) and (5.2) are used to compute the $L^{\infty}-$normal and corresponding convergence rate of time for some conservative schemes at $T = 1$. The obtained numerical results are given in Table 3. From Table 3, we can find the MQAV-RK4 scheme has fourth-order convergence. By using the denoted formulas in (5.3), we test the conservative properties of four conservative schemes at $T = 50$. Figure 3 reveals relative errors of conservation laws. As is shown above, the proposed method can both inherit the mass and original energy exactly. Finally, we define $\Psi = |\boldsymbol \varphi^n|+|\boldsymbol \phi^n|$ and study the relationship between the evolution of numerical solution and the fractional $s$. As is shown in Figure 4, the $s$ affects the evolution of solitons, and the shape of the solitons change quickly as $s$ changes.

Example 5.3. In this example, we choose the initial condition

where $(x, y)\in [-8, 8]^2$, and, some parameters are given as follows

In Table 4, we give the accuracy of the proposed scheme at $s = 1.3, 1.7, 2$. As is shown in the table, the scheme has fourth-order accuracy in the time direction for different $s$. In addition, the numerical errors also increase with the increase of $s$. Figure 5 implies that the scheme can inherit the original energy and mass of the system for different $s$.

6.   Conclusions

In this paper, we first propose the Hamiltonian structure of the fractional CNLS equation, and then develop the MQAV approach for the fractional CNLS equation with a complex energy function that can not be handled by using the original QAV approach. Then, a class of MQAV-RK schemes is proposed by combining with the symplectic Runge-Kutta method. The constructed schemes have high-accuracy and inherit the mass and Hamiltonian energy. Numerical results also demonstrate that the schemes have good numerical stability. In addition, the approaches presented in this work can be used to solve other conservative differential equations.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The research is supported by the Science Foundation of Henan (No. 212300410323).

Conflict of interest

The authors declare there is no conflict of interest.