On convergence of a structure preserving difference scheme for two-dimensional space-fractional nonlinear Schrödinger equation and its fast implementation

Abstract

In this paper we intend to construct a structure preserving difference scheme for two-dimensional space-fractional nonlinear Schrödinger (2D SFNS) equation with the integral fractional Laplacian. The temporal direction is discretized by the modified Crank-Nicolson method, and the spatial variable is approximated by a novel fractional central difference method. The mass and energy conservations and the convergence are rigorously proved for the proposed scheme. For 1D SFNS equation, the convergence relies heavily on the $L^{\infty }$-norm boundness of the numerical solution of the proposed scheme. However, we cannot obtain the $L^{\infty }$-norm boundness of the numerical solution by using the similar process for the 2D SFNS equation. One of the major significance of this paper is that we first obtain the $L^{\infty }$-norm boundness of the numerical solution and $L^2$-norm error estimate via the popular “cut-off” function for the 2D SFNS equation. Further, we reveal that the spatial discretization generates a block-Toeplitz coefficient matrix, and it will be ill-conditioned as the spatial grid mesh number M and the fractional order α increase. Thus, we exploit an linearized iteration algorithm for the nonlinear system, so that it can be efficiently solved by the Krylov subspace solver with a suitable preconditioner, where the 2D fast Fourier transform (2D FFT) is applied in the solver to accelerate the matrix-vector product, and the standard orthogonal projection approach is used to eliminate the drift of mass and energy. Extensive numerical results are reported to confirm the theoretical analysis and high efficiency of the proposed algorithm.

Introduction

Consider the two-dimensional space-fractional nonlinear Schrödinger equation

$$\mathbf{\text{i}}{u}_t-{(-\mathrm{\Delta })}^{\frac{\alpha }{2}}u+\beta |u{|}^{2}u=0,(x,y,t)\in {\mathbb{R}}^2\times (0,T],$$

$$u(x,y,0)=\phi (x,y),(x,y)\in {\mathbb{R}}^2,$$

where $1<\alpha \le 2$, ${\mathbf{\text{i}}}^2=-1$, and β is the nonlinearity strength (positive for a repulsive/defocusing interaction and negative for an attractive/focusing interaction) [7]. We assume that the functions u and φ are spatially compactly supported on the bounded domain Ω. The integral fractional Laplacian is defined as [25], [28]

$${(-\mathrm{\Delta })}^{\frac{\alpha }{2}}u(\mathit{x})=c_{2,\alpha }\mathrm{\text{P.V}}.\underset{{\mathbb{R}}^2}{\iint }\frac{u(\mathit{x})-u(\mathit{\xi })}{|\mathit{x}-\mathit{\xi }{|}^{2+\alpha }}d\mathit{\xi },c_{2,\alpha }=\frac{2^{\alpha }\mathrm{\Gamma }(\frac{\alpha +2}{2})}{\pi |\mathrm{\Gamma }(-\alpha /2)|},$$

where P.V. stands for Cauchy principle value, Γ is the Gamma function, $|\mathit{x}-\mathit{\xi }|$ denotes the Euclidean distance between points $\mathit{x}=(x,y)$ and $\mathit{\xi }=({\xi }_1,{\xi }_2)$. Besides, we postulate that the functions u and φ are smooth enough, so that our proposed scheme can attain the desired convergence order.

Schrödinger equation is famous for its applications in quantum mechanics, it is derived by using Feynman path integrals over Brownian process. Recently, Laskin [23] indicated that the path integral over the Lévy process allows us to develop the generalization of the quantum mechanics. Obviously, the Brownian process leads to the classical Schrödinger equation, and the symmetric α-stable Lévy process [26] leads to the space-fractional Schrödinger equation. Such important attempts have attracted growing attentions in many fields, such as nonlinear optics [24], propagation dynamics [42], and water wave dynamics [21]. The theoretical consequences of space-fractional Schrödinger equation, including the well-posedness of the smooth solution, as well as conservations of invariants have been studied, readers can refer [14], [15] and references therein.

It is known to us, the 2D SFNS equation has some physical invariants over the time, such as the mass and energy [22], i.e.,

$$\mathcal{M}(t)=\underset{{\mathbb{R}}^2}{\iint }|u(x,y){|}^{2}dxdy=\mathcal{M}(0),\mathcal{E}(t)=\underset{{\mathbb{R}}^2}{\iint }[\overline{u}{(-\mathrm{\Delta })}^{\frac{\alpha }{2}}u(x,y,t)-\frac{\beta }{2}|u(x,y,t){|}^4]dxdy=\mathcal{E}(0),$$

where $\overline{u}$ is the complex conjugate of u. Numerous theoretical and experimental results [13], [16] show that the structure preserving algorithms can preserve the intrinsic physical invariants of original conservative system and often have good numerical properties, such as the excellent long-time behavior, the linear error growth, and smaller amplitude and so on. In the past decades, the invariant-preserving integrator has attracted growing interests due to it preserves as much as possible the invariants of underlying differential systems. The well-known invariant-preserving integrators for conservative systems are about the averaged vector field (AVF) methods [29], the discrete variational derivative (DVD) methods [11] and the Hamiltonian boundary value methods (HBVMs) [2], [3], [4] and so on.

For SFNS equation in 1D space, the convergence analysis strictly relies on the $L^{\infty }$-norm boundness of numerical solution of the numerical scheme [34], [38]. In fact the proof of $L^{\infty }$-norm boundness of numerical solution relies heavily on not only the discrete conservative property but also the discrete version of the Sobolev inequality in 1D space

$${\Vert u\Vert }_{\infty }\le \mathcal{C}{\Vert u\Vert }_{H^{\frac{\alpha }{2}}},\forall u\in H_0^{\frac{\alpha }{2}}(\mathrm{\Omega }),\mathrm{\Omega }\subset \mathbb{R},$$

which immediately implies an a priori uniform $L^{\infty }$-norm boundness for numerical solution from the discrete energy formulation. However, we cannot obtain the $L^{\infty }$-norm boundness of numerical solution for the 2D SFNS equation by using the similar ways.

At present, there are some numerical schemes have been studied for the 2D SFNS equation. Zhao et al. [41] proposed a compact alternating direction implicit difference scheme and obtained the associated convergence result. Wang and Huang [35] constructed a split-step alternating direction implicit difference scheme, and the authors discussed the convergence analysis. Wang et al. [37] established a split-step spectral Galerkin method and analyzed the convergence of the proposed scheme. Zhang et al. [43] considered a compact implicit integration factor method. Nevertheless, the aforementioned numerical schemes cannot prove to be conservative for the 2D SFNS equation in theoretical analysis.

Recently, Wang and Huang [33] proposed a structure preserving difference scheme for the 1D SFNS equation, the authors adopted a “cut-off” function to truncate the nonlinearity to a global Lipschitz function with compact support on domain Ω for 1D SFNS equation, and authors obtained the $L^2$-norm convergence consequence by assuming $\tau \lesssim h$, where τ and h is the time and space step size, respectively. More recently, Yin et al. [39] obtained the same theoretical consequence by extending this technique to the 2D SFNS equation with Riesz derivative. It is worth noting that Wang et al. [36] successful eliminated the time-space step ratio $\tau \lesssim h$ by using the “lifting” technique (i.e., by taking the discrete $L^2$-norm of both sides of the “error” equations, the estimate of the higher order difference quotient of the error functions can be obtained from the estimate of the lower order difference quotient of the error functions and local truncation error functions), the authors obtained a unconditional convergence result for the 2D classical nonlinear Schrödinger equation. As a matter of fact, the “lifting” technique relies on a discrete Sobolev inequality in 2D space

$${\Vert u\Vert }_{\infty }\le \mathcal{C}{\Vert u\Vert }^{\frac{1}{2}}(|u{|}_{H^1}+\Vert u{\Vert )}^{\frac{1}{2}},\forall u\in H_0^1(\mathrm{\Omega }),\mathrm{\Omega }\subset {\mathbb{R}}^2.$$

However, we cannot apply the “lifting” technique to further obtain the unconditional convergence for 2D SFNS equation, due to the lack of a similar inequality of (1.4).

In the current paper, we construct a novel structure preserving difference scheme for the 2D SFNS equation, where the temporal derivative is approximated by the modified Crank-Nicolson method, and the 2D fractional Laplacian is discretized by a second-order fractional central difference operator [18]. With the help of a “cut-off” function, we prove the $L^2$-norm convergence consequence and the $L^{\infty }$-norm boundness of numerical solution. In numerical implementations, the preconditioned conjugated gradient is chosen as a linear solver of the fixed iteration for the nonlinear discrete system, where the 2D FFT algorithm is used in the linear solver to accelerate the matrix-vector product. For invariants preservation, we adopt the standard projection method to eliminate the drift of mass and energy of the conservative scheme. Extensive numerical comparisons are provided to verify the feasibility and efficiency of the proposed scheme.

The rest of the paper is organized as follows. In Section 2, we introduce the spatial discretization and associated fractional Sobolev inequality. Next, we construct a structure preserving difference scheme for the 2D SFNS equation in Section 3. In Section 5, we proposed some fast algorithms to improve our computation efficiency. In Section 4, we are given some useful lemmas and the theoretical consequence, including the conservation laws and convergence analysis. In Section 6, we provide some numerical example to verify our numerical scheme. Finally, we draw some conclusions in Section 7.