A linearly implicit energy-preserving exponential time differencing scheme for the fractional nonlinear Schrödinger equation

1.   Introduction

In 2000, Laskin, a Canadian scholar, explored non-Brownian motion as a possible approach to deriving the path integral. He expanded on Feynmann's path integral by incorporating Levy paths, as documented in his works ^[11,12], which allowed for the modification of the Schrödinger equation. By utilizing the non-local Laplacian operator $(-\Delta)^{\alpha/2} (1 < \alpha\leq2)$, he established the Schrödinger model ^[13], which accurately describes the quantum state variations of non-local physical systems over time, surpassing the classical Schrödinger equation. The model has been successfully applied to various fields ^[9,15,23]. This paper focuses on the fractional nonlinear Schrödinger (NLS) equation with the form

where $1 < \alpha \leq 2$, $\text{i}^2 = -1$, $u_0({\mathit{\boldsymbol{x}}})$ is a given smooth function, $\beta$ is a real positive constant. The fractional Laplacian $-(-\Delta)^{\frac{\alpha}{2}}$ with period $2{L}$ can be defined by ^[22]

where $\upsilon_{k_1} = \frac{k_1\pi}{L}$, $\upsilon_{k_2} = \frac{k_2\pi}{L}$, $\boldsymbol{\upsilon}_{\textbf{k}}\cdot({\mathit{\boldsymbol{x}}}+L) = \upsilon_{k_1}(x+L)+\upsilon_{k_2}(y+L)$. By setting $u = p+\text{i}q$, we can rewrite system (1.1) as a pair of real-valued equations

Guo et al. derive the equation possessing the energy conservation law ^[6], namely

where $\mathcal{H}$ is energy function. From variational derivative formula ^[4], systems (1.4) and (1.5) can be rewritten as as an Hamiltonian system

It is widely acknowledged that conserving energy is crucial in establishing the existence and uniqueness of solutions for PDEs ^[3,19]. Traditional methods fail to maintain these conservation laws, making structure-preserving algorithms essential for the numerical analysis of classical differential equations. Recently, researchers have been exploring structure-preserving numerical schemes for fractional equations ^[18], with a focus on the fractional NLS equation. For example, Li et al. ^[16] developed a fast conservative finite element scheme for the strongly coupled fractional NLS equation. Wang et al. ^[20,21] also presented different schemes that can exactly conserve energy for various fractional NLS equations. Additional related research can be found in ^[4] and the references therein. However, most structure-preserving schemes for the fractional NLS equation are fully implicit, requiring iterations to solve algebraic systems at each time step. This leads to a substantial number of calculations.

In ^[2], the exponential integrators scheme was studied in depth. This scheme has a unique advantage in that it accurately evaluates the contribution of the linear part, allowing for stable and precise computations even when dealing with highly rigid linear terms. As a result, this scheme has been successfully applied to Hamiltonian PDEs ^[17] and has been found to permit larger step sizes and achieve greater accuracy than non-exponential schemes. Additionally, the invariant energy quadratization (IEQ) method, developed by Yang and colleagues, has been used to create efficient and accurate linearly-implicit schemes for certain gradient models ^[5] and conservative differential equations ^[1,10]. For example, in ^[1], a linearly implicit scheme for the fractional NLS equation was developed using the IEQ method. However, practical computations using these schemes still require smaller step sizes, making them computationally expensive for long-time numerical simulations. To address these challenges, we propose a novel approach that combines the exponential time differencing method with the IEQ approach to develop an efficient linearly-implicit energy-preserving scheme for the fractional NLS equation. Surprisingly, there has been little research in this area to date. The main goal of this paper is to develop a linearly-implicit exponential integrator scheme for this equation, using the IEQ method. Our contributions can be summarized as follows: (i) we introduce the Hamiltonian form of the equation; (ii) we combine the IEQ method and the exponential integrators technique to develop a linearly implicit scheme that conserves energy; (iii) we demonstrate that the construction process of the proposed scheme can be extended to develop conservative schemes for other differential equations. Moreover, our numerical results show that the scheme inherits the modified energy and has better numerical stability than traditional energy-preserving schemes in practical computations.

The paper is arranged as follows. In Section 2, we obtain a new equivalent system with modified energy via the IEQ method. In Section 3, we construct a linearly-implicit exponential integrator scheme for solving the fractional NLS equation, and prove that the proposed scheme can preserve the discrete energy. A numerical example is presented in Section 4 to illustrate the theoretical results. We draw some conclusions in Section 5.

2.   IEQ method for the equation

In this Section, we utilize the IEQ idea to transform Eqs (1.4) and (1.5) into an equivalent system. Following the idea of the IEQ approach, we introduce an auxiliary variable $w: = w(t) = \sqrt{G(p, q)+C_0}$, where $G(p, q) = \frac{\beta}{4}(p^2+q^2)^2$, by choosing a sufficiently large constant $C_0$, the well-posedness of $w$ is ensured. This results in a modified formalism for the Hamiltonian energy (1.6).

and the corresponding fractional NLS systems (1.4) and (1.5) can be written as

with the consistent initial condition $p_0 = p_0({\mathit{\boldsymbol{x}}}, 0) = \text{Re}\big(u_0({\mathit{\boldsymbol{x}}})\big), q_0 = q_0({\mathit{\boldsymbol{x}}}, 0) = \text{Im}\big(u_0({\mathit{\boldsymbol{x}}})\big), w_0 = \sqrt{G(p_0, q_0)+C_0},$ and the periodic boundary condition. It is easy to prove the systems (2.2)–(2.4) possesses the modified energy

3.   Construction of conservative scheme

3.1.   Spatial discretization

For positive even integers $N_x$, $N_y$ and $N$, we define: $h_x: = \frac{2L}{N_x}, h_y: = \frac{2L}{N_y}$, and set time step $\tau = \frac{T}{{N}}$. The spatial grid points are given: $\Omega_{h} = \{(x_i, y_j) |\ i = 0, 1, \cdots, N_x-1; j = 0, 1, \cdots, N_y-1\}$, where $x_i = -L+ih_x, 0\leq i \leq N_x-1, y_j = -L+jh_y, 0\leq j \leq N_y-1$. Let $\mathcal{U}_h = \{{\mathit{\boldsymbol{u}}}|{\mathit{\boldsymbol{u}}} = (u_{0, 0}, \cdots, u_{N_{x}-1, 0}, u_{0, 1}, \cdots, u_{N_{x}-1, 1}$, $\cdots, u_{0, N_y-1}, \cdots, u_{N_x-1, N_y-1})^T\}$. For ${\mathit{\boldsymbol{u}}}, {\mathit{\boldsymbol{v}}} \in \mathcal{U}_h$, we define

where ${\bar{v}}_{i, j}$ is the conjugate of $v_{i, j}$, and we introduce some operators

Let $(x_{i}, y_{j})\in \Omega_{h}$ be the Fourier collocation points. Using the interpolation polynomial $I_{\textbf{N}}{p}(x, y)$ to approximate ${p}_{\textbf{N}}(x, y)$ of the function ${p}(x, y)$, and we have

where $\mu = \frac{\pi}{L}$, and the coefficient

with $c_{k_1} = 1$ for $|k_1| < N_x/2$, $c_{k_2} = 1$ for $|k_2| < N_y/2$, $c_{k_1} = 2$ for $k_1 = \pm N_x/2$, and $c_{k_2} = 2$ for $k_2 = \pm N_y/2$.

Then the $-(-\Delta)^{\frac{\alpha}{2}}$ on $p(x_i, y_j)$ can be discretized by ^[8]

where

Especially for $N = N_x = N_y$, the negative definite and symmetric matrix $D^{s}$ is given

where $\textbf{F}$ is a orthogonal matrix with the entries

One can observe that the expression $(\textbf{F}\otimes \textbf{F})^{-1} \Lambda (\textbf{F}\otimes \textbf{F}) \textbf{p}$ can be efficiently calculated in $\mathcal{O}(N^2\log N)$ operations using the built-in function $\tt ifft2(\Lambda.* fft2(\textbf{p}))$ in Matlab.

3.2.   Construction of the conservative exponential scheme

By applying the Fourier pseudo-spectral to Eqs (2.2)–(2.4), we obtain

where $g_1({\mathit{\boldsymbol{p}}}, {\mathit{\boldsymbol{q}}}) = \beta({{\mathit{\boldsymbol{p}}}}^2+{{\mathit{\boldsymbol{q}}}}^2){{\mathit{\boldsymbol{p}}}}$, $g_2({\mathit{\boldsymbol{p}}}, {\mathit{\boldsymbol{q}}}) = \beta({{\mathit{\boldsymbol{p}}}}^2+{{\mathit{\boldsymbol{q}}}}^2){{\mathit{\boldsymbol{q}}}}$.

Let

then, systems (3.2)–(3.4) can be rewritten as

Let $\mathcal{V} = \tau \mathcal{S} \mathcal{D}$, then integrating systems (3.5) and (3.6) from $t^{n}$ to $t^{n+1}$, we can obtain

Then, we can use $f\big({\mathit{\boldsymbol{p}}}(t^n+\xi \tau), {\mathit{\boldsymbol{q}}}(t^n+\xi \tau), {\mathit{\boldsymbol{w}}}(t^n+\xi \tau)\big)$ and

to approximate $f(\hat{{\mathit{\boldsymbol{p}}}}^{n+\frac{1}{2}}, \hat{{\mathit{\boldsymbol{q}}}}^{n+\frac{1}{2}}, {{\mathit{\boldsymbol{w}}}}^{n+\frac{1}{2}})$ and

respectively, and obtain the following new scheme

where

Remark 3.1. It should be noted that the proposed schemes (3.9) and (3.10) are a three level method, therefore, we can calculate the ${\mathit{\boldsymbol{z}}}^{1}$ and ${\mathit{\boldsymbol{w}}}^{1}$ by using the ${\mathit{\boldsymbol{z}}}^0$ instead of ${ \hat{{\mathit{\boldsymbol{z}}}}}^{\frac{1}{2}}$ for the first step.

Theorem 3.1. Schemes (3.4) and (3.5) conserve the modified energy

To this end, we first give some lemmas.

Lemma 3.1. ^[7] Assuming that $g(x)$ is a sufficiently smooth function, in the vicinity of zero, with $0$ being a removable singularity, we can define $g(0)$ as the limit of $g(x)$ as $x$ approaches zero.

and for a matrix $\mathcal{A}$, the matrix valued function is defined by $g(\mathcal{A}) = \sum\limits_{k = 0}^{\infty}\frac{g^{k}(0)}{k!}\mathcal{A}^k.$

Lemma 3.2. ^[17] If $\mathcal{D}$ is a symmetric matrix and $\tau$ is a non-negative scalar, then the matrix $\mathcal{A}$ defined by $\mathcal{A} = \exp(\mathcal{V})^T \mathcal{D} \exp(\mathcal{V})-\mathcal{D}$ is a nilpotent matrix, where $\mathcal{V} = \tau \mathcal{S} \mathcal{D}$ and $\mathcal{S}$ is a skew-symmetric matrix.

Next, we show that the proposed scheme can preserve the modified energy.

Proof. If $\mathcal{D}$ is a singular matrix, let $\{\mathcal{D}_\epsilon\}$ be a series of symmetric and nonsingular matrices, which converge to $\mathcal{D}$ when $\epsilon\rightarrow 0$. We assume ${\mathit{\boldsymbol{z}}}^{n}_{\epsilon}, {\mathit{\boldsymbol{w}}}^n_{\epsilon}$ satisfy the perturbed scheme

where $\mathcal{V}_\epsilon = \tau \mathcal{S} \mathcal{D}_\epsilon$, $n = 0, 1, \cdots$. Then, we denote

and rewrite $E^n$ as

According to Eq (3.4), we have

On the other hand, from Eq (3.12), we have

This together with Eq (3.9), we can deduce

where

The last equality is based on Lemma 3.1 and the skew symmetry of the matrix $\mathcal{B}$. Thus, from Eq (3.9), when $\epsilon \rightarrow 0, {{\mathit{\boldsymbol{z}}}}_{\epsilon}^{n}\rightarrow {{\mathit{\boldsymbol{z}}}}^{n}, {{\mathit{\boldsymbol{w}}}}_{\epsilon}^{n}\rightarrow {{\mathit{\boldsymbol{w}}}}^{n}$, we obtain $E^{n+1} = E^{n}$.

If $\mathcal{D}$ is a nonsingular matrix, the proof is similar to above process, for brevity, we omit it. The proof is completed.

4.   Numerical example

Two examples are given to confirm the efficiency, conservation and stability of the proposed scheme. We use the formula

to evaluate the convergence rate, where $\tau_{i}, \text {error}_{i}, ({i} = 1, 2)$ are step sizes and $L^{\infty}$-norm errors with the step size $\tau_{i}$, respectively. The relative energy error is defined as

where $E^{n}$ and $E^{n}$ denote the conservation laws at $t_n = n\tau$.

Example 4.1. Considering the equation with

The exact solution is not given, we use

to obtain numerical errors.

In the follow-up of the article, the E-IEQ is the developed scheme in our work, CN-IEQ is a linearly implicit Crank-Nicolson IEQ scheme ^[1], E-AVF scheme is a fully-implicit exponential AVF scheme, and the E-SAV scheme ^[14] represents a linearly implicit exponential scheme based on the SAV approach. Yet the general, we first take $\alpha = 1.8$ to demonstrate the accuracy, efficiency and conservation properties of E-IEQ scheme. In Table 1, numerical errors and convergence orders implies that four schemes have second-order accuracy in the time. Noting that 'NaN' implies the CN-IEQ scheme cannot be implemented with a large time step, but the E-IEQ scheme can. Figure 1 shows that linearly implicit schemes are more efficient than the E-AVF scheme, and the E-IEQ scheme is most efficient of the four schemes.

In Figure 2, the relative errors of the modified energy and original energy of the equation are presented. From the results, we can see three linearly implicit schemes only can conserve the modified energy, the E-AVF scheme can preserve the inherent energy of the equation. Finally, we focus on the correlation between the fractional order $\alpha$ and the soliton's shape. Figure 3 displays numerical outcomes for various $\alpha$ values. It is evident that $\alpha$ has a significant impact on the soliton's form. As $\alpha$ shifts from 1.3 to 2, the soliton's appearance changes dramatically. The findings reveal that the solitary wave can endure stable propagation under finite initial circumstances, even after long periods of propagation.

Example 4.2. We consider the equation on $(x, y) = (-\pi, \pi)\times(-\pi, \pi)$ and has the form

with

where the amplitude $A = 1$, the wave numbers $\lambda_1 = \lambda_2 = 2$ and $\beta = 2$.

This example takes $\alpha = 1.4, 1.7, 2$. Figure 4 shows that the proposed scheme has second-order in time. We then run a long-time simulation till $T = 50$ and plot the relative energy deviation in Figure 5 which indicates that the new scheme can preserve the modified energy exactly in discrete scene.

5.   Conclusions

In this paper, we introduce a unique technique to solve the fractional nonlinear Schrödinger equation while conserving energy. We achieve this through a combination of the invariant energy quadratization method and the exponential time differencing method. Our proposed scheme is highly efficient and performs well even with large time steps during long-time computations. In discrete settings, the modified energy is precisely preserved. To demonstrate our approach, we present a numerical example. Furthermore, this technique can be applied to other fractional Hamiltonian PDEs, including the fractional nonlinear wave equation and the fractional Klein-Gordon-Schrödinger equation, by following a similar process.

Acknowledgments

The research is supported by the National Natural Science Foundation of China (No. 12001471).

Conflict of interest

The authors declare there is no conflict of interest.