Linearly Implicit High-Order Exponential Integrators Conservative Runge–Kutta Schemes for the Fractional Schrödinger Equation

: In this paper, a family of high-order linearly implicit exponential integrators conservative schemes is constructed for solving the multi-dimensional nonlinear fractional Schrödinger equation. By virtue of the Lawson transformation and the generalized scalar auxiliary variable approach, the equation is ﬁrst reformulated to an exponential equivalent system with a modiﬁed energy. Then, we construct a semi-discrete conservative scheme by using the Fourier pseudo-spectral method to discretize the exponential system in space direction. After that, linearly implicit energy-preserving schemes which have high accuracy are given by applying the Runge–Kutta method to approximate the semi-discrete system in temporal direction and using the extrapolation method to the nonlinear term. As expected, the constructed schemes can preserve the energy exactly and implement efﬁciently with a large time step. Numerical examples conﬁrm the constructed schemes have high accuracy, energy-preserving, and effectiveness in long-time simulation. 65M06; 65M70


Introduction
In this paper, we mainly consider the fractional nonlinear fractional Schrödinger (NLS) equation in the following form [1,2] where i is the imaginary unit root, u(x, t) satisfy the periodic boundary, x = [−L, L] d ⊂ R d (d = 2 or 3), u 0 (x) is the smooth initial function, the parameter γ is a dimensionless interaction constant which describes the strength of short-range (or local) nonlinear interactions. The fractional Laplacian operator (−L) where u(ξ) = Ω u(x)e −iξx dx denotes the Fourier transform of u(x).
Based on the conjugacy of (−L) α 2 , the system (1) with periodic boundary has two invariants [3], that is the Hamiltonian energy H(u) = ((−L) α 2 u, u) − γ 2 (|u| 2 , |u| 2 ) = cosnt. (3) and the mass where (φ, ψ) = Ω φψdx,ψ is the conjugate of ψ. By setting u = φ + iψ, system (1) can be reformulated as an infinite-dimension canonical Hamiltonian system [4] dz dt = S −1 δH δz , with z = (φ, ψ) T , S = 0 1 −1 0 , where δH/δz is the vector of variational derivatives with respect to z. The equation is first introduced by the physicist, it is an important model in quantum mechanics and has been implemented in some scientific engineering fields [5][6][7]. We all know that it is difficult obtain exact solution of the fractional NLS equation for the reason that the solution contains special function. Therefore, different numerical methods were proposed and analyzed for solving the equation and approximating fractional derivative in recent decades, such as finite difference methods [8][9][10], spectral methods [4,11,12], finite element methods [13,14] and so on. The prior research substantiated the finding that energy conservation plays a crucial role in studying the solution properties for conservative system. The method can conserve some intrinsic properties of the dynamical system, which we called the structure-preserving algorithm, and scholars have developed structure-preserving schemes for system (1). For example, Wang et al. developed some conservative finite difference schemes for solving the fractional NLS-type Equation [15,16], and then some scholars further discussed unconditional L ∞ -norm convergence error estimates [17,18]. In Ref. [19], Duo et al. constructed an efficient Fourier pseudo-spectral scheme which can inherit the mass and energy for system (1) in two-dimensional case. There are many works related to it, the readers can find in Refs. [20,21].
Unfortunately, a majority of conservative schemes only achieve second order temporal accuracy, these schemes can not get satisfactory numerical solutions with a large time step. Over the past decade, many numerical methods were developed and can be extend to construct high-order numerical schemes to inherit the energy of the fractional NLS equation, such as the Hamiltonian Boundary Value Method (HBVM) [22,23], the averaged vector flied (AVF) method [24,25], and the Runge-Kutta (RK) method [26]. However, these schemes are fully implicit, one has to use the nonlinear iterative to solve algebraic systems and thus it is time consuming in numerical modeling. To reduce the computational complexity, explicit schemes that can preserve a modified energy were constructed to solve the fractional NLS system based on the energy quadratization method [27,28]. Such schemes are very efficient but have poor stability in practical computation.
In the past years, the exponential integrators method was used to construct structurepreserving schemes for conservative (dissipation) Equations [29][30][31]. The proposed schemes can approximate linear part of the system, and provide satisfactory stability and high accuracy for highly oscillatory systems. Recently, scholars developed conservative exponential integrators scheme for the NLS equation based on the scalar auxiliary variable (SAV) approach [32,33]. The energy-preserving algorithms mentioned above are either fully-implicit or linearly-implicit so that nonlinear or linear solvers are required to obtain the conservative solutions. Over all, developing high order and efficiency conservative exponential integrators schemes for solving Equation (1) is still be considered.
The generalized scalar auxiliary variable (GSAV) approach is build upon the SAV approach [34][35][36][37], and has been used to construct efficient energy stable methods for solving some gradient flow models [38,39]. The auxiliary variables of the newly method not just limited to the square root function and can be solved more efficiently than the traditional SAV approach. The method is also extended to construct linearly implicit highorder conservative schemes for Hamiltonian partial differential Equations [40], the resulted scheme needs very small time step and it is still a computationally expensive procedure in long time numerical simulation. Whether the GSAV approach and the exponential integrators can be extended to construct efficient conservative schemes that conserve the energy of the fractional NLS equation has not been studied. With this aim in mind, the paper presents a new method for constructing a high-order linearly implicit conservative scheme for system (1) by combing the GSAV method and the RK method. Unlike the traditional GSAV schemes, such a reformation allows the approximation of entire nonlinear terms by extrapolation and retains the energy conservation, whereas the auxiliary variable in GSAV schemes has to be discretized implicitly. The proposed linearly-implicit schemes with high accuracy in the time direction can be solved efficiently, inherit modified energy of the system, and have better stability results than the non-exponential scheme in the practical numerical simulation.
The outline of this article is organized as follows. In Section 2, we derive an equivalent exponential system with a modified energy for the fractional NLS equation based on the idea of the Lawson transformation and GSAV method. In Section 3, a class of linearlyimplicit high-order exponential integrators schemes are constructed, and the conservation property and fast solver of the developed schemes are discussed. In Section 4, Some numerical examples are given to confirm the theoretical results. We draw some conclusions in Section 5.

The GSAV Approach for the NLS Equation
According to the basic idea of the GSAV method, we define an auxiliary variable as follows where the function G is invertible, and Formula (3) can be rewritten as Then the nonlinear functional can be transformed into the following equivalent formulation By taking the derivative of (5) with respect to t, we obtain where Re(•) denotes the real part of •, and we can deduce the following equality Then, system (1) can be reformulated as with the consistent initial condition u 0 = u(x, 0), r 0 = G(u 0 ). Noting that (−L) α 2 u, u ∈ R, we can prove system (10) and (11) conserve the mass and the modified energy, namely Theorem 1. The continuous system (10) and (11) can preserve the mass and the modified energy, namely where the the mass functional M(t) and the modified energy functional H(t) are given by Proof. Based on the definition of the mass functional, we can derive where Im(•) are the imaginary parts of •. Similarly, noting that (−L) The proof is completed.
Furthermore, based on the Lawson transformation method [41], by setting u = exp(−i(−L) α 2 t)v, and multiplying both sides of (10) by the operator exp(i(−L) The derivation above means that system (10) and (11) and we can prove system (17) and (18) also preserve the modified energy.

Remark 1.
By taking the different function G(x), we can derive various SAV approaches [38]. Without losing generality, we set G( the newly developed GSAV approach is equivalent to the SAV approach [36,37,42].

Construction of the Energy-Preserving Schemes
In this section, we will present a detailed construction of linearly implicit exponential integrators conservative schemes for system (17) and (18) in two-dimensional case. Similarly, it can be extended to three-dimensional case.

Fourier Pseudo-Spectral Approximation of Spatial Derivatives
Under the homogeneous Dirichlet boundary condition, the fractional Laplacian operator is equivalent to the Riesz derivatives, the spatial discretization of operators can be used by finite difference methods [16,18]. The paper considers the fractional NLS problem with periodic boundary conditions. Therefore, the Fourier-pseudo spectral method with high order accuracy and fast calculation is employed for the spatial discretization to approximate the fractional Laplacian operator.
For positive even integers N = N x = N y and N, the step sizes in space are defined by h = h x = h y := 2L N , and the time step τ := T N . Then, we set The discrete inner product and L ∞ -norm are given as In practical computation, the Laplace operator (−L) α 2 can be approximated by Fourierpseudo spectral method, namely [4] (−L) where µ x = µ y = π L , F N is Fourier transform in discrete scene, and the corresponding transform and the spectral differentiation matrix can be expressed as By setting Λ = I N ⊗ Λ 1 + Λ 2 ⊗ I N , we can get therefore, it implies that the exponential factor can be implemented by using the FFT technique. Then, the Fourier pseudo-spectral method is applied to system (17) and (18) Noting that u = exp(iDt)v, we can derive a conservation property as follows. (21) and (22) can inherit the mass and modified energy in semi-discrete scene, namely

Theorem 2. System
where the the semi-discrete mass functional M(t) and modified energy functional E (t) are given by Proof. The proof is similar to Theorem 1, we omit it.

Fully Discrete Energy-Preserving Schemes
In this subsection, we construct fully discrete schemes by using the symplectic RK method to approximate system (21) and (22) in time.
First, we define t n = nτ, n = 0, 1, · · · , N. Let a ij , b i and c i = m ∑ i=1 a ij be real numbers.
For one-step interval [t n , t n+1 ], the m-stage RK method for system (21) and (22), and further apply the extrapolation technique to the nonlinear term, we can derive the following fully discrete scheme whereū in is an extrapolation approximation to u(t n + c i τ) (i = 1, 2, · · · , m) and can reach order of O(τ m+1 ) [43,44]. Then u n+1 , r n+1 can be updated by Theorem 3. The fully-discrete system (26)- (29) can preserve a modified energy if the coefficients of the RK method satisfy a ij b i + a ji b j = b i b j , namely Proof. According to (28), we have From (26)-(28), we derive Then, we can obtain from (26)-(28) that in G ū in r i were used, and we can derive The proof is completed.
Remark 2. Theorem 2 shows the semi-discrete scheme can preserve the mass, but the developed fully-discrete scheme (26)- (29) can not conserve the mass for the reason that the nonlinear terms of system (21) and (22) are explicitly discretized.

Remark 3.
According to Remark 1, G(x) can be taken different function, therefor we can obtain a class of energy-preserving schemes in practical numerical simulation.

•
The tanh SAV scheme: In this scheme, we select G(x) = tanh( x C ), where C is a positive constant to make r = tanh Ω |u(x, t)| 4 dx not too close to ±1 numerically since tanh −1 (r = ±1) → ∞. Thus, we set and prove that the corresponding scheme can preserve the following discrete energy • The exponential SAV scheme: The exponential function is a special function that can keep the range constant positive. Thus, we define an exponential scalar auxiliary variable Then, we derive the following discrete energy conservation law E n = Du n , u n h + γC 2 ln(r n ).

Fast Solvers for the Proposed Schemes
In particular, for i, j = 1, 2, · · · , m, if the coefficients of an m-stage RK method satisfy the RK method is symplectic, and we called the formula (32) is the RK symplecticity conditions. Without losing generality, we choose the 2m-4th and 2m-6th symplectic RK method, which are displayed by the following Butcher tabular [45] Then, taking m = 2 as an example, we propose a fast solver for the third order linearly implicit exponential integrators energy-preserving (LI-EI3) scheme. By setting in G ū in r in , (i = 1, 2), and rewrite and the extrapolation technique for the nonlinear terms can be expressed as [43,44] u n1 = (2 Then, we can derive Based on the equality of (28), (34) and (35), we have where Multiplying both sides of (38) and (39) with D and computing discrete inner product with α * 1 and α * 2 , respectively, we can deduce with We can solve Re Du n1 , α * 1 h , Re Du n2 , α * 2 h T based on the linear system (40) and (41), and l i , k i and u ni , i = 1, 2 can be updated from (33)-(37), respectively. Subsequently, we can compute u n+1 and r n+1 by using (28). For m = 3, the extrapolation technique for the nonlinear terms can be expressed as follows [43,44] Similarly argument, we can derive a fast solver for the fourth order linearly implicit exponential integrators energy-preserving (LI-EI4) scheme. For simplicity, we omit it.

Numerical Experiments
In this section, we implement our simulations by using Matlab R2018a software on a computer which composed of Intel(R) Core(TM) i7-9750H, 2.6 GHz CPU machine with 16 GB RAM and display some numerical experiments to verify the conservation property and efficiency of the constructed schemes in this section. For simplicity, the relative modified energy error can be defined by where E n represents the energy at t = nτ. The convergence rate of the proposed schemes can be obtained by the formula Rate = ln(error 1 /error 2 )/ln(τ 1 /τ 2 ), where τ i is the time step, error i , (i = 1, 2) represents the L ∞ -norm errors at τ i , and we use the '*' to represent that the convergence order cannot be obtained. In addition, we also compare the proposed schemes with existing scheme in computing efficiency, accuracy and conservation. Thus, we define • LI-EI-i (i = 3 or 4): The paper constructs third and fourth order energy-preserving schemes by using the Runge-Kutta methods shown in Table 1. • LI-4: A fourth order linearly implicit conservative RK method is based on the GSAV approach [40]. • LI-EI-2: A second order linearly-implicit exponential time differencing conservative scheme is developed in Ref. [33]. • FI-EI-4: A fourth order fully-implicit conservative exponential time differencing method is presented in Ref. [32].
In practical calculation, we take h = π 8 such that the spatial discretization errors are negligible to test the temporal discretization errors of the proposed energy-preserving methods. The discrete L ∞ -norm errors and the corresponding convergence orders at T = 1 are displayed in Table 2.
Numerical results indicate that the LI-EI-2 scheme only has second order accuracy in time. The proposed LI-EI-3 scheme has third order accuracy in time, and the LI-EI-4 scheme, the LI-4 scheme and the FI-EI-4 scheme has fourth order accuracy in time. We should note that 'NaN' represents the LI-4 scheme can not be implemented with τ = 1 50 , but the LI-EI-4 scheme can run, which implies the proposed exponential integrators schemes can be implemented with a large time step, and the step ratio is smaller than in the non-exponential scheme, and the stability result is better than the LI-4 scheme. Then, we compare five schemes in computational costs. The CPU time of five scheme using the different time step with T = 50 is shown in Figure 1. We note that the FI-EI-4 scheme is fully-implici and needs nonlinear iterations to solve the equations, therefor, it consumes the most CPU time. The proposed schemes and the LI-EI-2 scheme are linearly implicit schemes and can be solved efficiently, they enjoy the same computational advantages as the non exponential LI-4 scheme. In summary, it is preferable to develop efficient exponential integrator schemes which have better numerical stability and conserve discrete energy for the fractional NLS equation. Last but not least, we plot the modified energy deviation in a long time simulation, corresponding to all five schemes. As is shown in Figure 2, the proposed schemes also can conserve the modified energy exactly.
First, we test the accuracy of the proposed schemes in time, the mesh size is taken as h = 2π 16 so that the spatial error is negligible. The errors and convergence rates in L ∞ -norm are presented in Table 3, the results indicate that the LI-EI-3 scheme and LI-EI-4 scheme have third and fourth order convergence rates in time, respectively. In addition, numerical results also demonstrate that the parameter α affects the numerical solution error, namely, the error of numerical solution increases with the increase of α. Then, we study the conservation properties of two schemes at a large time T = 20. The relative modified energy errors are shown in Figure 3, which demonstrates the newly developed schemes can preserve the energy exactly. LI-EI-3 scheme

Conclusions
This paper presents a class of novel exponential integrators Runge-Kutta schemes for solving the nonlinear Schrödinger equation. The energy conservation property and high efficiency of the proposed schemes are supported by theoretical analysis and numerical results. Similar conservative schemes are also constructed to solve others Hamiltonian partial differential equations.