Energy preserving methods for nonlinear Schr\"odinger equations

TL;DR

This paper analyzes energy-preserving numerical methods for nonlinear Schrödinger equations, proving the order of the relaxation method and introducing a generalized version for various nonlinearities, supported by numerical simulations.

Contribution

It provides a rigorous proof of the relaxation method's order and proposes a generalized energy-preserving method for different nonlinearities.

Findings

Relaxation method has order 2 for cubic nonlinear Schrödinger equations.

Generalized relaxation method effectively handles various power law nonlinearities.

Numerical simulations demonstrate the efficiency of the proposed methods.

Abstract

This paper is concerned with the numerical integration in time of nonlinear Schr\"odinger equations using different methods preserving the energy or a discrete analog of it. The Crank-Nicolson method is a well known method of order 2 but is fully implicit and one may prefer a linearly implicit method like the relaxation method introduced in [10] for the cubic nonlinear Schr{\"o}dinger equation. This method is also an energy preserving method and numerical simulations have shown that its order is 2. In this paper we give a rigorous proof of the order of this relaxation method and propose a generalized version that allows to deal with general power law nonlinearites. Numerical simulations for different physical models show the efficiency of these methods.

Equations

$$i{\partial }_t\phi (t,x)=\left(-\frac{1}{2}\Delta +V(x)+\beta |\phi {|}^{2\sigma }(t,x)+\lambda \left(U\ast |\phi (t,\cdot ){|}^2\right)(x)-\Omega .R\right)\phi (t,x),$$

$$i{\partial }_t\phi (t,x)=\left(-\frac{1}{2}\Delta +V(x)+\sum _{k=1}^K{\beta }_k|\phi {|}^{2{\sigma }_k}(t,x)+\lambda \left(U\ast |\phi (t,\cdot ){|}^2\right)(x)-\Omega .R\right)\phi (t,x),$$

$$E(\phi )={\int }_{{\mathbb{R}}^d}\left(\frac{1}{4}\Vert \nabla \phi {\Vert }^2+\frac{1}{2}V|\phi {|}^2+\frac{\beta }{2\sigma +2}|\phi {|}^{2\sigma +2}+\frac{\lambda }{4}(U\ast |\phi {|}^2)|\phi {|}^2-\frac{\Omega }{2}\overline{\phi }R\phi \right)dx,$$

$$\mathrm{Re}{\int }_{{\mathbb{R}}^d}|\phi {|}^{2\sigma }\phi \overline{{\partial }_t\phi }dx=\frac{1}{2\sigma +2}\frac{d}{dt}{\int }_{{\mathbb{R}}^d}|\phi {|}^{2\sigma +2}dx,\forall \sigma \ge 0.$$

$$i\frac{{\phi }_{n+1}-{\phi }_n}{\delta t}=-\frac{1}{2}\Delta \frac{{\phi }_{n+1}+{\phi }_n}{2},$$

$$i\frac{{\phi }_{n+1}-{\phi }_n}{\delta t}=\left(-\frac{1}{2}\Delta +V+\beta g({\phi }_n,{\phi }_{n+1})+\lambda \left(U\ast \left(\frac{|{\phi }_{n+1}{|}^2+|{\phi }_n{|}^2}{2}\right)\right)-\Omega .R\right)\frac{{\phi }_{n+1}+{\phi }_n}{2},$$

$$\beta {\int }_{{\mathbb{R}}^d}g({\phi }_n,{\phi }_{n+1})\left(|{\phi }_{n+1}{|}^2-|{\phi }_n{|}^2\right)dx,$$

$$g({\phi }_n,{\phi }_{n+1})\left(|{\phi }_{n+1}{|}^2-|{\phi }_n{|}^2\right)=G({\phi }_{n+1})-G({\phi }_n).$$

$$i\frac{{\phi }_{n+1}-{\phi }_n}{\delta t}$$

$$=$$

$${\phi }_{n+1}={\Phi }_{\delta t}^{\mathrm{CN}}({\phi }_n),$$

$$\frac{\beta }{\sigma +1}\frac{|{\phi }_{n+1}{|}^{2\sigma +2}-|{\phi }_n{|}^{2\sigma +2}}{|{\phi }_{n+1}{|}^2-|{\phi }_n{|}^2}=\frac{\beta }{\sigma +1}\sum _{k=0}^{\sigma }|{\phi }_{n+1}{|}^{2k}|{\phi }_n{|}^{2(\sigma -k)},$$

$$E_{\mathrm{CN}}(\phi )=E(\phi ),$$

$$i{\partial }_t\phi (t,x)=-\frac{1}{2}\Delta \phi (t,x)+\beta |\phi (t,x){|}^2\phi (t,x),$$

$$\left\{\begin{array}[]{l}\Upsilon=|\varphi(t,x)|^{2},\\ \displaystyle i\partial_{t}\varphi(t,x)=-\frac{1}{2}\Delta\varphi(t,x)+\beta\Upsilon\varphi(t,x).\end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle\frac{\Upsilon_{n+1/2}+\Upsilon_{n-1/2}}{2}=|\varphi_{n}|^{2},\\ \displaystyle i\frac{\varphi_{n+1}-\varphi_{n}}{\delta t}=\left(-\frac{1}{2}\Delta+\beta\Upsilon_{n+1/2}\right)\frac{\varphi_{n+1}+\varphi_{n}}{2},\end{array}\right.$$

$$E_{\mathrm{rlx}}(\phi ,{\rm Y})=\frac{1}{4}{\int }_{{\mathbb{R}}^d}\Vert \nabla \phi {\Vert }^{2}dx+\frac{\beta }{2}{\int }_{{\mathbb{R}}^d}{\rm Y}|\phi {|}^{2}dx-\frac{\beta }{4}{\int }_{{\mathbb{R}}^d}{{\rm Y}}^{2}dx.$$

$$\displaystyle\mathrm{Re}\left(\Upsilon_{n+1/2}\frac{\varphi_{n+1}+\varphi_{n}}{2}\overline{\varphi_{n+1}-\varphi_{n}}\right)=\frac{\Upsilon_{n+1/2}}{2}\Big{(}|\varphi_{n+1}|^{2}-|\varphi_{n}|^{2}\Big{)}.$$

$$\begin{array}[]{ll}\displaystyle\Upsilon_{n+1/2}\Big{(}|\varphi_{n+1}|^{2}-|\varphi_{n}|^{2}\Big{)}&\displaystyle=\Upsilon_{n+1/2}\Big{(}|\varphi_{n+1}|^{2}-|\varphi_{n}|^{2}\Big{)}+\Upsilon_{n-1/2}\Big{(}|\varphi_{n}|^{2}-|\varphi_{n}|^{2}\Big{)}\\ &=\displaystyle\Big{(}\Upsilon_{n+1/2}|\varphi_{n+1}|^{2}-\Upsilon_{n-1/2}|\varphi_{n}|^{2}\Big{)}-{\Big{(}\Upsilon_{n+1/2}-\Upsilon_{n-1/2}\Big{)}|\varphi_{n}|^{2}}.\end{array}$$

$$\displaystyle{\Big{(}\Upsilon_{n+1/2}-\Upsilon_{n-1/2}\Big{)}|\varphi_{n}|^{2}}=\frac{\Big{(}\Upsilon_{n+1/2}\Big{)}^{2}-\Big{(}\Upsilon_{n-1/2}\Big{)}^{2}}{2}.$$

$$\mathrm{Re}\left(\Upsilon_{n+1/2}\frac{\varphi_{n+1}+\varphi_{n}}{2}\overline{\varphi_{n+1}-\varphi_{n}}\right)=\Big{(}\Upsilon_{n+1/2}|\varphi_{n+1}|^{2}-\Upsilon_{n-1/2}|\varphi_{n}|^{2}\Big{)}+\frac{\Big{(}\Upsilon_{n+1/2}\Big{)}^{2}-\Big{(}\Upsilon_{n-1/2}\Big{)}^{2}}{2},$$

$$E_{\mathrm{rlx}}(\phi ,|\phi {|}^2)=E(\phi ).$$

$$\left\{\begin{array}[]{l}i\partial_{t}\varphi+\dfrac{1}{2}\Delta\varphi=\beta\Upsilon\varphi,\\ \partial_{t}\Upsilon=2\mathrm{Re}(\overline{v}\varphi),\\ i\partial_{t}v+\dfrac{1}{2}\Delta v=\beta(\partial_{t}\Upsilon\varphi+\Upsilon\partial_{t}\varphi).\end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle\Phi_{n+\frac{1}{2}}=\Upsilon_{n+\frac{1}{2}}\left(\frac{\varphi_{n+1}+\varphi_{n}}{2}\right),\\[5.69046pt] \displaystyle\Xi_{n+\frac{1}{2}}=2\mathrm{Re}\left(v_{n+\frac{1}{2}}\left(\overline{\frac{\varphi_{n+1}+\varphi_{n}}{2}}\right)\right),\\[5.69046pt] \displaystyle V_{n+\frac{1}{2}}=\left(\frac{\Upsilon_{n+\frac{3}{2}}+\Upsilon_{n-\frac{1}{2}}}{2}\right)\left(\frac{v_{n+\frac{3}{2}}+2v_{n+\frac{1}{2}}+v_{n-\frac{1}{2}}}{4}\right)\\[8.5359pt] \displaystyle\qquad\qquad+2\mathrm{Re}\left(v_{n+\frac{1}{2}}\left(\overline{\frac{\varphi_{n+1}+\varphi_{n}}{2}}\right)\right)\left(\frac{\varphi_{n+2}+\varphi_{n+1}+\varphi_{n}+\varphi_{n-1}}{4}\right).\end{array}\right.$$

$$\left\{\begin{array}[]{lr}\displaystyle i\frac{\varphi_{n+2}-\varphi_{n+1}}{\delta t}+\dfrac{1}{2}\Delta\left(\frac{\varphi_{n+2}+\varphi_{n+1}}{2}\right)=\beta\Phi_{n+\frac{3}{2}},&(\ref{systran}.b)\\[8.5359pt] \displaystyle\frac{\Upsilon_{n+\frac{3}{2}}-\Upsilon_{n-\frac{1}{2}}}{2\delta t}=\Xi_{n+\frac{1}{2}},&(\ref{systran}.a)\\[8.5359pt] \displaystyle i\frac{v_{n+\frac{3}{2}}-v_{n-\frac{1}{2}}}{2\delta t}+\dfrac{1}{2}\Delta\left(\frac{v_{n+\frac{3}{2}}+2v_{n+\frac{1}{2}}+v_{n-\frac{1}{2}}}{4}\right)=\beta V_{n+\frac{1}{2}}.&(\ref{systran}.c)\end{array}\right.$$

$$X_n:=\left({{\rm Y}}_{n-\frac{1}{2}},{{\rm Y}}_{n+\frac{1}{2}},{\phi }_{n-1},{\phi }_n,{\phi }_{n+1},v_{n-\frac{1}{2}},v_{n+\frac{1}{2}}\right).$$

$$X_n\mapsto X_{n+1}.$$

$$\left\{\begin{array}[]{l}\displaystyle v_{n-\frac{1}{2}}=\frac{\varphi_{n}-\varphi_{n-1}}{\delta t},\\ \displaystyle v_{n+\frac{1}{2}}=\frac{\varphi_{n+1}-\varphi_{n}}{\delta t},\\ \displaystyle\Upsilon_{n+\frac{1}{2}}+\Upsilon_{n-\frac{1}{2}}=2|\varphi_{n}|^{2},\\ \displaystyle i\frac{\varphi_{n}-\varphi_{n-1}}{\delta t}=\left(-\dfrac{1}{2}\Delta+\beta\Upsilon_{n-1/2}\right)\frac{\varphi_{n}+\varphi_{n-1}}{2},\\ \displaystyle i\frac{\varphi_{n+1}-\varphi_{n}}{\delta t}=\left(-\dfrac{1}{2}\Delta+\beta\Upsilon_{n+1/2}\right)\frac{\varphi_{n+1}+\varphi_{n}}{2},\end{array}\right.$$

$$\left\{\begin{array}[]{ll}\displaystyle\Upsilon_{\frac{1}{2}}&\displaystyle=2|\varphi_{0}|^{2}-\Upsilon_{-\frac{1}{2}},\\ \displaystyle\varphi_{-1}&\displaystyle=\left(2-i\delta t\left(-\dfrac{1}{2}\Delta+\beta\Upsilon_{-1/2}\right)\right)^{-1}\left(2+i\delta t\left(-\dfrac{1}{2}\Delta+\beta\Upsilon_{-1/2}\right)\right)\varphi_{0},\\ \displaystyle\varphi_{1}&\displaystyle=\left(2+i\delta t\left(-\dfrac{1}{2}\Delta+\beta\Upsilon_{1/2}\right)\right)^{-1}\left(2-i\delta t\left(-\dfrac{1}{2}\Delta+\beta\Upsilon_{1/2}\right)\right)\varphi_{0},\\ \displaystyle v_{-\frac{1}{2}}&\displaystyle=(\varphi_{0}-\varphi_{-1})/\delta t,\\ \displaystyle v_{\frac{1}{2}}&\displaystyle=(\varphi_{1}-\varphi_{0})/\delta t.\end{array}\right.$$

$$A=(i-\delta t\Delta /4{)}^{-1}(i+\delta t\Delta /4)\text{ and }B=(i-\delta t\Delta /4{)}^{-1},$$

$$\mathcal{C}=\left(\begin{array}{ccccccc}I& 0& 0& 0& 0& 0& 0\\ 0& I& 0& 0& 0& 0& 0\\ 0& 0& B& 0& 0& 0& 0\\ 0& 0& 0& B& 0& 0& 0\\ 0& 0& 0& 0& B& 0& 0\\ 0& 0& 0& 0& 0& B& 0\\ 0& 0& 0& 0& 0& 0& B\end{array}\right).$$

Full text

AMS Classification: 35Q41, 81Q05, 65M70

Keywords. Nonlinear Schrödinger equation, Gross-Pitaevskii equation, numerical methods, relaxation methods.

1. Introduction

The nonlinear Schrödinger equation (NLSE) is a fairly general dispersive partial differential equation arising in many areas of physics and chemistry [1, 2, 16, 23, 12]. One of the most important application of the NLSE is for laser beam propagation in nonlinear and/or quantum optics and there it is also known as parabolic/paraxial approximation of the Helmholtz or time-independent Maxwell equations [1, 2, 16, 12]. In the context of the modeling of the Bose-Einstein condensation (BEC), the nonlinear Schrödinger equation is known as the Gross-Pitaevskii equation (GPE), which is a widespread model that describes the averaged dynamics of the condensate [23, 6]. Depending on the physical situation that one considers, several terms in the right-hand side of the NLSE appear. Our goal is to develop numerical methods for the time integration of a fairly general NLSE including realistic physical situations, that have high-order in time and have good qualitative properties (preservation of mass, energy, etc) over finite times. We develop our analysis in spatial dimension $d\in\{1,2,3\}$ because it fits the physical framework, even if most methods and results naturally extend to higher dimensions. The space variable $x$ will sometimes lie in ${\mathbb{R}}^{d}$, sometimes on the $d$-dimensional torus ${\mathbb{T}}_{\delta}^{d}=({\mathbb{R}}/(\delta{\mathbb{Z}}))^{d}$ (for some $\delta>0$). In this paper, we consider a NLSE of the form

[TABLE]

where $\varphi$ is an unknown function from ${\mathbb{R}}\times{\mathbb{R}}^{d}$ or ${\mathbb{R}}\times{\mathbb{T}}_{\delta}^{d}$ to ${\mathbb{C}}$, $\Delta$ is the Laplace operator, $V$ is some real-valued potential function, $\beta\in{\mathbb{R}}$ is a parameter that measures the local nonlinearity strength, $\lambda\in{\mathbb{R}}$ is a parameter that measures the nonlocal nonlinearity strength with convolution kernel $U$, $\Omega\in{\mathbb{R}}^{d}$ is a vector encoding the direction and the speed of a rotation, and $R$ is a rotation operator that is important in the modeling of rotating BEC (for example, $R=x\wedge(-i\nabla)$ when $d=3$). The NLSE is supplemented with an initial datum $\varphi_{\text{in}}$. The results presented in this paper extend to more general power law nonlinearities such as

[TABLE]

but we restrict ourselves to $K=1$ for the sake of simplicity. Equation (1) is hamiltonian for the energy functional

[TABLE]

provided $U$ is a real-valued convolution kernel, symmetric with respect to the origin 111For real-valued functions, symmetry with respect to the origin is equivalent to real-valued Fourier transform.. In practice the convolution kernel $U$ may for example correspond to a Poisson equation ($U(x)=1/(4\pi|x|)$ in dimension $d=3$) or it may represent dipole-dipole interactions (see [9]).

The main goal of this paper is the analysis of numerical methods for the time integration of (1) that preserve the energy (2) or a discretized analogue of it. In particular, we are interested in the order (in time) of such methods. A well known method is the Crank-Nicolson method introduced in [14] for parabolic problems (see for example a posteriori error estimates in [3]) and applied in [17] to Schrödinger equations. For all nonlinearities, these methods are fully implicit. However, they have second order in time (see [24] for the case of the cubic NLS equation and [28] for the case of a system with possibly fractional derivatives) and preserve discrete analogues of the energy (2) as well as the total mass (squared $L^{2}$-norm) of the solutions. Unfortunately, since they are fully implicit, these methods are costly. To work around this problem, the methods introduced and analyzed in this paper belong to the family of relaxation methods.

Relaxation methods for Schrödinger equations were introduced in [10, 11]. They have been applied to different NLS equation for example in the context of plasma physics [22]. For cubic nonlinearities ($\sigma=1$ in (1)), they are linearly implicit hence very popular [4, 15, 6, 18, 20]. They preserve the $L^{2}$-norm and a discrete analogue of (2). It is well-known that they have numerical order 2 but up to our knowledge there is no proof of order 2 in the literature. This paper presents two new results with respect to relaxation methods for (1). First, we prove rigorously that the classical relaxation method, applied to the classical cubic NLS equation (i.e. (1) with $V\equiv 0$, $\beta=1$, $\sigma=1$, $\lambda=0$ and $\Omega=0$) is of order 2. Second, we present a generalized relaxation method that allows to deal with general power law nonlinearities ($\sigma\neq 1$ in (1)) and with the full GPE. The generalized relaxation methods that we introduce in this context are implicit (actually explicit for $\sigma\leq 4$), have numerical order 2 and we show that they preserve an energy which is also a discretized analogue of (2).

This paper is organized as follows. In Section 2 we recall the definition of the Crank-Nicolson method and give a short explanation of its energy preserving property. Section 3 is devoted to relaxation methods applied to (1): In a first part we recall the method introduced in [10] for the cubic Schrödinger equation and in a second part we give a proof of the optimal order of convergence for an initial datum belonging to $H^{s+4}({\mathbb{R}}^{d})$, $d$ in $\{1,2,3\}$ and $s>d/2$. In section 4 we propose a generalized relaxation method that allows to deal with general nonlinearites and we prove that this method is also an energy preserving method. Section 5 deals with numerical results in different physical models showing the efficiency of the methods.

2. Preservation of energy and Crank-Nicolson scheme

A usual way to prove the conservation of the energy (2) consists in multiplying the equation (1) by $\overline{\partial_{t}\varphi}(t,x)$, where $\bar{z}$ denotes the conjugate value of a complex $z$, integrating over ${\mathbb{R}}^{d}$ and taking the real part of the result. This computation relies on the identity which holds for all smooth functions $\varphi$ of time with values into a space of sufficiently integrable functions

[TABLE]

A possible way to derive numerical schemes that preserve an energy functional is therefore to mimic this identity at the discrete level.

In 1981, Delfour, Fortin and Payre [17], following an idea of Strauss and Vasquez [26], proposed a way to deal with the nonlinear term $|\varphi|^{2\sigma}\varphi$ for the Crank-Nicolson scheme. This method generalizes the second order mid-point scheme for the linear Schrödinger equation

[TABLE]

where $\varphi_{n}(x)$ denotes an approximation of $\varphi(t_{n},x)$ with the discrete time $t_{n}=n\delta t$ defined with the time step $\delta t$. Their approach can be explained as follows. If one looks for a real-valued function $g:{\mathbb{C}}^{2}\rightarrow{\mathbb{R}}$ such that the scheme takes the form

[TABLE]

then, multiplying this relation by $i\overline{\varphi_{n+1}-\varphi_{n}}$, integrating overs ${\mathbb{R}}^{d}$ and taking the real part, as we did in the time-continuous setting above, yields to [math] in the left-hand side and several terms in the right-hand side. Amongst these terms, those involving $g$ are equal to

[TABLE]

since $g$ is real-valued. Let us denote by $G$ the function $v\mapsto|v|^{2\sigma+2}/(2\sigma+2)$. A sufficient condition for the method (4) to preserve an energy of the form (2) is therefore to have

[TABLE]

This is exactly the definition of $g$ chosen in [17].

In the following, the Crank-Nicolson method for the GPE (1) is therefore defined using the formula

[TABLE]

We shall use the notation

[TABLE]

for the Crank-Nicolson method (LABEL:eq:CNmethod). In the expression above, the term corresponding to the nonlinearity should be understood as

[TABLE]

so that it is indeed non-singular and it is consistent with the non-linear term $\beta|\varphi|^{2\sigma}$. The Crank-Nicolson method is fully implicit. It is known to have order two for the cubic NLS equation [24]. Moreover it preserves exactly the $L^{2}$-norm of the solution as well as the following energy:

[TABLE]

with $E$ defined by (2).

In Section 4, we shall use similar ideas to derive energy-preserving relaxation methods for general power laws nonlinearities. Before doing so, we first deal with the classical relaxation method in Section 3.

3. The classical relaxation method

3.1. An energy preserving method

In [11], Besse introduced the usual relaxation method (10) applied to the nonlinear Schrödinger equation (1) with $V=0$, $\lambda=0$, $\Omega=0$ and $\sigma=1$ that is known as the cubic nonlinear Schrödigner equation

[TABLE]

with $\varphi(0,x)=\varphi_{\text{in}}(x)$. The idea of the relaxation method is to add to (8) a new unknown $\Upsilon=|\varphi|^{2}$ and the equation (8) is transformed in

[TABLE]

The relaxation method then consists in discretizing both equations respectively at discrete times $t_{n}$ and $t_{n+1/2}$ and to solve iteratively

[TABLE]

to compute approximations $\varphi_{n}$ of $\varphi(n\delta t)$. This system is usually initialized with $\Upsilon_{-1/2}=|\varphi(-\delta t/2)|^{2}$ or by second order approximation of $|\varphi(-\delta t/2)|^{2}$. This method is linearly implicit (recall that $\sigma=1$). Moreover, it is known to preserve exactly the $L^{2}$-norm and the discrete energy [11]:

[TABLE]

Indeed

[TABLE]

But

[TABLE]

Using the definition of $\Upsilon_{\cdot+1/2}$ in (10), a simple computation leads to

[TABLE]

We therefore conclude that

[TABLE]

which allows to prove the conservation of the discrete energy (11).

It is interesting to note the consistency of the energy associated to relaxation scheme with the energy (2) for cubic nonlinear Schrödinger equation

[TABLE]

The relaxation method was proved to converge in [11] but consistency analysis was missing. We present it in the next subsection.

3.2. Consistency analysis for NLS equation with cubic nonlinearity

The aim of this subsection is to prove that this method has temporal order 2 under fairly general assumptions. This fact is supported by numerical evidences in the literature for years. We provide the first rigorous proof below.

The first equation in (10) is the discrete equivalent of the continuous constraint $\Upsilon=|\varphi|^{2}$. In particular, this constraint is not an evolution equation. Therefore, we use the ideas introduced in [11] and rewrite the continuous equation (1) (recall that $V=0$, $\lambda=0$ and $\Omega=0$ and $\sigma=1$) as the system

[TABLE]

The discrete system (10) has a discrete augmented equivalent (see [10, 11]). Let us denote by $v_{n+\frac{1}{2}}=\dfrac{\varphi_{n+1}-\varphi_{n}}{\delta t}$ the discrete time derivative of $\varphi_{n}$ and define the nonlinearities as

[TABLE]

The augmented system writes

[TABLE]

This system allows to compute $(\varphi_{n+2},\Upsilon_{n+\frac{3}{2}},v_{n+\frac{3}{2}})$ from

[TABLE]

We consider the system (14) as the mapping

[TABLE]

The seven variables involved in $X_{n}$ are not independent. If they satisfiy the five relations

[TABLE]

then the seven variables in $X_{n+1}$ satisfy the same five relations with $n$ replaced by $n+1$. This fact is proved in [11]. We describe now how to build the seven initial data in $X_{0}$ from $\varphi_{0}=\varphi_{\text{in}}$ and $\Upsilon_{-\frac{1}{2}}$ so that they satisfy the five relations above:

[TABLE]

Let us define the operators

[TABLE]

and the matrix of operators

[TABLE]

The mapping $X_{n}\mapsto X_{n+1}$ reads

[TABLE]

In a more compact form, we define $\mathcal{B}$ and $\mathcal{M}$ so that the mapping (16) reads

[TABLE]

We introduce the hypotheses that will allow us to prove our consistency and convergence result for the classical relaxation method (10) in Theorem 8.

Remark 1**.**

The results below extend to more general cases. In particular, one may treat the case where $V$ is non zero smooth autonomous potential such that the multiplication by $V$ is a bounded linear operator between Sobolev spaces and the case where $V$ is a nonautonomous such operator with sufficient regularity with respect to time.

Hypotheses 2**.**

We fix $d\in\{1,2,3\}$ and $s>d/2$. We assume $\varphi_{0}\in H^{s+4}({\mathbb{R}}^{d})$ is given. We denote by $T^{*}>0$ the existence time of the maximal solution $\varphi$ of the Cauchy problem (1) (with $V=0$, $\lambda=0$, $\Omega=0$ and $\sigma=1$) in $H^{s+4}({\mathbb{R}}^{d})$. We assume there exists $\delta t_{0}>0$ such that $\tau\mapsto\varphi(\tau,\cdot)$ is a smooth map from $(-\delta t_{0},T^{*})$ to $H^{s+4}({\mathbb{R}}^{d})$. Moreover, we assume that there exists $R_{1},R_{2}>0$ such that for all $\Upsilon_{-1/2}\in H^{s+4}({\mathbb{R}}^{d})$ with $\|\Upsilon_{-1/2}\|_{H^{s+4}}\leq R_{1}$, the numerical solution $X_{n}$ (with initial datum (15)) is uniquely determined by (16) for all $\delta t\in(0,\delta t_{0})$ and all $n\in{\mathbb{N}}$ such that $nh\leq T$, and it satisfies $\|X_{n}\|_{(H^{s+2}({\mathbb{R}}^{d}))^{7}}\leq R_{2}$ for all such $n$.

Remark 3**.**

The hypotheses above on the exact solution are fullfilled in several cases. For example, for the exact solution $\varphi$, it is well known (see [19]) that $T^{*}=+\infty$ in at least two cases:

• •

if $\varphi_{\text{in}}$ has small $H^{s+4}$-norm and $\beta<0$

• •

if $\beta>0$ and $\varphi_{\text{in}}\in H^{s+4}({\mathbb{R}}^{d})$.

For the numerical solution, they are fullfilled provided $T^{*}<+\infty$ (see [11]) and also when $T^{*}=+\infty$ and $\beta>0$ (see [10]).

Let $t_{n}=n\delta t$ denote the discrete times and $t\mapsto X(t)$ the vector

[TABLE]

Using the definition of $\mathcal{M}$ (see (13) and (17)), the fact that $H^{s}({\mathbb{R}}^{d})$, is an algebra since $s>d/2$, and the fact that the exact and numerical solutions stay in a bounded set of $H^{s}({\mathbb{R}}^{d})$, it is easy to prove the following lemma.

Lemma 4**.**

Assume Hypotheses 2 is satisfied. There exists $C>0$ such that for all $\delta t\in(0,\delta t_{0})$, all $n\in\mathbb{N}$ such that $(n+1)\delta t\leq T$

[TABLE]

Note that the constant $C$ depends only on the initial data.

Before starting the proof of our main result of this section (see Theorem 8), we state and prove another lemma.

Lemma 5**.**

There exists a constant $b>0$ such that the operator $\mathcal{B}$ defined in (16) and (17) satisfies for all $\delta t>0$, and $n\in{\mathbb{N}}$,

[TABLE]

where $\vvvert\cdot\vvvert$ is the norm of linear continuous operators from $(H^{s}({\mathbb{R}}^{d}))^{7}$ to itself.

Proof.

The operator $\mathcal{B}$ is defined by three diagonal blocks

[TABLE]

The first block $\mathcal{B}_{1}$ is an isometry from $(H^{s}({\mathbb{R}}^{d}))^{2}$ to itself and so are all its powers. The powers of the second block $\mathcal{B}_{2}$ read for all $n\geq 2$

[TABLE]

Since all the powers of $A$ are of norm less than one, we infer that the norm of $\mathcal{B}_{2}^{n}$ is less than $\sqrt{3}$. Since the norm of $B$ is less than one, we infer that the norm of

[TABLE]

is less than $\sqrt{3}$.

The block $\mathcal{B}_{3}$ can be diagonalized by blocks as

[TABLE]

so that for all $n\in{\mathbb{N}}$,

[TABLE]

Therefore, the last $2\times 2$ block of $\mathcal{B}^{n}\mathcal{C}$ reads

[TABLE]

Since $B(I+A)^{-1}=(1/2)I$, we infer

[TABLE]

Since all the powers of $A$ are of norm less than $1$, we infer that

[TABLE]

This proves the result. ∎

Remark 6**.**

The powers of the operator $\mathcal{B}$ are not uniformly bounded. However, the powers of $\mathcal{B}$ mutiplied by $\mathcal{C}$ are uniformly bounded as shown above. The main reason is that the third matrix in the right hand side of (20) becomes singular at the end of the spectrum of $\Delta$.

Lemma 7**.**

Assume $d,s,R_{1},\delta t_{0}$ and $\varphi_{\text{in}}$ are given as in Hypotheses 2. There exists $c>0$ such that for all $\Upsilon_{-1/2}\in H^{s+4}$ with $\|\Upsilon_{-1/2}\|_{H^{s+4}}\leq R_{1}$ and all $\delta t\in(0,\delta t_{0})$,

[TABLE]

Proof.

The $H^{s+2}$-norm of each of the seven components of the vector $X_{0}-X(0)$ is estimated separately. The $H^{s+2}$-norm of the first component $\Upsilon_{-1/2}-|\varphi(-\delta t/2)|^{2}$ is bounded by the $H^{s+4}$-norm of the same quantity. For the $H^{s+2}$-norm of the second component, we define $f(t)=|\varphi(t)|^{2}$, which is a smooth function from $(-\delta t_{0},\delta t_{0})$ to $H^{s+4}$, thanks to Hypotheses 2. We may write using a Taylor formula at [math]

[TABLE]

and similarly

[TABLE]

Therefore, one has

[TABLE]

By triangle inequality, We infer that

[TABLE]

where $c=\max(1,\sup_{\sigma\in(-\delta t_{0}),\delta t_{0}}\|f^{\prime\prime}(\sigma)\|_{H^{s+2}})$ only depends on the exact solution of (1). We now estimate the $H^{s+2}$-norm of the fifth component of $X_{0}-X(0)$. Note that the $H^{s+2}$-norm of the third component can be estimated the very same way and that the $H^{s+2}$-norm of the fourth component is zero. We start with the identity

[TABLE]

and we denote by $r(\delta t)$ the consistency error defined by

[TABLE]

Using the fact that $t\mapsto\varphi(t,\cdot)$ is a smooth function from $(-\delta t_{0},\delta t_{0})$ to $H^{s+4}$, we may write another Taylor expansion to obtain

[TABLE]

Using

[TABLE]

we obtain

[TABLE]

Note that we have

[TABLE]

[TABLE]

and

[TABLE]

where $c$ doesn’t depend on $\delta t$. Then we infer from the estimates above that

[TABLE]

Now, we denote by $e(\delta t)$ the fifth component $\varphi_{1}-\varphi(\delta t)$ of the vector $X_{0}-X(0)$ and we have

[TABLE]

We want to estimate the $H^{s+2}$-norm of $e(\delta t)$ using this relation and the estimate (23). To this aim, we take $\alpha\in{\mathbb{N}}^{d}$ with $|\alpha|=\alpha_{1}+\dots+\alpha_{d}\leq s+2$ and differentiate the relation above to obtain that

[TABLE]

Multiplying this relation by $\overline{\partial_{x}^{\alpha}e(\delta t)}$ integrating over ${\mathbb{R}}^{d}$, and taking the real part, we obtain

[TABLE]

When $\alpha=0_{{\mathbb{N}}^{d}}$, the first term on the right hand side in the equation above vanishes and

[TABLE]

using Cauchy-Schwartz inequality. We infer,

[TABLE]

Now, when $\alpha\neq 0_{{\mathbb{N}}^{d}}$, the Leibniz’ rule applied in (24) provides us with:

[TABLE]

where $c(k,\alpha)$ are integers. Note that since $\Upsilon_{1/2}$ is real valued, the term corresponding to $k=0_{{\mathbb{N}}^{d}}$ in the sum vanishes. Then, with Cauchy-Schwarz inequality, we infer

[TABLE]

Since $k\leq\alpha$ in the sum above, we have $|k|\leq|\alpha|$. For such $k$, one has

[TABLE]

where $c$ is the Sobolev constant of the injection from $H^{(d+1)/2}({\mathbb{R}}^{d})$ to $L^{\infty}({\mathbb{R}}^{d})$ and where we have used the fact that $d\in\{1,2,3\}$. Since the $H^{s+4}$-norm of $\Upsilon_{1/2}$ is controlled by (22), we have

[TABLE]

where $c_{\alpha}$ does not depend on $\delta t\in(0,\delta t_{0})$. Then, by induction on $|\alpha|\in\{0,\dots,s+2\}$ starting with (25) for $\alpha=0_{{\mathbb{N}}^{d}}$, we have for some positive constant $c_{s+2}$:

[TABLE]

for all $\delta t\in(0,\delta t_{0})$ and all $\alpha\in{\mathbb{N}}^{d}$ with $|\alpha|\leq s+2$. Therefore, we have $\|e(\delta t)\|_{H^{s+2}}$ is controlled by $\|r(\delta t)\|_{H^{s+2}}$ and the conclusion for the fifth term follows using (23):

[TABLE]

where $c$ does not depend on $\delta t\in(0,\delta t_{0})$.

It remains to estimate the $H^{s+2}$-norm of the sixth and seventh components of the vector $X_{0}-X(0)$. We only give the details for the seventh term since the computation is similar and even simpler for the sixth term. Let us denote by $p(\delta t)$ the consistency error defined as

[TABLE]

Using a Taylor expansion, since the exact solution is a smooth function from $(-\delta t_{0},\delta t_{0})$ to $H^{s+2}$, we have

[TABLE]

Let us denote by $q(\delta t)=v_{1/2}-\partial_{t}\varphi\left(\dfrac{\delta t}{2}\right)$ the seventh component of $X_{0}-X(0)$. We have

[TABLE]

Using estimates (26) and (27) we have by triangle inequality,

[TABLE]

This concludes the proof of the lemma. ∎

We prove below that the relaxation method (10) is of order 2.

Theorem 8**.**

Assume $d$, $s$, $R_{1}>0$, $\varphi_{\text{in}}\in H^{s+4}({\mathbb{R}}^{d})$, $T<T^{\star}$ are given and satisfy Hypotheses 2. There exists $\mathscr{C}>0$ and $\delta t_{0}>0$ (smaller than the one in the hypotheses) such that for all $\Upsilon_{-1/2}\in H^{s+4}({\mathbb{R}}^{d})$ with $\|\Upsilon_{-1/2}\|_{H^{s+4}({\mathbb{R}}^{d})}\leq R_{1}$, all $n\in\mathbb{N}$ and all $\delta t\in(0,\delta t_{0})$ with $n\delta t\leq T$,

[TABLE]

Proof.

First, the initial datum $X_{0}$ is computed from $\varphi_{0}=\varphi_{\text{in}}$ and $\Upsilon_{-1/2}$ using (15). Therefore, using Lemma 7, there exists a constant $c>0$ such that for all $\delta t\in(0,\delta t_{0})$ and all $\|\Upsilon_{-1/2}\|_{H^{s+4}({\mathbb{R}}^{d})}\leq R_{1}$, we have the estimate (21)

[TABLE]

Second, we define the consistency error of the relaxation scheme at time $t_{n}=n\delta t\leq T$ by the formula

[TABLE]

Substracting this definition from (17), we obtain

[TABLE]

From now on, we set for all $n\in{\mathbb{N}}$ such that $n\delta t\leq T$, $e_{n}=X_{n}-X(t_{n})$. Iterating the relation above, we obtain, as long as $n\delta t\leq T$,

[TABLE]

This implies

[TABLE]

using Lemmas 4 and 5. Then

[TABLE]

Let $\delta t$ be small enough to ensure that $1-C\delta tb\geq\dfrac{1}{2}$. This implies

[TABLE]

Taylor expansions and the fact that one can differentiate the last two lines of (12) with respect to time show that there exists a constant $Q$ such that for all $n$ and $\delta t$ with $n\delta t\leq T$

[TABLE]

This implies that $2b\sum_{k=0}^{n-1}\|R_{k}(\delta t)\|_{(H^{s+2}({\mathbb{R}}^{d}))^{7}}\leq 2bQT\delta t^{2}$. Then we obtain

[TABLE]

Using a discrete Gronwall Lemma (see section 5 in [21]), we get

[TABLE]

This estimate and (21) prove the result. ∎

4. The generalized relaxation method for general nonlinearities

4.1. Generalized relaxation method for NLS equation

We start by considering the simplified equation (1) with zero potential ($V=0$), no convolution operator ($U=0$) and without rotation ($\Omega=0$). Assuming $\sigma\in{\mathbb{N}}^{\star}$, we are therefore dealing with the classical nonlinear Schrödinger equation

[TABLE]

with initial datum $\varphi_{\text{in}}$.

The original relaxation method applied to (30) would consists in adding the variable $\Upsilon$ to (30) and discretizing the following continuous system as discrete times $t_{n}$ and $t_{n+1/2}$

[TABLE]

It is however known to not conserve energy functional.

As a generalization of the classical relaxation method (10), which is designed for the special cubic case ($\sigma=1$), we propose the following method, which allows for a general nonlinearity exponent $\sigma\in{\mathbb{N}}^{\star}$. We propose to substitute system (31) by

[TABLE]

The modification seems ligth but allows to build an energy preserving scheme. The second equation is approximated to second order at time $t_{n+1/2}$ by

[TABLE]

We now want to find an approximation of $\gamma^{\sigma}=|\varphi|^{2\sigma}=\gamma^{\sigma-1}|\varphi|^{2}$ that allow energy conservation following the ideas that were presented in section 2. As for the classical relaxation method, we note that

[TABLE]

The last term also reads

[TABLE]

The only choice that allow to preserve energy is to choose

[TABLE]

Moreover, we remark that

[TABLE]

is a second order approximation of $\gamma^{\sigma}=\gamma^{\sigma-1}|\varphi|^{2}$ at time $t=t_{n}$.

This method is therefore designed so that it preserves exactly the following energy

[TABLE]

Since at continuous level $\gamma^{\sigma}=|\varphi|^{2\sigma}$, $E_{\mathrm{rlx}}(\varphi,\gamma)$ reduces to the true energy

[TABLE]

Moreover, the generalized relaxation method preserves the $L^{2}$-norm of the solution. We present in Theorem 9 a more general method, which includes possibly non zero other terms in the equation (see Section 4.2).

Starting with $\varphi_{0}=\varphi_{\text{in}}$ and $\gamma_{-1/2}$ approximating $|\varphi(-\delta t/2)|^{2}$, the generalized relaxation method is

[TABLE]

Like for the Crank-Nicolson scheme and equation (6), the first equation also reads

[TABLE]

so we also have the other version of generalized relaxation method

[TABLE]

Note that, when $\sigma=1$, the generalized relaxation method (34) reduces to the classical relaxation method (10). In contrast to the classical relaxation method (10), when $\sigma\geq 2$, the generalized relaxation method (34) is fully implicit on its first stage, and linearly implicit in its second stage. For small values of $\sigma$ ($\sigma=1,2,3,4$) the first stage of (34) is polynomial of degree $\sigma$ and hence we can use explicit formulas for the computation of the solution $\gamma_{n+1/2}$. For example, for quintic nonlinearity and $\sigma=2$, we have explicit solutions to the quadratic equation

[TABLE]

where $k_{1}=(3/2)|\varphi_{n}|^{2}-\gamma_{n-1/2}$.

For higher values of $\sigma$, one can use the following iterative fixed-point procedure, starting with $\gamma_{n+1/2,0}=|\varphi_{n}|^{2}$ :

[TABLE]

which one stops when $\|\gamma_{n+1/2,p+1}-\gamma_{n+1/2,p}\|_{L^{2}}$ is below some small tolerance parameter, and for this index $p$, one sets $\gamma_{n+1/2}=\gamma_{n+1/2,p+1}$.

Numerically, this generalized relaxation method has order 2, as we will see in the numerical experiments section 5. However, we do not address this theoretical question in this paper. In the next subsection, we show how one can design a generalized relaxation method similar to (34) in order to treat the cases with non-zero potential $V$, non-zero convolution operator $U$ or non-zero rotation $\Omega$.

4.2. A generalized relaxation method for GPE

We propose the following generalization of the relaxation method introduced in [11], which uses two additional unknowns $\gamma$ and $\Upsilon$: starting from $\varphi_{0}=\varphi_{\text{in}}$, we initialize $\Upsilon_{-1/2}$ and $\gamma_{-1/2}$ with approximations of $|\varphi(-\delta t/2)|^{2}$ and compute for $n\in{\mathbb{N}}$, $(\varphi_{n+1},\Upsilon_{n+1/2},\gamma_{n+1/2})$ from $(\varphi_{n},\Upsilon_{n-1/2},\gamma_{n-1/2})$ using the formulae

[TABLE]

The first equation of (36) is implicit and local. In order to solve it, we propose two different approaches. The first one, for small integer values of $\sigma$, consists in using exact formulae, since the equation is polynomial of low degree in these cases. The second approach, is to use a fixed point method, as described in (35). The second equation is explicit, as in the classical relaxation method (10). The third equation of (36) is implicit and non-local. We propose again two different approaches to solve it. The first approach, presented here for $x\in{\mathbb{T}}_{\delta}^{d}$, consists in using another fixed point iteration method: one starts from $\varphi_{n+1}^{0}=\varphi_{n}$ and computes $\varphi_{n+1}^{p+1}$ from $\varphi_{n+1}^{p}$ using the iterative procedure

[TABLE]

where $\mathcal{F}$ stands for the Fourier transform in space as defined in Appendix A, and we set $\varphi_{n+1}=\lim_{p\to\infty}\varphi_{n+1}^{p}$. In practice, the iterative procedure $\varphi_{n+1}^{p}\rightarrow\varphi_{n+1}^{p+1}$ stops whenever the $L^{2}$-norm of the difference between two consecutive steps is below some small tolerance parameter. Note that we decided to implicit the Laplace operator in the iterative procedure (37) in order to ensure that the Sobolev regularity of $\varphi_{n+1}^{p}$ does no decrease a priori when $p$ increases. Alternatively, the other approach to solve the last equation of (36) consists in following [5], and using the linearity of the equation. Let us introduce the new unknown $\varphi_{n+1/2}=(\varphi_{n+1}+\varphi_{n})/2$, so that the equation reads

[TABLE]

This equation can also be written as

[TABLE]

where

[TABLE]

and

[TABLE]

The precondionning operator $P$ can be easily inverted in Fourier space and the solution of (38) can be obtained by a Krylov method. Note that other choices of preconditionning operator are possible (see [5]).

In the following, we shall use the notation

[TABLE]

for the method (36) above.

4.2.1. Energy preservation property for generalized relaxation methods

The generalized relaxation method (36) is designed to preserve exactly the following energy:

[TABLE]

as we prove in Theorem 9. Note that this energy is consistent with the energy (2) of equation (1) in the sense that

[TABLE]

Theorem 9**.**

The generalized relaxation method (36) applied to the equation (1) with initial data $\varphi_{0}=\varphi_{\text{in}}$, $\gamma_{-1/2}$ and $\Upsilon_{-1/2}$ preserves exactly the $L^{2}$ norm and the energy functional $E_{\rm rlx}$ defined in (4.2.1) in the sense that for all $n\in{\mathbb{N}}$ such that a solution of (36) is defined, we have

[TABLE]

Proof.

Multiplying the last equation of (36) by $\overline{\varphi_{n+1}-\varphi_{n}}$, integrating in space, and taking the real part, one finds that a sum of 5 terms is equal to

[TABLE]

The first term reads

[TABLE]

The second term reads

[TABLE]

Using the first equation of (36), the third term reads

[TABLE]

Using the fact that the convolution kernel $U$ is real-valued and symmetric with respect to the origin, and also using the second equation of (36), the fourth term reads

[TABLE]

Finally, thanks to the fact that the operator $R$ is symmetric, the last term reads

[TABLE]

The proof of the preservation of the energy (40) is completed by adding altogether (42), (43), (44), (44), and (46), and using relation (41). The preservation of the $L^{2}$-norm follows from multiplying the last relation in (36) by $\overline{\varphi_{n+1}+\varphi_{n}}$, integrating in space, and taking the imaginary part. ∎

5. Numerical experiments

In this section, we make some numerical experiments and show that the classical and generalized relaxation methods are efficient methods that preserve mass and energy to machine epsilon.

5.1. One dimensional example: the one-dimensional quintic and septic NLS equation

We present in this subsection some numerical experiments to show the efficiency of the generalized relaxation (34) compared with the Crank-Nicolson scheme (LABEL:eq:CNmethod) when considering

[TABLE]

with $\varphi(0,x)=\varphi_{\mathrm{in}}(x)$ and $(t,x)\in[0,T]\times(x_{\ell},x_{r})$. To deal with the space variable, we discretize the space operators using Fourier spectral approximation and consider periodic boundary conditions. The spatial mesh size is $\delta x>0$ with $\delta x=(x_{r}-x_{\ell})/J$ with $J=2^{P}$, $P\in{\mathbb{N}}^{*}$. The time step is $\delta t=T/N$ for some $N\in{\mathbb{N}}^{\star}$. The grid points and the discrete times are

[TABLE]

Let $\varphi_{j,n}$ be the approximation of $\varphi(t_{n},x_{j})$ satisfying

[TABLE]

where $\hat{\varphi}_{k,n}$ denotes the discrete Fourier transform of sequence $(\varphi_{j,n})_{j}$ given by

[TABLE]

where $\omega_{J}=\exp{(-2i\pi/J)}$. We also apply the discrete Fourier transform to the approximation $\gamma_{j,{n+1/2}}$ of $|\varphi(t_{n+1/2},x_{j})|^{2}$. Let us define the discrete gradient operator $\nabla_{d}$

[TABLE]

Let us denote by $\Pi_{d}$ the projection operator

[TABLE]

We define the discrete $\ell^{r}$ norm on ${\mathbb{C}}^{J}$ as

[TABLE]

the mean

[TABLE]

and the discrete energies:

[TABLE]

and

[TABLE]

As in the continuous case, we have $E_{d}(v)=E_{\mathrm{rlx},d}(v,v)$ for any $v\in{\mathbb{C}}^{J}$.

Using these definitions, the energy conservation for the relaxation scheme is built through the following relative error

[TABLE]

For the Crank-Nicolson scheme, it is defined by

[TABLE]

We present on Figure 1 the evolution of $\mathcal{E}_{E,\delta t}$ for various $\delta t$ when $\sigma=2$ and $\sigma=3$ for both classical and generalized relaxation scheme compared to Crank-Nicolson scheme. The initial datum is chosen to be $\varphi_{\mathrm{in}}(x)=\exp(-x^{2})$ and $\beta=-1$. The time-space domain is $[0,1/2]\times[-30,30]$. We consider periodic boundary conditions and the interval $[-30,30]$ is meshed with $2^{13}+1$ nodes. As it is known, the standard relaxation scheme does not preserve the energy when $\sigma=2$ or $\sigma=3$ and the error curve show a second order convergence. On the other hand, both generalized relaxation and Crank-Nicolson preserve energy to epsilon machine.

5.2. A nonlocal Schrödinger equation with cubic-quintic nonlinearities

In [13], Chen et al. investigate the interactions of dark solitons under competing nonlocal cubic and local quintic nonlinearities. They consider a 1D optical beam with an amplitude $\varphi(t,x)$ with competing nonlocal cubic and local quintic nonlinearities. Their model is given by the following nonlocal nonlinear Schrödinger equation

[TABLE]

where $U(x)$ is the nonlocal kernel and $\alpha_{1}$ and $\alpha_{2}$ are real parameters. When $\alpha_{j}>0$, $j=1,2$, we are considering focusing nonlinearities, whereas when $\alpha_{j}<0$, the nonlinearities are defocusing. This model is generalized in two dimensions in [25] for the study of vortex solitons where $\partial_{x}^{2}$ is replaced by the two-dimensional Laplace operator. Typically, the kernel is regular and is given by

[TABLE]

or

[TABLE]

the parameter $\mu$ allowing to control the width of the kernel, or in other words, the nonlocality strength. At the limit $\mu\to 0$, the kernel $U_{j}$, $j=1,2$, tends to a Dirac distribution and we recover a local nonlinear model. The equation (50) is associated to the energy

[TABLE]

We reproduce the numerical experiments presented in [13, 25] and show the ability of the generalized relaxation method to preserve the energy (4.2.1) in the form (56) below. The generalized relaxation method for (50) consists in approximating the system of equations

[TABLE]

and the numerical scheme reads

[TABLE]

with $\varphi_{0}(x)=\varphi_{\mathrm{in}}(x))$ and $\Upsilon_{-1/2}(x)=\gamma_{-1/2}(x)$ is some second order approximation of $\varphi(-\delta t/2,x)$. In our numerical experiments, this approximation is obtained by applying the Crank-Nicolson scheme starting from $\varphi_{0}$ on reverse time step $-\delta t/2$. The energy associated to (54) is

[TABLE]

and we have the conservation property (40).

We are first interested in the one-dimensional case with kernel $U_{1}$ and we choose a defocusing nonlocal nonlinearity by considering $\alpha_{1}=-1$. Like in the previous subsection, the space variable is discretized using Fourier spectral approximation and we consider periodic boundary conditions. In order to avoid any interaction between the nonlocal kernel and the boundaries, we take a very large domain, typically $x\in[-256\pi,256\pi]$ discretized with $J=2^{14}+1$ nodes. The time step is $\delta t=5\cdot 10^{-3}$ and the final time is $T=30$. The initial datum is made of two solitons at a relative distance of $2x_{0}$, where $x_{0}=1$. Following [13], we choose

[TABLE]

where $D$ is the positive root of the equation

[TABLE]

We present results in Figure 2 and 3 for a defocusing cubic nonlinearity $\alpha_{2}=-0.5$ and a focusing one where $\alpha_{2}=0.1$. Two values for $\mu$ are proposed $\mu=0.5$ and $\mu=2.5$. The behavior of the defocusing-defocusing case is surprising for “strong” nonlocality $\mu=2.5$ since eventually the two solitons breathe.

The evolution of the relative energy error

[TABLE]

is presented on Figure 4. The energy is clearly very well preserved.

Concerning the two dimensionsal case, we consider the kernel $U_{2}$, $\alpha_{1}=1$ and $\alpha_{2}=-0.02$. The initial datum is chosen as a vortex beam with angular momentum

[TABLE]

where $r=\sqrt{x_{1}^{2}+x_{2}^{2}}$, $m=1$ is the topological charge, $\phi$ is such that $x_{1}+ix_{2}=r\exp(i\phi)$ and $A=5.8$ is the amplitude. The computational domain is $[-8,8]^{2}$ and we use $256$ Fourier modes. The final time of simulation is $T=10$ and the time step is $5\cdot 10^{-3}$. The width of the Gaussian kernel $U_{2}$ is determined by $\mu=0.4$. Once again, the energy is very well preserved, since the relative error energy is bounded by $10^{-11}$ at the end of the simulation (see Figure 5). We present in Figure 6 the evolution of the solution with respect to time. The 3D representation of the time evolution is presented in figure 6a. A 2D projection of it is displayed in Fig. 6b and the final solution in Fig. 6c. It is interesting to note that the same kind of breathing behaviour observed in the one dimensional setting is also present in 2D experiment.

5.3. Two dimensional rotating dipolar Bose-Einstein condensate

In this final subsection, we present results of simulations with the relaxation scheme for a rotating Bose-Einstein condensate subject to long range dipole-dipole interaction (DDI). The model is (see Eq. (8.15)-(8.17) in [8])

[TABLE]

where $L_{x_{3}}=-i(x_{1}\partial_{x_{2}}-x_{2}\partial_{x_{1}})$ is the $x_{3}$-component of the angular momentum and $\Omega$ represents the rotating frequency. The nonlinear parameters $\beta$ and $\lambda$ respectively describe the strength of the short-range two-body interactions in the condensate and the strength of the dipolar interaction modeled with a Coulomb potential. The real-valued external trapping potential $V$ is chosen as $V(x)=(\gamma_{x_{1}}^{2}{x_{1}}^{2}+\gamma_{x_{2}}^{2}{x_{2}}^{2})/2$, where $\gamma_{x_{1}}>0$, $\gamma_{x_{2}}>0$ are dimensionless constants proportional to the trapping frequencies in the both directions. The normal $\mathbf{n}=(n_{1},n_{2},n_{3})^{T}$ represents the dipole axis and we define $\mathbf{n}_{\perp}=(n_{1},n_{2})^{T}$ and $\partial_{\mathbf{n}_{\perp}}=\mathbf{n}_{\perp}\cdot\nabla$. The energy is given by

[TABLE]

The initial datum is computed as a minizer of the energy on the intersection of the unit sphere of $L^{2}$ with the energy space

[TABLE]

We use a preconditioned nonlinear conjugate gradient method developed in [7]. Contrary to the previous subsection where the convolution kernel was regular, the computation of the nonlocal term $\psi$ with the Coulomb potential is known to be a costly task. We have chosen to apply the technique developed in [27], which allows fast convolution using truncated Green’s functions. The parameters of the simulation are $\Omega=0.97$, $\gamma_{x_{1}}=\gamma_{x_{2}}=1$, $\lambda=175$, $\beta=(250-\lambda)\sqrt{5/\pi}$. The computational domain is $(-16,16)^{2}$ with $2^{8}=256$ Fourier modes in each direction, and the computational time is set to $T=15$ with the time step $\delta t=10^{-3}$. Initially, the dipole axis is $\mathbf{n}=(1,0,0)^{T}$ and we change it to $\mathbf{n}=(\cos(\pi/3),\sin(\pi/3),0)^{t}$. We plot on Figure 7 the evolution of the solution from time $t=0$ to time $t=20$. We add the direction of the dipole axis on each frame. The axis are removed to make the presentation clearer. As it can be seen on Figure 8, the energy is very well preserved as expected.

6. Conclusions and future works

In this paper we have given for the first time a proof of the second order of the relaxation method introduce in [10] for the cubic nonlinear Schrödinger equation. We also have extended the previous method to deal with general power law nonlinearites showing that this method is still an energy preserving method. Note that, since Crank-Nicolson methods are one-step methods, one may use composition techniques to achieve higher orders (4,6,8, etc) while preserving both an energy and the $L^{2}$-norm, and this is not possible so straightforwardly for relaxation methods. Therefore, in future works we want to focus on energy preserving relaxation methods which have higher order.

Appendix A Notation

In this section, the symbol $\mathcal{F}$ will denote different operators depending on whether $x$ belongs to ${\mathbb{R}}^{d}$ or $x$ belongs to ${\mathbb{T}}_{\delta}^{d}$. When $x\in{\mathbb{R}}^{d}$, $\mathcal{F}$ denotes the Fourier transform ${\mathbb{R}}^{d}$ defined for $\varphi\in L^{1}({\mathbb{R}}^{d})$ via the formula

[TABLE]

and its inverse ${\mathcal{F}}^{-1}$ is defined for $\phi\in L^{1}({\mathbb{R}}^{d})$ through the formula

[TABLE]

When $x\in{\mathbb{T}}_{\delta}^{d}=({\mathbb{R}}/(\delta{\mathbb{Z}}))^{d}$, the $\mathcal{F}$ denotes the Fourier transform defined for $\varphi\in L^{1}({\mathbb{T}}_{\delta}^{d})$ via the formula

[TABLE]

and its inverse ${\mathcal{F}}^{-1}$ is defined for $\phi\in\ell^{1}(2\pi\delta^{-1}{\mathbb{Z}}^{d})$ through the formula

[TABLE]