Stochastic symplectic and multi-symplectic methods for stochastic NLS equation with quadratic potential

# Stochastic symplectic and multi-symplectic methods for stochastic NLS equation with quadratic potential

Jialin Hong Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing, China Lijun Miao Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing, China  and  Liying Zhang Department of Mathematics, College of Sciences, China University of Mining and Technology, Beijing 100083, China.
September 23, 2019
###### Abstract.

Stochastic nonlinear Schrödinger equation with quadratic potential models Bose–Einstein condensations under a magnetic trap when , which is not only an infinite-dimensional stochastic Hamiltonian system but also a stochastic Hamiltonian partial differential equation essentially. We firstly indicate that this equation possesses stochastic symplectic structure and stochastic multi-symplectic conservation law. It is also shown that the charge and energy have the evolution laws of linear growth with respect to time in the sense of expectation. Moreover, we propose a stochastic symplectic scheme in the temporal discretization, and give the error estimate of this scheme in probability. In order to simulate the evolution laws of charge and energy numerically, we also give a stochastic multi-symplectic scheme, of which the corresponding discrete charge and energy are in accordance with the continuous case. Numerical experiments are performed to test the proposed numerical scheme.

###### Key words and phrases:
stochastic nonlinear Schrödinger equation, quadratic potential, additive noise, stochastic symplectic scheme, stochastic multi-symplectic scheme
###### 2010 Mathematics Subject Classification:
60H35, 65C30, 65P10

## 1. Introduction

In this paper, we are interested in the numerical methods for the following stochastic nonlinear Schrödinger equation with quadratic potential and additive noise

 (1.1) idu+(Δu+θ|x|2u+λ|u|2σu)dt=dW,(t,x)∈(0,T]×D,u(t,0)=u(t,1)=0,t∈(0,T],u(0,x)=u0(x),x∈D,

where , , , , and denotes a Wiener process. This equation models Bose–Einstein condensations under a magnetic trap when , the quadratic potential describes the magnetic filed whose role is to confine the movements of particles (see e.g., ), and the additive noise expresses the random perturbations (). If the noise vanishes in Eq. (1.1), then it is a deterministic nonlinear Schrödinger equation with quadratic potential, which possesses the charge and energy conservation laws.

If , it reduces to the well-known stochastic nonlinear Schrödinger equation without potential.  studies the well-posedness for it and gets the evolution laws of charge and energy. In order to understand the behavior of solution, numerical schemes are proposed to approximate the solution of this equation and estimate the order of the error, (see e.g.,[2, 4, 6, 10]). Furthermore, this equation is a stochastic Hamiltonian system, many authors focus on the stochastic symplectic and multi-symplectic methods for it. For example,  proposes the stochastic symplectic structure and constructs a stochastic symplectic scheme for it.  establishes the theory about stochastic multi-symplectic conservation law for stochastic Hamiltonian partial differential equations and proposes a stochastic multi-symplectic method. For more results on this subject, we refer the reader to [4, 5, 7, 12, 13, 14].

If , the existence of the unique solution of (1.1) has been established in . To our knowledge, there has been no work concerning structure-preserving numerical schemes for this equation. It is well-known that this kind of schemes are shown to be superior to conventional ones especially in long time computation for Hamiltonian system.This advantage motivates us to consider the symplectic and multi-symplectic structures for (1.1) and design the corresponding numerical schemes to preserve these structures and properties of charge and energy.

In this paper, we first indicate that (1.1) is not only an infinite-dimensional stochastic Hamiltonian system and possesses the symplectic structure, but also a Hamiltonian partial differential equation and possesses stochastic multi-symplectic conservation law. To inherit stochastic symplectic structure, we propose the temporal semi-discretization via the mid-point scheme. Due to the fact that the nonlinear term of (1.1) is not global Lipschitz, it is difficult to analyze the convergence order of the proposed scheme. Motivated by , we use the truncated argument to get the truncated equation whose nonlinear term is global Lipschitz. Then we discretize the truncated equation using mid-point scheme and deduce the order one of this scheme in mean square sense. Furthermore, it yields order one in probability for the proposed scheme approximating the original equation. To preserve the stochastic multi-symplectic conservation law, we obtain a stochastic multi-symplectic scheme via the classical mid-point scheme both in temporal and spatial directions. The evolution laws of discrete charge and energy of the scheme are given, and they satisfy the property of linear growth with respect to time in the sense of expectation, which coincide with the theoretical results.

The rest of the paper is organized as follows. Some properties of the solution, including the evolution law of charge, the uniform boundedness of energy and solution, are given in Section 2. In Section 3, we first show that (1.1) owns the stochastic symplectic structure, then we construct a stochastic symplectic scheme and prove that its temporal order of convergence is one in probability. Section 4 is contributed to proving that (1.1) possesses the stochastic multi-symplectic conservation law and we give a stochastic multi-symplectic scheme. Both charge and energy evolution laws of this scheme are analyzed. In Section 5, we perform numerical experiments to verify our theoretical results.

In the remainder of the article, is a generic constant whose value may vary in different occurrences, denotes the constant depending on some parameters.

## 2. Stochastic NLS equation with quadratic potential

In order to state precisely Eq. (1.1), we consider the probability space endowed with a normal filtration . Let be a sequence of independent real valued Brownian motions on , and be an orthonormal basis of some Hilbert space . We consider the complex-valued Wiener process

 W=∑k∈Nβk(t,ω)ϕek(x),t∈[0,T],x∈D,ω∈Ω,

where the space of Hilbert-Schmidt operators from to another Hilbert space . The corresponding norm is then given by

 ∥ϕ∥2L2(U,H)=tr(ϕ∗ϕ)=∑k∈N∥ϕek∥2H.

In addition, denotes Hilbert space with inner product for . is denoted by , where is Sobolev space consisting of functions such that exist and are square integrable for all is positive integer. Throughout the paper, we assume that for a certain parameter , and Sobolev space will be used .

Now, we recall the mild solution of Eq. (1.1) from .

An -valued -adapted process is called a mild solution of (1.1) if for every holds -a.s.

 u(t)= S(t)u0+iθ∫t0S(t−r)|x|2u(r)dr (2.1) +iλ∫t0S(t−r)|u(r)|2σu(r)dr−i∫t0S(t−r)dW(r),

where denotes the semigroup of solution operator of the deterministic linear differential equation

 idu+Δudt=0  in  D,  u=0  %on  ∂D×(0,T),  u(0)=u0  in  D.

If the noise term is eliminated in (1.1), then it is the deterministic nonlinear Schrödinger equation with the quadratic potential

 idu+(Δu+θ|x|2u+λ|u|2σu)dt=0,

it possesses the charge conservation law

 (2.2) M(u(t))=∫D|u(t,x)|2dx=∫D|u(t0,x)|2dx=M(u0),

and energy conservation law

 (2.3) H(u(t))= 12∫D|∇u(t,x)|2dx−θ2∫D|xu(t,x)|2dx−λ2σ+2∫D|u(t,x)|2σ+2dx = H(u0).

But these conservation laws of Eq. (1.1) are invalid. Now, we state its charge evolution law and energy evolution law, respectively.

###### Lemma 2.1.

Eq. (1.1) has the following global charge evolution law a.s.

 (2.4) M(u(t))=M(u0)−2I∑k∈N∫t0∫Du(¯¯¯¯¯¯¯¯ϕek)dxdβk(s)+t∥ϕ∥2L02,t≥0.

Moreover, for any there exists a constant

 (2.5) Esupt∈[0,T]Mm(u(t))≤ CmEMm(u0).
###### Proof.

The proof is based on the application of Itô formula to functional Since is Fréchet derivable, the derivatives of along directions and are as follows:

 DM(u)(φ)=2R∫Du¯¯¯¯φdx,D2M(u)(φ,ψ)=2R∫Dφ¯¯¯¯ψdx.

From Itô formula, we have

 M(u)= M(u0)+∫t0DM(u(s))du+12∫t0tr[D2M(u)(−iϕ)(−iϕ)∗]ds = M(u0)+2∫t0R∫Du(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯iΔu+iθ|x|2u+iλ|u|2σu)dxds +2∫t0R∫Du(−i)dxd¯¯¯¯¯¯¯¯¯¯¯¯¯W(s)+∫t0tr(ϕϕ∗)ds = M(u0)−2I∫t0∫Dudxd¯¯¯¯¯¯¯¯¯¯¯¯¯W(s)+t∥ϕ∥2L02 = M(u0)−2I∑k∈N∫t0∫Du¯¯¯¯¯¯¯¯ϕekdxdβk(s)+t∥ϕ∥2L02.

To prove (2.5), we apply Itô formula to

 Mm(u)= Mm(u0)−2mI∑k∈N∫t0Mm−1(u)∫Du¯¯¯¯¯¯¯¯ϕekdxdβk(s)+m∥ϕ∥2L02∫t0Mm−1(u)ds +2m(m−1)∫t0Mm−2(u)R∑k∈N(∫Du¯¯¯¯¯¯¯¯ϕekdx)2ds,t≥0.

The inequality (2.5) follows from Proposition 3.2 of . ∎

Taking expectation in both sides of (2.4), we have that

 (2.6) EM(u)=EM(u0)+t∥ϕ∥2L02,t≥0.

The formula (2.6) indicates that the average charge is linear growth with respect to time Next, we present the energy evolution law of (1.1), which has established in .

###### Lemma 2.2.

Eq. (1.1) has the following global energy evolution law a.s.

 H(u(t))= H(u0)+I∫D∫t0(Δu+θ|x|2u+λ|u|2σu)d¯¯¯¯¯¯¯¯¯¯¯¯¯W(s)dx+12t∑k∈N∥∇ϕek∥2L2 −θ2t∑k∈N∫D|x|2|ϕek|2dx−λ2∑k∈N∫t0∫D|u|2σ|ϕek|2dxds (2.7) −σλ∑k∈N∫t0∫D|u|2σ−2(I(¯¯¯u(ϕek)))2dxds,t≥0.

The average energy is also linear with respect to time at most if we take expectation in both sides of the above formula (2.2). Under the assumptions that and is a bounded interval, we have the following inequalities related to the energy , the proof of this lemma is an immediate consequence of Gagliardo–Nirenberg’ inequality (See e.g. ).

###### Lemma 2.3.

Assume that . For , there exist constants and such that
(i) if , then
(ii) if , then

Using this lemma, we get the uniform boundedness of .

###### Lemma 2.4.

Let , , . There exists a constant such that

 (i)  sup0≤t≤TE(H(u(t)))p≤C, (ii)  Esup0≤t≤T(H(u(t)))p≤C.
###### Proof.

We first consider the case of . If applying the expectation to (2.2), we have

 EH(u(t))≤ EH(u0)+12tE∑k∈N|∇ϕek|2L2+|θ|2tE∑k∈N∫D|x|2|ϕek|2dx ≤

the assertion (i) holds. If

 EH(u(t))≤ EH(u0)+12tE∑k∈N∥∇ϕek∥2L2+|θ|2tE∑k∈N∫D|x|2|ϕek|2dx +|λ|E∑k∈N∫t0∫D(12|u|2σ|ϕek|2+σ|u|2σ−2(I(¯¯¯u(ϕek)))2)dxds ≤

Since Hölder inequality, Young’s inequality and Gagliardo–Nirenberg’ inequality, we have

 ∑k∈N∫D|u|2σ|ϕek|2dx ≤∥u∥2σL2σ+2∥ϕ∥2L2(U,L2σ+2)≤σσ+1∥u∥2σ+2L2σ+2+1σ+1∥ϕ∥2σ+2L2(U,L2σ+2) ≤σ2∥∇u∥2L2+C(σ)∥u∥4+2σ2−σL2+1σ+1∥ϕ∥2σ+2L2(U,L2σ+2) ≤σH(u)+C(θ)∥u∥2L2+C∥u∥4+2σ2−σL2+1σ+1∥ϕ∥2σ+2L12,

the last inequality follows from Lemma 2.3 and is embedded into . Therefore,

 EH(u(t))≤Ct+CE∫t0H(u(s))ds.

Then we use Gronwall’s inequality to get the assertion (i).

If we apply Itô formula to then

 (H(u))p= +12∫t0p(p−1)(H(u))p−2∑k∈N(I∫D(Δu+θ|x|2u+λ|u|2σu)¯¯¯¯¯¯¯¯ϕekdx)2ds +12∫t0p(H(u))p−1(∑k∈N∥∇ϕek∥2L2−θ∑k∈N∫D|x|2|ϕek|2dx (2.8) −λ∑k∈N∫D|u|2σ|ϕek|2dx−2σλ∑k∈N∫D|u|2σ−2(I(¯¯¯u(ϕek)))2dx)ds.

Since the second term on the right-hand side vanishes after taking expectation, there remains to estimate the third term and the last term. For the third term, there exists a constant such that

 (2.9) ∑k∈N(I∫D((Δu+θ|x|2u+λ|u|2σu)¯¯¯¯¯¯¯¯ϕek)dx)2≤C∥ϕ∗(Δu+θ|x|2u+λ|u|2σu)∥2L2.

The operator is bounded from into . Furthermore, is embedded into , is also bounded from into and Young’s inequality, so we obtain

 ∥ϕ∗(Δu+θ|x|2u+λ|u|2σu)∥2L2≤∥ϕ∥2L12(∥∇u∥L2+|θ||x|2∥u∥L2+|λ|∥u∥2σ+1L2σ+2)2 ≤3∥ϕ∥2L12(∥∇u∥2L2+θ2|x|4∥u∥2L2+λ2∥u∥4σ+2L2σ+2) (2.10) ≤3(∥ϕ∥4L12+λ22σ+2∥ϕ∥4σ+4L12+12∥∇u∥4L2+θ4|x|82∥u∥4L2+λ2(2σ+1)(2σ+2)∥u∥4σ+4L2σ+2).

From Gagliardo–Nirenberg’ inequality and

 (2.11) ∥u∥4σ+4L2σ+2≤C∥∇u∥2σL2∥u∥2σ+4L2≤C(σ)(∥∇u∥4L2+∥u∥4σ+82−σL2).

Substituting (2.11) into (2), and combining with (2.5) and Lemma 2.3, we have

 ∑k∈N(I∫D(Δu+θ|x|2u+λ|u|2σu)¯¯¯¯¯¯¯¯ϕekdx)2 ≤C(|θ|,|λ|,σ)(∥ϕ∥4L12+∥ϕ∥4σ+4L12+∥∇u∥4L2+∥u0∥4L2+∥u0∥4σ+82−σL2) (2.12) ≤C(|θ|,|λ|,σ)((H(u))2+∥ϕ∥4L12+∥ϕ∥4σ+4L12+∥u0∥4L2+∥u0∥4σ+82−σL2).

For the last term in (2), due to Young inequality, Lemma 2.3 and is embedded into we get

 ∑k∈N∥∇ϕek∥2L2−θ∑k∈N∫D|x|2|ϕek|2dx −λ∑k∈N∫D|u|2σ|ϕek|2dx−2σλ∑k∈N∫D|u|2σ−2(I(¯¯¯u(ϕek)))2dx ≤∥ϕ∥2L12+C(|θ|)∥ϕ∥2L12+|λ|(2σ+1)∥u∥2σL2σ+2∥ϕ∥2L2σ+2 (2.13) ≤C(|θ|,|λ|,σ)(H(u)+∥ϕ∥2L12+∥ϕ∥2σ+2L12+∥u0∥2L2+∥u0∥2σ+42−σL2).

Because of (2), (2), and Hölder inequality, we deduce

 E(H(u(t)))p≤Ct+CE∫t0(H(u(s)))pds,

we apply Gronwall’s inequality to obtain the estimate (i).

The assertion (ii) uses the similar arguments with (i) and Theorem 3.4 of , so we skip the details here. ∎

In , Theorem 4.6, a uniform boundedness for the energy is used to construct a unique mild solution with continuous -valued paths for stochastic nonlinear Schrödinger equation. We can follow the same strategy to construct the unique global mild solution with continuous-valued paths with under . Moreover, we have the following uniform boundedness of the solution under . Similar to Theorem 2.1 of , the proof of stability of the solution in Sobolev space is directly obtained by analyzing the functional

 f(u)=∥∇4u∥2L2−λR∫D((−Δ)3u)(|u|2¯u)dx.
###### Lemma 2.5.

Let , , and . There exists a constant such that

 Esup0≤t≤T∥u∥2pH4≤C.

## 3. Stochastic symplectic scheme

As we all know, the stochastic Schrödinger equation without quadratic potential is an infinite-dimensional stochastic Hamiltonian system, which characterize the geometric invariants of the phase flow and contribute to constructing the numerical schemes for long time computation. In the following, we show that Eq. (1.1) possesses stochastic symplectic structure.

Denote by and the real and imaginary parts of , respectively. Set with and being real–valued functions. Then Eq. (1.1) is equivalent to

 dtp+Δqdt+θ|x|2qdt+λ|p2+q2|σqdt=dW2, (3.1) dtq−Δpdt−θ|x|2pdt−λ|p2+q2|σpdt=−dW1

with initial datum Set

 H1=−12∫D|∇u|2dx+θ2∫D|xu|2dx+λ2σ+2∫D|u|2σ+2dx.

Then the equation (3) can be rewritten as

 dtp=−δH1δqdt+dW2, (3.2) dtq=δH1δpdt−dW1,

where and denote the variational derivative of with respect to and , respectively. In fact, (3) is an infinite-dimensional stochastic Hamiltonian system. Using the same procedure as , one can derive that (3) possesses the symplectic structure

 (3.3) ¯¯¯ω(t):=∫x1x0dp∧dqdx.
###### Theorem 3.1.

The phase flow of Eq. (1.1) preserves the symplectic structure (3.3).

In order to preserve the stochastic symplectic structure, we consider the mid-point scheme of the temporal discretization for (1.1),

 (3.4) iun+1−unτ+Δun+12+θ|x|2un+12+λ|un+12|2σun+12=△nWτ,  n=0,⋯,N−1,

where is time step size, Denote

 ¯¯¯ωn:=∫x1x0dpn∧dqndx,

we have the following result.

###### Theorem 3.2.

The scheme (3.4) possesses the discrete symplectic structure, i.e.

 ¯¯¯ωn+1=¯¯¯ωn.
###### Proof.

For convenience, denote then (3.4) can be rewritten as

 pn+1=pn−Δqn+12τ−θ|x|2qn+12τ−λ∂Φ∂q(pn+12,qn+12)τ+△nW2, qn+1=qn+Δpn+12τ+θ|x|2pn+12τ+λ∂Φ∂p(pn+12,qn+12)τ−△nW1.

Differentiating the above equation on the phase space, we obtain

 dpn+1=dpn−Δdqn+12τ−θ|x|2dqn+12τ−λ∂2Φ∂q∂pdpn+12τ−λ∂2Φ∂q2dqn+12τ, dqn+1=dqn+Δdpn+12τ+θ|x|2dpn+12τ+λ∂2Φ∂p∂qdqn+12τ+λ∂2Φ∂p2dpn+12τ.

Then

 dpn+1∧dqn+1 =dpn∧dqn+(−Δdqn+12τ−θ|x|2dqn+12τ−λ∂2Φ∂q∂pdpn+12τ−λ∂2Φ∂q2dqn+12τ)∧dqn +dpn∧(Δdpn+12τ+θ|x|2dpn+12τ+λ∂2Φ∂p∂qdqn+12τ+λ∂2Φ∂p2dpn+12τ) +(−Δdqn+12τ−θ|x|2dqn+12τ−λ∂2Φ∂q∂pdpn+12τ−λ∂2Φ∂q2dqn+12τ) ∧(Δdpn+12τ+θ|x|2dpn+12τ+λ∂2Φ∂p∂qdqn+12τ+λ∂2Φ∂p2dpn+12τ).

Substituting

 dpn=dpn+12+Δdqn+12τ+θ|x|2dqn+12τ+λ∂2Φ∂q∂pdpn+12τ+λ∂2Φ∂q2dqn+12τ, dqn=dqn+12−Δdpn+12τ−θ|x|2dpn+12τ−λ∂2Φ∂p∂qdqn+12τ−λ∂2Φ∂p2dpn+12τ

into the above equality, we have

 ∫x1x0dpn+1∧dqn+1dx=∫x1x0dpn∧dqndx,  a.s.

In fact, the existence of solution of the scheme (3.4) follows from a standard Galerkin method and Brouwer theorem (see ).

###### Proposition 3.1.

Let and be -measurable with values in , then for sufficiently small , there exists an -valued -adapted solution of (3.4) .

###### Proof.

Fix a family of deterministic functions in , we also fix the existence of solution of

 (3.5) i˜un+1−˜unτ+Δ˜un+12+θ|x|2˜un+12+λ|˜un+12|2σ˜un+12=ξn

follows from a standard Galerkin method and Brouwer theorem, see Lemma 3.1 in  and Lemma 3.1 in . Assuming that is a solution of (3.5), multiplying (3.5) by the complex conjugate of