Mass- and energy-conserved numerical schemes for nonlinear Schrödinger equations
Xiaobing Feng Hailiang Liu Shu Ma























shapes.geometric




thmTheorem























schemeScheme
remarkRemark




 §.§ 
[1]#1 



  §.§.§ 
[1]#1 


C[1]>
#1

m#1<




graphics/	
 



Mass- and energy-conserved numerical schemes for nonlinear Schrödinger equations
    
Xiaobing FengDepartment of Mathematics, The University of
Tennessee, Knoxville, TN 37996, U.S.A.  (xfeng@math.utk.edu)



Hailiang LiuDepartment of Mathematics, Iowa State University, Ames, IA 50011, U.S.A.
(hliu@iastate.edu).
Shu MaDepartment of Applied Mathematics, Northwestern Polytechnical University, 
Xian, Shaanxi, 710065, P. R. China.
(mashu@mail.nwpu.edu.cn). 

    December 30, 2023
========================================================================================================================================================================================================================================================================================================================================================================================


In this paper, we propose a family of time-stepping schemes for approximating 
general nonlinear Schrödinger equations. 
The proposed schemes all satisfy both mass and energy conservation (in a modified form for the latter).
Truncation and dispersion error analyses are provided for four proposed schemes. 
Efficient fixed-point iterative solvers are also constructed to solve 
the resulting nonlinear discrete problems. As a byproduct, an efficient one-step implementation 
of the BDF schemes is obtained as well. 

Extensive numerical experiments are presented to demonstrate 
the convergence and the capability of capturing the blow-up time of the proposed schemes.

	

Nonlinear Schrödinger equations, mass conservation and energy conservation,
BDF schemes, finite element methods, finite time blow-ups.

	

65M06, 65M12


myheadings
plain
XIAOBING FENG, HAILIANG LIU and SHU MA
ENERGY-CONSERVED NUMERICAL METHODS FOR NLS EQUATIONS 




§ INTRODUCTION.

In this paper we consider the following nonlinear 
Schrödinger (NLS) equation:   


i u_t    = -Δu +λf(|u|^2) u            _T: = ×(0,T),




u (0)    = u _0           , 

where ⊂𝐑^d  (d=1,2,3) is a bounded domain, λ=±1, T>0 and 
i=√(-1) stands for the imaginary unit.  
u=u(x,t): _T → is  a complex-valued function. 
f: 𝐑_+→𝐑_+ is a given real-valued function,
which could be different in different applications, e.g., see
<cit.> and the references therein. 
The best known f is f(s)=s, which leads to the well-known nonlinear 
Schrödinger equation with cubic nonlinearity. To close the system, we 
also need to specify a boundary condition for u. In this paper we consider both  
homogeneous Dirichlet and periodic boundary condition (see Sections <ref> and <ref>), 
although most of our derivations and proofs are independent of the boundary condition. 
 
The above Schrödinger equation  
describes many physical phenomena in optics, mechanics, and plasma physics.
Mathematically, the NLS equation is a prototypical dispersive wave equation, 
its solutions exhibit some intriguing properties such as energy conservation, 
soliton wave, and possible blow-ups <cit.>. 
In particular, the equation preserves both the mass and the Hamiltonian energy, 
that is, the following quantities are constants in time:

    (u)(t)    := u(t)_L^2^2=∫_ |u(t)|^2  dx, 
    (u)(t)    :=∫_(|∇ u(t)|^2+λ F(|u(t)|^2))  dx,
       F(s)=∫_0^s f(μ)  dμ.

Here the dependence of u on x variable is suppressed for notational brevity. 
The case with positive λ is called defocusing and with negative λ is called focusing 
which allows for bright soliton solutions as well as breather solutions. 

Dispersion and nonlinearity can interact to produce permanent and localized wave forms 
in nonlinear dispersive wave equations such as the Korteweg-de Vries (KdV) equation
<cit.> 
and the cubic Schrödinger equation <cit.>. 
A distinct feature of these equations is the infinite many conservation laws 
(conserved integrals as invariants), allowing for soliton solutions which emerge from 
collision unchanged over time. The quality of the numerical approximation hence hinges 
on how well the conserved integrals can be preserved at the discrete level. 
Numerical methods without this property may result in substantial phase and shape 
errors after long time integration.
Indeed for some wave equations the invariant preserving high order numerical methods 
have been shown more accurate than non-conservative methods after long-time numerical 
integration (see  <cit.>). 

For the nonlinear Schrödinger equation considered in this paper,  
a natural question is whether it is possible to design numerical schemes which conserve the mass  and 
energy  simultaneously.  A lot of effort has been made to preserve the mass by high order 
spatial discretization such as spectral methods <cit.>,   
and discontinuous Galerkin methods <cit.>.  
A modified numerical energy may be preserved approximately by the corresponding spatial discretization 
(see <cit.>). However, since those methods are based on a time-splitting technique <cit.>, 
except the Crank-Nickson scheme, they are only mass-conserved. The objective of this work is to provide an attempt to address the above  question. Specifically, in this paper we develop and analyze a family of  
mass- and energy-conserved time-stepping schemes for approximating the cubic and general 
nonlinear Schrödinger equations. It should be noted that the energy conservation is also achieved for a 
modified energy, instead of the original energy. 

The rest of this paper is organized as follows. In Section <ref> we present a general framework 
involving two sequences of time-stepping schemes, which is shown to preserve both mass and 
energy for arbitrary time-step sizes, for the cubic nonlinear Schrödinger equation. 
In Section <ref> we present many specific examples of mass- and (modified) energy-conserved time 
stepping schemes which fit into the general framework, and derive the truncation errors for four of these schemes.
In Section <ref> we present an efficient iterative algorithm to solve the resulting 
nonlinear equations. 
In Section <ref> we extend the framework and examples to the general Schrödinger equations 
with arbitrary nonlinearity. 
In Section <ref> we present a dispersion error analysis and derive the convergence rates 
for the dispersion errors.
In Section <ref> we present numerical experiments to validate the theoretical results and 
to gauge the performance of the proposed schemes, especially the sharpness of the convergence rates. 
In Section 8 we present additional numerical experiments to demonstrate the capability of the proposed 
numerical schemes for resolving the blow-up
phenomenon. The paper is completed with some concluding remarks and comments 
given in Section <ref>.
 



§ SEMI-DISCRETIZATION IN TIME: A GENERAL FRAMEWORK.

In this section we propose a family of energy-conserved time-stepping schemes 
for approximating the cubic nonlinear Schrödinger equation.

Let τ>0 and t=t_n=nτ for n=0,1,2,⋯, N be a uniform mesh for [0,T]. 
Let k be a positive integer. We propose the following general k-step time-stepping scheme 
for problem (<ref>)–(<ref>): 
Seeking {R^n,u^n} for n=k,k+1,⋯, N such that 

    id_t R^n+1 =   -Δ R^n+1/2 + λ/2(|R^n|^2+|R^n+1|^2)R^n+1/2,
    
    u^n+1=   β_0^-1 (R^n+1-∑_j=1^k-1β_j u^n+1-j)         k>1,  
    β_0^-1 R^n+1        k=1,

where we use notation 

    d_t R^n+1=R^n+1-R^n/τ,      R^n+1/2=R^n+1+R^n/2.


Note that from (<ref>) we see that 
R^n is a linear combination of u^n,u^n-1,⋯,u^n-k+1 given as follows:

    R^n =∑_j= 0^ k-1β_j u^n-j.


As expected, choosing the parameters {β_j} is a delicate issue. The guideline we use (see the details in
the next section) is to choose them such that d_t R^n+1 has a specific order of accuracy for approximating  u_t(t_n+1).  In Table <ref> we list several sets of parameters {β_j} to be used in  
scheme (<ref>)–(<ref>), which results in various specific schemes for the nonlinear 
Schrödinger equation (again, see the details in the subsequent section). 



We also remark that scheme (<ref>)–(<ref>) produces two sequences, 
namely {R^n} and {u^n}.
The first sequence can be regarded as the auxiliary quantities 
which are generated by solving nonlinear equation (<ref>), 
while the second sequence are obtained 
as linear combinations of the first one. 

To prove a key mass- and energy-conservation property of scheme (<ref>)–(<ref>),
we define the following discrete mass and energy

    ^n:=R^n_L^2^2,    ^n:= 1/2∇ R^n_L^2^2 + λ/4R^n_L^4^4.


We start with establishing the following conservation results for the solution of scheme (<ref>)–(<ref>).

  
The solution to scheme (<ref>) and (<ref>) satisfies
^n = ^0 and ^n = ^0 for all n ≥ 1.



We multiply equation (<ref>) by R̅^n+1/2,
integrate it over  and use the integration by parts to get, 

    i/2τ∫_(R^n+1-R^n)
    (R̅^n+1+R̅^n)  dx 
            =∇ R^n+1/2_L^2^2 
    + λ/2∫_(|R^n|^2+|R^n+1|^2 )|R^n+1/2|^2   dx.

Taking the imaginary part of the resulting equation to get

    1/2τ[∫_(R^n+1-R^n)
    (R̅^n+1+R̅^n)  dx]=0.

It follows from the identity [(a-b)(a̅+b̅)] =|a|^2-|b|^2 
for a,b ∈ that

    ^n+1 = ^n,   ∀  n≥ 0.


To show the second conservation property, we multiply (<ref>) by d_t R̅^n+1,
integrate the equation over  and use the integration by parts to get

    id_tR^n+1_L^2^2
    =   1/2τ∫_(∇ R^n+1-∇ R^n)
    (∇R̅^n+1+∇R̅^n)  dx
       + λ/4τ∫_(|R^n|^2+|R^n+1|^2)(R^n+1-R^n)
    (R̅^n+1+R̅^n)  dx.

Taking the real part on the equation and applying identity 
[(|a|^2+|b|^2)(a-b)(a̅+b̅)] =|a|^4-|b|^4 for a,b ∈ yield 

    ^n+1 = ^n,   ∀  n≥ 0.
 
The proof is completed.





§ SPECIFIC SCHEMES AND THEIR TRUNCATION ERROR ANALYSIS.

In this section we propose a number of specific schemes by defining R^n in terms of 
u^n, u^n-1,⋯, u^n-k+1. In other words, we shall specify the choice of 
parameters {β_j}_j=0^k-1 for each scheme. 



 §.§ A modified Crank-Nicolson scheme (k=1).
 
By setting R^n=u^n, the modified Crank-Nicolson scheme is defined as

    i/τ(u^n+1-u^n) = -1/2Δ(u^n+1+u^n) 
    +λ/4(|u^n|^2 + |u^n+1|^2)
    (u^n+1 + u^n).


The local truncation error (LTE) of the modified  Crank-Nicolson scheme 
(<ref>) is defined by

    TE^n =    i1/τ(u(t_n+1)-u(t_n))
    + 1/2Δ(u(t_n+1) + u(t_n))
       -λ/4(|u(t_n)|^2 + |u(t_n+1)|^2)
    (u(t_n+1) + u(t_n))
    
    
    =   :_1+_2+_3,

where u(t) is the true solution of the nonlinear Schrödinger equation as follows:

    i u_t(t) +Δ u(t) -λ |u(t)|^2 u(t) = 0.

Here and below the dependence of u on x variable is suppressed for notational brevity.

The first lemma establishes the local truncation error for the modified Crank-Nicolson scheme.


The local truncation error of the modified Crank-Nicolson scheme (<ref>) (k=1) is O(τ^2).



By using a Taylor series expansion of u(t_n+1) and u(t_n) about u(t_n+1/2),
the accurate order of derivative term _1 and Laplace term _2 are O(τ^2).

    _1
    =   i/τu(t_n+1/2) + i/2u_t(t_n+1/2)
    +iτ/8u_tt(t_n+1/2) + iτ^2/2^33!u^(3)(t_n+1/2)+ O(τ^3)
    
       -i/τu(t_n+1/2) + i/2u_t(t_n+1/2)
    -iτ/8u_tt(t_n+1/2) + iτ^2/2^33!u^(3)(t_n+1/2)+ O(τ^3)
    
    
    =     iu_t(t_n+1/2)+O(τ^2).
    _2
    =   1/2Δ u(t_n+1/2) + τ/4(Δ u)_t (t_n+1/2)
    +τ^2/16(Δ u)_tt (t_n+1/2) + O(τ^3)
    
       +1/2Δ u(t_n+1/2) - τ/4(Δ u)_t (t_n+1/2)
    +τ^2/16(Δ u)_tt(t_n+1/2)+ O(τ^3) 
    
    =   Δ u(t_n+1/2)+O(τ^2).

Note that

    |u(t_n)|^2
    =   (u(t_n+1/2)-τ/2u_t(t_n+1/2)
    + O(τ^2))
    (u̅(t_n+1/2)-τ/2u̅_t(t_n+1/2)
    + O(τ^2))
    
    
    =   |u(t_n+1/2)|^2
    -τ/2u̅(t_n+1/2)u_t(t_n+1/2)
    -τ/2u̅_t(t_n+1/2)u(t_n+1/2)+O(τ^2),
    
    |u(t_n+1)|^2
    =   (u(t_n+1/2)+τ/2u_t(t_n+1/2)
    + O(τ^2))
    (u̅(t_n+1/2)+τ/2u̅_t(t_n+1/2)
    + O(τ^2)), 
    
    
    =   |u(t_n+1/2)|^2
    +τ/2u̅(t_n+1/2)u_t(t_n+1/2)
    +τ/2u̅_t(t_n+1/2)u(t_n+1/2)+O(τ^2).

which leads to 

    _3=   -λ/4(2|u(t_n+1/2)|^2
    +O(τ^2))
    (2u(t_n+1/2)+O(τ^2))
    
    =   -λ |u(t_n+1/2)|^2 u(t_n+1/2) +O(τ^2).


Recall the definition of the true solution u(t) 
and using (<ref>)-(<ref>) in (<ref>) gives

    TE^n=   [iu_t(t_n+1/2)
    +Δ u(t_n+1/2)
    -λ |u(t_n+1/2)|^2 u(t_n+1/2) ]
    +O(τ^2)
    =O(τ^2).

Thus the local truncation error of above Crank-Nicolson scheme is O(τ^2). The proof is completed. 




 §.§ A Leapfrog scheme (k=2).
 
For k=2 we set β_0=β_1=1/2 in (<ref>) so that 
R^n=1/2(u^n+u^n-1),
which inserted into (<ref>)) leads to the following  Leapfrog scheme:

    i/2τ(u^n+1-u^n-1)
    =-1/4(Δ u^n+1
    +2Δ u^n+Δ u^n-1) 
           + λ/8(|u^n+u^n-1/2|^2
    +|u^n+1+u^n/2|^2)
    (u^n+1 + 2u^n+ u^n-1).


The local truncation error of the leapfrog scheme  (<ref>) is defined by 

    TE^n =    i/2τ(u(t_n+1) - u(t_n-1))
    + 1/4Δ( u(t_n+1)
    + 2 u(t_n) + u(t_n-1))
       -λ/8[|u(t_n) + u(t_n-1)/2|^2
    +|u(t_n+1) + u(t_n)/2|^2]
    (u(t_n+1) + 2u(t_n)+u(t_n-1))
    
    =:   _1+_2+_3.


The next lemma establishes the local truncation error for the above Leapfrog scheme.

The local truncation error of the Leapfrog scheme (k=2) (<ref>) is O(τ^2).



We formally apply the Taylor series expansions of u(t_n+1) and u(t_n-1)
about u(t_n) for _1 to get

    _1
    =   i/2τu(t_n) + i/2u_t(t_n)
    +iτ/4u_tt(t_n) + iτ^2/12u^(3)(t_n) + O(τ^3)
    
       -i/2τu(t_n) + i/2u_t(t_n)
    -iτ/4u_tt(t_n) + iτ^2/12u^(3)(t_n)+ O(τ^3) 
    
    =    iu_t(t_n) + O(τ^2).


Applying the Taylor series expansions of u(t_n+1) and u(t_n-1)
about u(t_n) again for _2, we have

    _2
    =   1/4Δ u(t_n) + τ/4(Δ u)^'(t_n)
    +τ^2/8(Δ u)^”(t_n)
    +τ^3/24(Δ u)^(3)(t_n) + O(τ^4)
       + 1/2Δ u(t_n)
    +1/4Δ u(t_n)
    -τ/4(Δ u)_t(t_n)
    +τ^2/8(Δ u)_tt(t_n)
       -τ^3/24(Δ u)^(3)(t_n) + O(τ^4)
    
    
    =   Δ u(t_n) + O(τ^2).


By using the following facts for _3

    |u(t_n)+u(t_n-1)/2|^2
    =   1/4(u(t_n)+u(t_n-1))(u̅(t_n)+u̅(t_n-1)) 
    
    =   1/4(2u(t_n)-u_t(t_n)τ
    + O(τ^2))
    (2u̅(t_n)-u̅_t(t_n)τ
    + O(τ^2)) 
    
    
    =   |u(t_n)|^2
    -1/2τu̅(t_n)u_t(t_n)
    -1/2τu̅_t(t_n)u(t_n) + O(τ^2),
    
    |u(t_n)+u(t_n+1)/2|^2
    =   1/4(u(t_n)+u(t_n+1))(u̅(t_n)+u̅(t_n+1)) 
    
    =   |u(t_n)|^2
    +1/2τu̅(t_n)u_t(t_n)
    +1/2τu̅_t(t_n)u(t_n) + O(τ^2),

the order of nonlinear term can be estimated as

    _3
    =   -λ/8(2|u(t_n)|^2
    +O(τ^2))
    (4u(t_n) + O(τ^2))
    ≈    -λ |u(t_n)|^2u(t_n) + O(τ^2).


Combing (<ref>)-(<ref>) together in (<ref>) leads to

    TE^n = [iu_t(t_n) +Δ u(t_n) -λ |u(t_n)|^2 u(t_n)] + O(τ^2)
    = O(τ^2).

Thus the local truncation error of above Leapfrog scheme is O(τ^2). The proof is completed.





 §.§ Modified BDF schemes (k=s).
 
Let s>1. We recall that the s-step BDF scheme approximates the time derivative 
u_t (t_n+1) as follows: 

    u_t(t_n+1) ≈   1/τ∑_j = 0^s  u^n+1-j,

where  are given in Table <ref>.   
Our idea is to rewrite the above BDF expression as a first order backward difference, that is,

    1/τ∑_j = 0^s  u^n+1-j = 1/τ(R^n+1-R^n) 
          
    R^n=∑_j = 0^s-1 u^n-j,

which turns out is possible. Somehow this simple reformulation has not been seen 
in the literature before. 
  
In order to determine the coefficients , we solve them using the following identity:

    β_s,0u^n+1 + ∑_j = 1^s-1 (-β_s,j-1)u^n+1-j
    -β_s,s-1u^n+1-s
    = ∑_j = 0^s  u^n+1-j.

This holds true for all u^j as long as the following matrix equation is satisfied:

    [ 
    [  1; -1  1;     ⋱  ⋱;    -1  1;     ⋱  ⋱;    -1  1;    -1 ]]_(s+1)× s[ [   β_s,0;   β_s,1;       ⋮;        ;       ⋮; β_s,s-1 ]] 
    = [ [   α_s,0;   α_s,1;       ⋮;        ;       ⋮; α_s,s-1;   α_s,s ]].


It is easy to check that 

    = ∑_ℓ = 0^j α_s,ℓ,      j=0,1,2,⋯, s-1.

Thus our modified BDF schemes are defined as

    i/τ(∑_j = 0^s-1    u^n+1-j
    -∑_j = 0^s-1 u^n-j)
    
    =   -Δ R^n+1/2 
    + λ/2(|∑_j = 0^s-1 u^n+1-j|^2
    +|∑_j = 0^s-1 u^n-j|^2 )R^n+1/2,

where 

    R^n+1/2=
    1/2β_s,0  u^n+1
    +1/2(∑_j = 1^s-1
    +∑_j = 1^sβ_s,j-1) u^n+1-j,

and  are given in the following table (note that since BDF methods with s>6 
are not zero-stable, so we only present s-step BDF with s ≤ 6 here).  


Set :=∑_j = 0^s-1 u(t_n-j) as a linear combination of exact solution values, 
then the local truncation error of the modified BDF schemes (<ref>) is defined by

    TE^n=   i/τ(-))
    +1/2(Δ+ Δ)
       -λ/4(||^2
    +||^2)
    (+ )
    
    
    =   _1+_2+_3.



The local truncation error of the modified BDF schemes (<ref>) (k=s) is O(τ^2).



Recall the approximation of the derivative yields _1=iu_t(t_n+1)+O(τ^s).
Using the Taylor series expansion of u(t_n+1-j) about u(t_n) we have  

    u(t_n+1-j)
    =   u(t_n+1) -jτ u_t(t_n+1) + (jτ)^2/2 u_tt(t_n+1) + O(τ^3),

which implies 

    _2
    =   1/2β_s,0Δ u(t_n+1)
    +1/2∑_j = 1^s-1 (+β_s,j-1)
    Δ u(t_n+1-j)
    +1/2β_s,s-1Δ u(t_n+1-s)
    
    
    =   (∑_j = 0^s-1)
    Δ u(t_n+1)
    -
    τ/2(∑_j = 1^s-1j +∑_j = 1^s jβ_s,j-1)
    (Δ u)_t(t_n+1)
    
       +τ^2/4(∑_j = 1^s-1j^2 
    + ∑_j = 1^s j^2β_s,j-1)
    (Δ u)_tt(t_n+1)
    + O(τ^3)
    
    
    =   Δ u(t_n+1) + O(τ^2),

where we have used the facts that

    ∑_j = 0^s-1 = 1
        ∑_j = 1^s-1j+∑_j = 1^sjβ_s,j-1
    =∑_j = 0^s-1(2j+1) = 0.


For the nonlinear term, we have the following estimates:

    ||^2
    =   | ∑_j = 1^s β_s,j-1 u(t_n+1-j)|^2
    
    =   |∑_j = 1^s β_s,j-1 
    u(t_n+1)
    -τ∑_j = 1^s jβ_s,j-1 
    u_t(t_n+1) 
    +O(τ^2)|^2
    
    
    =   |u(t_n+1)|^2
    -τ∑_j = 1^s jβ_s,j-1 
    u_t(t_n+1) u̅(t_n+1) 
       -τ∑_j = 1^s jβ_s,j-1 
    u(t_n+1) u̅_t(t_n+1)
    +O(τ^2),


    ||^2
    =   |∑_j = 0^s-1 
    u(t_n+1)
    -τ∑_j = 1^s-1 j 
    u_t(t_n+1) 
    +O(τ^2)|^2
    
    
    =   |u(t_n+1)|^2
    - τ∑_j = 1^s-1 j 
    u_t(t_n+1) u̅(t_n+1)
       -τ∑_j = 1^s-1 j 
    u(t_n+1) u̅_t(t_n+1)
    +O(τ^2),

where we use ∑_j = 0^s-1=∑_j = 1^sβ_s,j-1=1, 
∑_j = 1^s-1 j≠ 0  
and ∑_j = 1^s jβ_s,j-1≠ 0.
Since ∑_j = 1^s-1 j + ∑_j = 1^s jβ_s,j-1 = 0, hence,

    ||^2
    +||^2
    =2|u(t_n+1)|^2+O(τ^2).


Similar to the estimates in (<ref>), we have
+  = 2u(t_n+1) + O(τ^2) and

    _3
    =   -λ/4(||^2
    +||^2)
    (+ )
    
    =   -λ/4(2|u(t_n+1)|^2
    +O(τ^2) ) (2u(t_n+1) +O(τ^2)) 
    
    =   -λ |u(t_n+1)|^2 u(t_n+1) +O(τ^2).


Combining (<ref>) and (<ref>) in (<ref>), we obtain

    TE^n=   [iu_t(t_n+1) + Δ u(t_n+1) - λ |u(t_n+1)|^2 u(t_n+1) ] 
    +O(τ^2)
    =O(τ^2).

Thus the local truncation error is O(τ^2).



As a by-product, the above construction also gives a (one-step) backward Euler reformulation 
for BDF schemes. Recall that the s-stage BDF scheme for  u'(t)=f(t, u(t)) is defined as

    ∑_j=0^s α_s,j u^n+1-j = f( t_n+1, u^n+1).

Since

    ∑_j=0^s α_s,j u^n+1-j = 1/τ( R^n+1- R^n )       
    R^n+1 = ∑_j=0^s-1β_j u^n+1-j,
 
then we can rewrite the BDF scheme as

    1/τ( R^n+1- R^n ) = f( t_n+1, R^n+1),

where

    f( t_n+1, R^n+1)
       = f( t_n+1,   β_0^-1( R^n+1- ∑_j=1^s-1β_s,j u^n+1-j)  ), 
    
    u^n+1   = β_0^-1( R^n+1- ∑_j=1^s-1β_s,j u^n+1-j).
 
Hence, each BDF scheme can be implemented as a one-step backward Euler scheme as (<ref>) shows. 





 §.§ A four-step symmetric scheme (k=4).
 
To define this scheme, we set 

    R^n+1=1/12(-u^n+1+7u^n+7u^n-1-u^n-2),

which fits (<ref>) with β=1/12(-1, 7, 7, -1)^⊤.   




The local truncation error of the four-step symmetric scheme (k=4) is O(τ^2).



The derivative term is O(τ^4) as follows

    1/τ(R^n+1-R^n)
    =   1/12τ(-u^n+1+8u^n-8u^n-2+u^n-3)
    ≈    u_t(t_n-1)+O(τ^4).


For the Laplace term, we obtain

    1/2Δ(R^n+1+R^n)
    =   1/24Δ(- u^n+1+6u^n+14u^n-1+6u^n-2-u^n-3)
    ≈   Δ u(t_n-1) + O(τ^2).


For the nonlinear term, we have the following estimates:

    -λ/2(|R^n|^2   +|R^n+1|^2 )R^n+1/2
    
    =   -λ/4(2|u(t_n-1)|^2
    +O(τ^2))
    (2u(t_n-1)+O(τ^2))
    ≈   -λ |u(t_n-1)|^2 u(t_n-1) + O(τ^2).

The proof is complete.




§ AN EFFICIENT FIXED-POINT NONLINEAR SOLVER.
 
To solve the nonlinear equation (<ref>), we adapt the fixed-point iterative algorithm of
<cit.> to the job. The proposed algorithm is defined below.

Algorithm 1


Step 1: Given u^l for l=0,1,...,n-k+1, set

    R^n=∑_j= 0^k-1β_j u^n-j.


Step 2:
Update R^n+1 as follows:  define { w^l}_l=0^L iteratively by solving 

    (iI+ τ/2Δ) w^l+1 
    -λτ/4(|R^n|^2+|2w^l-R^n|^2 )w^l+1
    =i R^n,    l=0, 1, ⋯, L,

such that w^L-w^L-1≤δ for some prescribed tolerance level δ,  then set

    R^n+1=2w-R^n.


Step 3:
Update u^n+1 from R^n+1 by (<ref>), that is, 

    u^n+1 = β_0^-1(R^n+1-∑_j=1^k-1β_j u^n+1-j)         k>1,  
    β_0^-1 R^n+1        k=1.


We note that any spatial discretization method, such as finite element, spectral and discontinuous Galerkin methods,
can be employed in combination with the above algorithm to solve the nonlinear Schrödinegr equation (cf.  <cit.>). 




§ EXTENSIONS TO SCHRÖDINGER EQUATIONS WITH ARBITRARY NONLINEARITY.
 
We consider the following initial-boundary value problem for the general nonlinear 
Schrödinger equation :   


i u_t    = -Δu +λf(|u|^2) u            _T: = ×(0,T),




u (0)    = u _0           , 


We extend scheme (<ref>)–(<ref>) as follows for problem 
(<ref>)–(<ref>):
Seeking {R^n,u^n} for n=k,k+1,⋯, N such that 

    id_t R^n+1 =   -Δ R^n+1/2 
    + λ G(|R^n+1|^2, |R^n|^2 ) R^n+1/2,
    
    u^n+1=   β_0^-1 (R^n+1-∑_j=1^k-1β_j u^n+1-j)         k>1,  
    β_0^-1 R^n+1        k=1,

where G(a,b) is the following two variable function: 

    G(a,b)=F(a)-F(b)/a-b,     F(s):=∫_0^s f(ν)  dν.


For example, G(a,b)= ∑_j=0^p a^p-j b^j if f(s)=s^p+1.  
Again, from (<ref>) we have  

    R^n=∑_j= 0^ k-1β_j u^n-j.

Define 

    ^n_g:=R^n_L^2^2,    ^n_g:= ∇ R^n_L^2^2 + λ F(|R^n|^2)_L^1.

we have the following mass- and energy-conservation property of 
scheme (<ref>)–(<ref>).
 
The solution to scheme (<ref>) and (<ref>) satisfies
^n_g = ^0_g and ^n_g = ^0_g for all n ≥ 1.


The nonlinear solver, Algorithm 1, now is replaced by the following modified algorithm.

Algorithm 2

Step 1: Given u^l for l=0,1,...,n-k+1, set

    R^n=∑_j= 0^k-1β_j u^n-j.


Step 2:
Update R^n+1 as follows:  
define { w^l}_l=0^L iteratively by solving 

    (i I+ τ/2Δ) w^l+1 
    -λτ/2 G(|2w^l-R^n|^2, |R^n|^2)w^l+1
    =i R^n,    l=0, 1, ⋯, L,

such that w^L-w^L-1≤δ for some prescribed tolerance level δ, then set

    R^n+1=2w-R^n.


Step 3:
Update u^n+1 from R^n+1 by (<ref>), that is, 

    u^n+1
    =β_0^-1(R^n+1-∑_j=1^k-1β_j u^n+1-j)         k>1,  
    β_0^-1 R^n+1        k=1.


Again, we remark that any spatial discretization method, such as finite element, spectral and discontinuous Galerkin methods,
can be employed in combination with the above algorithm to solve the nonlinear Schrödinegr equation. 





§ DISPERSION ERROR ANALYSIS.
 
In this section, we analyze the difference between the exact and
numerical dispersion relations for the nonlinear Schrödinger equation 
and investigate ways to reduce the dispersive error generated by 
our mass- and energy-conserved time-stepping scheme (<ref>)–(<ref>). 
To minimize the numerical phase error while accurately solving the Schrödinger equation, 
the idea of preserving dispersion relation equation (DRE), which was proposed earlier   
in the area of computational aeroacoustics by Tam and Webb <cit.>, 
is adopted. We refer the reader to <cit.> for a discussion of other 
structure-preserving algorithms for solving ordinary differential equations. 

Consider the cubic nonlinear Schrödinger equation,

    i u_t +Δ u = λ|u|^2 u.
 
Substituting the plane wave solution u = exp(i(kx-ω t)) into  
equation (<ref>), the relation between the angular frequency
ω and the wave number k is given by <cit.>

    ω = k^2 + λ.
 

To derive the numerical dispersion relation equation for scheme 
(<ref>)–(<ref>), the discrete plane wave solution

of the form u^n = e^i(kx- n ) is utilized,
where  is the numerical angular frequency.


The numerical dispersion relation of the Crank-Nickson scheme 
for the cubic nonlinear Schrödinger equation (<ref>) is given by 

    = 2/arctan((k^2 + λ)/2).


 

Substituting u^n = e^i(kx- n ) 
into the discrete Crank-Nicolson scheme (<ref>) we get

    i/(e^-i- 1)
    = 1/2(k^2 + λ)(e^-i+ 1).


Multiplying the both sides of (<ref>) by e^i/2 to obtain

    i/(e^-i/2 - e^i/2 )
    = 1/2 (k^2 + λ)(e^-i/2
    + e^i/2 ).


By using the identities e^ix-e^-ix = 2isin x and e^ix+e^-ix = 2cos x,  
it follows from (<ref>) that 

    2/sin(/2) = (k^2 + λ)cos(/2).

Hence, (<ref>) holds. The proof is complete.


To analyze the difference between the exact and
numerical dispersions for nonlinear Schrödinger equations,
we define the following dispersion error 

    ω_error := |ω-|/ω.


Table <ref> shows the computed dispersion errors and the convergence order 
for the modified Crank-Nicolson scheme. The numerical results indicate that 
this scheme has a second order dispersion error.




The numerical dispersion relation of the Leapfrog scheme 
for the cubic nonlinear Schrödinger equation (<ref>) is given by 

    2/sin(/2)
    = (k^2 + λcos^2(/2))cos(/2).


 


Setting u^n = e^i(kx -  n ) and using the identities e^ix+e^-ix = 2cos x and  
cos 2θ = 2cos^2θ -1 in |1/2(u^n+u^n-1)|^2 
and |1/2(u^n+u^n+1)|^2 we get

    |1/2(u^n   +u^n-1)|^2 + |1/2(u^n+u^n+1)|^2
    
    =   |1/2 e^ikx-i(n-1/2 ) (e^-i/2 
    + e^i/2)|^2 
       +|1/2 e^ikx-i(n+1/2) (e^i/2+ e^-i/2)|^2 
    
    =   |e^ikx-i(n-1/2) cos(/2)|^2
    +|e^ikx-i(n-1/2) cos(/2)|^2 
    
    =   2cos^2(/2).


Substituting u^n = e^i(kx -  n ) into the Leapfrog scheme (<ref>) 
and using (<ref>), we get the following equality:

    i/(e^-i- e^i)
    =(k^2 + λcos^2(/2))
    (e^-i+1/2 + e^i-1/2).


Using the identities e^ix-e^-ix = 2isin x, e^ix+e^-ix = 2cos x, 
sin 2θ = 2sinθcosθ, and cos 2θ = 2cos^2θ -1 in (<ref>), 
we then obtain the following numerical dispersion relation equation of the Leapfrog scheme
for the cubic nonlinear Schrödinger equation:

    2/sin(/2)
    = (k^2 + λcos^2(/2)) cos(/2).

Hence, (<ref>) holds. The proof is complete.


Table <ref> shows the computed dispersion errors and convergence rates 
for the Leapfrog scheme. The numerical results indicate that the Leapfrog scheme also has a second 
order dispersion error.
   
  



The numerical dispersion relations of the modified BDF schemes  
for the cubic nonlinear Schrödinger equation (<ref>) are given by 

    2/sin(/2)
    = (k^2 + λ H(,))cos(/2).

where H(,) is given in Table <ref> below.

 


We only present the proofs of the numerical dispersion relation equations 
for the modified BDF2 and  BDF3 schemes because the proofs for the remaining ones are similar.

We first consider the modified BDF2 scheme. Setting u^n = e^i(kx -  n ) 
and using the identities 
e^ix+e^-ix = 2cos x and cos2θ = 2cos^2θ -1 
in |1/2(3u^n-u^n-1)|^2, we get 

    | 3u^n-u^n-1/2|^2
    =   |1/2 e^ikx-i(n-1/2 ) (3e^-i/2 -e^i/2)|^2
    
    =   |e^ikx-i(n-1/2) (cos(/2) -2isin(/2))|^2 
    
    =   cos^2(/2) + 4sin^2(/2) 
    
    =   1/2(5 - 3cos()).

A similar calculation applies to |1/2(3u^n+1-u^n )|^2. Substituting 
u^n = e^i(kx -  n ) into the modified BDF2 scheme (<ref>)
and using (<ref>) yield 

    i/   (3u^n+1-4u^n+u^n-1/2)
    +3Δ u^n+1+ 2Δ u^n-Δ u^n-1/4
       =λ/2(|3u^n+1-u^n/2|^2 + |3u^n-u^n-12|^2)
    (3Δ u^n+1+ 2Δ u^n-Δ u^n-1/4), 
    
       =i/[3(e^-i-1) + (e^i-1)] 
       =1/2(k^2 + λ H_1(ω, ))
    [3 (e^-i+1) - (e^i+1)],

where H_1(ω,) =1/2(5-3 cos()).

It is easy to obtain the following equation from (<ref>):  

    2i/[3 (   e^-i/2 - e^i/2)e^-i/2
    +(e^i/2 - e^-i/2)e^i/2]
    
    =   3(k^2 + λ H_1(ω, ))(e^-i/2
    + e^i/2)e^-i/2
       -(k^2 + λ H_1(ω, ))(e^i/2
    +e^-i/2)e^i/2.

Using the identities e^ix-e^-ix = 2isin x and 
e^ix+e^-ix = 2cos x in (<ref>), we get 

    2/   sin(/2)
    (3e^-i/2 - e^i/2)
       =(k^2 + λ H_1(ω,  )) cos(/2)
    (3e^-i/2 - e^i/2).


Since 3e^-i/2 -e^i/2 0 in (<ref>), 
the numerical dispersion relation  of the modified BDF2 scheme can be written as

    2/sin(/2)
    =   (k^2 + λ H_1(ω, ))cos(/2)
    
    =   (k^2 +λ/2 (5 - 3cos()))cos(/2).

The desired equation (<ref>) holds by letting H(ω, ) = H_1(ω, ). 

Next, we consider the modified BDF3 scheme. Using the identities 
 e^ix+e^-ix = 2cos x and 
  cos2θ = 2cos^2θ -1 

we get

    |11u^n-7u^n-1 +2u^n-2/6|^2
    = |1/6 e^ikx-i(n-1)(11e^-i- 7 + 2e^i)|^2 
        1in = |1/6 e^ikx-i(n-1)(13cos()-
    7- 9isin())|^2 
        1in =1/18(44cos^2()-91cos^2() +65) 
        1in =1/36(83 - 47cos(2)).


Substituting u^n = e^i(kx -  n ) into the modified BDF3 scheme, we obtain 

    i/[
    11 (e^-i-1)+9(e^i-1)-2(e^2i-1)/6]
          =1/2(k^2 + λ H_2(ω, ))
    11 (e^-i-1)-5(e^i-1)+2(e^2i-1)-4/6,

where H_2(ω, ) = 1/36(83 - 47 cos(2)).

Similar to the derivation of (<ref>),    
using identities e^ix-e^-ix = 2isin x, e^ix+e^-ix = 2cos x, e^ix-e^-ix = 2isin x, 
and e^ix+e^-ix = 2cos x in (<ref>), we get 

    sin(/2)(11e^-i/2 - 9e^i/2)
    +2sin()e^i
            = cos(/2)(11e^-i/2 - 5e^i/2)
    +2cos()e^i - 2
       sin(/2)(11e^-i/2 - 9e^i/2
    + 4cos(/2)e^i) 
            =cos(/2)(11e^-i/2 - 9e^i/2
    + 4cos(/2)e^i),

where we have used the following identity in (<ref>):

    2cos()e^i-2
       =2(2cos^2(/2)-1)e^i
       =4cos^2(/2)e^i - 2(e^i+1)
       =4cos^2(/2)e^i - 4cos(/2)e^i/2.


Since (11e^-i/2 - 9e^i/2 + 4cos(/2)e^i) 0
in (<ref>), the numerical dispersion relation for the modified BDF3 can be written as

    2/sin(/2)
    =   (k^2 + λ H_2(ω,  ))
    cos(/2)
    
    =   (k^2 + λ/36(83 - 47cos(2)))
    cos(/2),

which gives (<ref>) after setting H(ω, ) = H_2(ω,). 
The proof is complete.





Table <ref> shows the computed dispersion errors and convergence rates 
for the modified BDF schemes. The numerical results indicate that these modified BDF schemes have a second 
order dispersion error.
 




§ NUMERICAL EXPERIMENTS: VALIDATING THE CONVERGENCE RATES.
 
In this section, we present several 1D numerical tests to illustrate our theoretical results, in particular, 
to verify the rates of convergence of the proposed time-stepping schemes. Our computations are done using the 
software package FEniCS and the linear finite element method is employed for the spatial discretization 
in all our numerical tests.

We consider the cubic nonlinear Schrödinger equation <cit.>
 (i.e., f(s) = s, λ = -2)

i u_t +Δu +2|u|^2 u    =0,           t>0,

u (0)    = u _0,        

where initial data  u_0 = (x)exp(2ix) is chosen so that the
exact solution is given by <cit.>

    u(x,t) = (x-4t)exp(i(2x-3t)).










We solve problem  (<ref>)–(<ref>) with the periodic boundary condition in [-20, 20]
by a few selected schemes from the family of  the time–stepping schemes  proposed in Section <ref>.  
In  Theorem <ref>, we proved the mass- 
and energy-conservation properties of scheme (<ref>)–(<ref>) 
without any restrictions on h and . 
To show the conservation properties of problem (<ref>)–(<ref>),
uniform spatial and temporal  meshes are used with h = 2^-5 and = 2^-6 
and the discrete mass and energy are defined in (<ref>).

The evolution of the L^2-norm of R^n and u^n for one trajectory of 
the modified Crank-Nicolson scheme and Leapfrog scheme are shown in Figure <ref>.
For the Crank-Nicolson scheme, we also show in Figure <ref> that the evolution  
of the L^2-norm is exactly conserved although at each time step, the nonlinear equation 
is not exactly solved. Meanwhile, the energy is also exactly conserved  as shown in Figure <ref>. 

It should be noted that the Leapfrog scheme is a multi-step method,  it requires two starting values, 
which are usually generated by a one-step method. 
In order to choose a suitable starting one-step scheme, we present computational results using different 
starting values generated by the exact solution and the Crank-Nicolson scheme.
The computed evolution of L^2-norm and the energy of  the Leapfrog scheme are shown in 
Figure <ref> and Figure <ref>, from which we observe that the  
evolution has a large oscillation at the early stage, and gradually becomes stable.

The evolution of the L^2-norm for one trajectory of the  modified BDF schemes 
are show in Figure <ref>,fig_bdf3_m_e:a and Figure <ref>,
fig_bdf5_m_e:a,fig_bdf6_m_e:a.
As expected, all of them gradually become stable and are exactly conserved. 
The differences of these oscillations (the zoom-in graphics at the beginning stage of the evolution) 
are obtained in Figure <ref>. We observe that all modified BDF schemes start with some 
oscillations and the amplitudes of the oscillations quickly diminish with time. 
In addition,  similar phenomenon for the energy are also seen in  
Figure <ref>,fig_bdf3_m_e:b, Figure <ref>,fig_bdf5_m_e:b,
fig_bdf6_m_e:b and Figure <ref>.


The accuracy of a proposed method is examined numerically by comparing the solution obtained on 
a sequence of coarse (time) meshes with the exact solution given in (<ref>).  
The computed errors 
and rates of the Crank-Nicolson scheme are shown in Table <ref>. We observe 
that the L^2-norm errors decrease by a factor 2 when the step-size  
is halved. Hence,  a second order convergence rate  is verified.  



Recall that the Leapfrog scheme is a multi-step numerical method. 
In order to choose a suitable starting (one-step) scheme,   
we present some convergence results using difference starting values generated by the exact 
solution and by the Crank-Nickson scheme in Table <ref>. 
It is clear to see that the L^2-norm error rate is 2, which confirms our theoretical result.



We run the same tests for the proposed modified BDF schemes, the computed results are shown  
in Table <ref>, we again observe that the L^2-norm error rate is 2.


 
We conclude this section by presenting a convergence and performance comparison of the 
Leapfrog scheme and the modified BDF schemes in Figure <ref>.
 





§ NUMERICAL EXPERIMENTS: CAPTURING THE BLOW-UP TIME.

Our aim in this section is to present a numerical study of the blow-up phenomenon for the 
quintic nonlinear Schrödinger equation, which is known to be very delicate to simulate  
in order to have an accurate understanding of this behavior. In particular, we focus 
on predicating and capturing the blow-up time using several proposed schemes. 

Existing numerical results have shown that starting with an initial condition of a given amplitude, 
one can claim that the solution has a singularity as soon as its amplitude becomes three (or more) 
times bigger than the initial amplitude <cit.>. 



We consider the following quintic nonlinear Schrödinger problem:

i u_t +Δu +|u|^4 u    =0,           t>0,  -10 ≤x ≤10,

u (0) = u    _0,        

with periodic boundary conditions <cit.>. The initial condition is chosen as 
u_0 = 1.6 e^-x^2. Since the initial energy is negative,  it is known that the
a blow-up in the solution must occur in finite time <cit.>.



 §.§ Comparison of -norm profiles of the computed solutions by different schemes.

In this subsection, we want to test whether all or which of our proposed time-stepping 
schemes will be able to capture the blow-up phenomenon.

Figure <ref> shows the simulation results of various schemes. 
We clearly see the formation of a singularity and that the used mesh size is small enough to
capture the essential feature of the blow-up by all but the Leapfrog scheme in 
Figure <ref>, As a result, we conclude that the Leapfrog scheme is not, 
but all other proposed schemes are, capable of capturing the blow-up phenomenon of the quintic 
nonlinear Schrödinger equation (<ref>). 



The further comparison of -norm profiles obtained by other schemes (excluding the Leapfrog scheme) 
are shown in Figure <ref>–<ref>. The tests in Figure <ref> 
indicate that the Crank-Nickson scheme and the modified BDF schemes are capable of 
capturing the blow-up phenomenon. 

In order to study whether the blow-up phenomenon will affect the mass- and 
energy-conservation results of Theorem <ref>, we present the time evolution 
of the mass and energy of R^n in Figure <ref>–<ref>. 
As expected, the mass of R^n is exactly conserved and the energy 
is also conserved before and after the blow-up time, in spite of a sharply 
increase in energy at the blow-up time. The similar behaviors of u^n are observed 
from Figure <ref>–<ref>.






 §.§ Comparison of –norm profiles obtained using different time-step sizes.

To better understand the sensitivity of the blow-up simulations to the time-step size,
we analyze the capability of different schemes for capturing the blow-up phenomenon 
and provide three criteria for the blow-up time in this subsection.

The comparison of –norm profiles obtained by different time-step sizes 
for the proposed time-stepping schemes are shown in Figure <ref>.
We observe similar behavior  for most schemes. 
Different simulation results for the Leapfrog scheme are shown in Figure <ref> 
and for the linearized scheme of <cit.> in Figure <ref>, 
although small enough time-step sizes are used.
As already mentioned earlier, our numerical tests show that the 
Crank-Nickson scheme (see Figure <ref>) and the modified BDF schemes 
(see Figure <ref>–fig_schemes_max:g) are capable of capturing the 
blow-up phenomenon. These results also consist with our previous results. 




In order to provide some criterion for identifying the blow-up time, we first propose three different 
criteria and then to present a comparison of them on different schemes in Table <ref>. 
The first criterion is to identify t^max corresponding to the time at which u^n_ takes 
its maximum as shown in Figure <ref>, where u^max represents the maximum value.  
One difficulty with this criterion is that we may get different t^max with different schemes. 
To be specific, the modified lower order BDF schemes (i.e., M–BDF2, M–BDF3 and M–BDF4) identify 
the same earliest t^max, while the modified higher order BDF schemes (i.e., M–BDF5 and M–BDF6) 
capture the same latest t^max. In addition, the t^max found by the Crank-Nickson scheme 
is in the middle of the above two values, so it is inconclusive that which t^max is the most accurate.

However, it should be noted that the Crank-Nickson scheme is a preferable scheme for capturing the blow-up 
phenomenon because it finds the largest u^max_ as shown in Table <ref>.

The other two criteria identify t_1^n and t_2^n
in Figure <ref>, where t_1^n represents the time  
point at which the energy of R^n is the smallest and 
t_2^n denotes the time point at which the energy of R^n has the maximum increase. 
From Table <ref> we observe that all time-stepping schemes identify the same blow-up time using 
these two criteria, which shows the robustness of both criteria. 






 



§ CONCLUSION.

In this paper we present a family of mass- and energy-conserved time-stepping 
schemes for general nonlinear Schrödinger equations. This includes the modified 
Crank-Nicolson scheme, Leapfrog scheme and modified BDF schemes as well  as 
a four-step symmetric scheme. We have shown that the proposed schemes have second-order 
convergence while preserving both mass and energy in the discrete setting without 
any mesh constraint. We also derive the dispersion relation equation for 
each of the proposed schemes and numerically show the convergence orders  
for the numerical dispersions. Extensive numerical experiments have been presented 
to illustrate the performance of the proposed schemes and to validate the theoretical 
results of the paper. Additional numerical experiments have also be provided 
to test the capability of the proposed schemes for capturing the blow-up phenomenon of the quintic 
nonlinear Schrödinger equation. Various criteria are proposed for identifying
the blow-up time and their effectiveness is also extensively examined. 

It is a bit disappointing that all proposed time-stepping schemes of this paper only have
second order accuracy and second order truncation errors. A very interesting question
is whether it is possible to improve these schemes into higher order schemes while 
still conserving both mass and (a modified) energy. Another challenging question
is whether it is possible to construct mass- and energy-conserved linear schemes (that is,  
only a linear problem needs to be solved at each time step). These 
open questions are worthy of further investigation and will be addressed in a further work. 

	
Acknowledgments. 
The work of the first author was partially supported by the NSF Grant 
DMS-1620168, and the work of the second author was partially supported by NSF Grant
DMS-1812666. Part of the third author's work was done during a recent visit to the 
University of Tennessee at Knoxville (UTK), the author would like to thank Department 
of Mathematics of UTK for the support and hospitality. 
The visit was financially supported by a scholarship from the author's home institution,
Northwestern Polytechnical University of China. 
    

abrv

plain
10

Antoine2013Computational
X. Antoine, W. Bao, and C. Besse, Computational methods for the
	dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations,
Computer Physics Communications, 184 (2013), pp. 2621–2633.

bao2013numerical
W. Bao, Q. Tang and Z. Xu,
Numerical methods and comparison for computing dark and bright solitons in the nonlinear Schrödinger equation,
Journal of Computational Physics,
235(2013), pp. 423–445.

bao2013Mathematical
W. Bao, Y. Cai, Mathematical theory and numerical methods for Bose-Einstein condensation, Kinetic & Related Models,
6(2013), pp. 1937-5093.
 
Bona2013Conservative
J. L. Bona, H. Chen, O. Karakashian, and Y. Xing, Conservative,
	discontinuous Galerkin-methods for the generalized Korteweg-de Vries
	equation, Mathematics of Computation, 82 (2013), pp. 1401–1432.

Bourgain1999Global
J. Bourgain, Global Solutions of Nonlinear Schrödinger
	Equations, vol. 46, American Mathematical Society, 1999.

Calogero1988Spectral
F. Calogero and A. Degasperis, Spectral Transform and Solitons:
	How to Solve and Investigate SNnlinear Evolution Equations,
Springer, New York, 1988.

Debussche2002Numerical
A. Debussche and L. D. Menza, Numerical simulation of focusing
	stochastic nonlinear Schrödinger equations, Physica D Nonlinear
Phenomena, 162 (2002), pp. 131–154.






Lax2010Integrals
P. D. Lax, Integrals of nonlinear equations of evolution and
	solitary waves, Communications on Pure and Applied Mathematics, 21 (2010),
pp. 467–490.

Lin2017Numerical
J. Lin, Y. Hong, L.-H. Kuo, and C.-S. Liu, Numerical simulation of
	3D nonlinear Schrödinger equations by using the localized method of
	approximate particular solutions, Engineering Analysis with Boundary
Elements, 78 (2017), pp. 20–25.

Liu2018On
H. Liu, Y. Huang, W. Lu, and N. Yi, On accuracy of the
	mass-preserving DG method to multi-dimensional Schrödinger equations,
IMA Journal of Numerical Analysis,  (2018).

Liu2016A
H. Liu and N. Yi, A Hamiltonian preserving discontinuous Galerkin
	method for the generalized Korteweg-de Vries equation, Journal of
Computational Physics, 321 (2016), pp. 776–796.

Lu2015Mass
W. Lu, Y. Huang, and H. Liu, Mass preserving discontinuous
	Galerkin methods for Schrödinger equations, Journal of Computational
Physics, 282 (2015), pp. 210–226.

Mclachlan2003Featured
R. Mclachlan, Featured review: Geometric numerical integration:
	Structure-preserving algorithms for ordinary differential equations, Siam
Review, 45 (2003), pp. 817–821.

Miura1968Korteweg
R. M. Miura, Korteweg-de Vries equation and generalizations.
	II. existence of conservation laws and constants of motion., Journal of
Mathematical Physics, 9 (1968), pp. 1204–1209.

Pelinovsky1996Nonlinear
D. E. Pelinovsky, V. V. Afanasjev, and Y. S. Kivshar, Nonlinear
	theory of oscillating, decaying, and collapsing solitons in the generalized
	nonlinear Schrödinger equation, Physical Review E Statistical Physics
Plasmas Fluids and Related Interdisciplinary Topics, 53 (1996), p. 1940.

Sch1996Traveling
H. W. Schürmann, Traveling-wave solutions of the cubic-quintic
	nonlinear Schrödinger equation, Physical Review E Stat Phys Plasmas
Fluids Relat Interdiscip Topics, 54 (1996), pp. 4312–4320.

Sheu2015Dispersion
T. W. Sheu and L. Lin, Dispersion relation equation preserving
	FDTD method for nonlinear cubic Schrödinger equation, Journal of
Computational Physics, 299 (2015), pp. 1–21.

Taghizadeh2011Exact
N. Taghizadeh, M. Mirzazadeh, and F. Farahrooz, Exact solutions of
	the nonlinear Schrödinger equation by the first integral method,
Journal of Mathematical Analysis and Applications, 374 (2011), pp. 549–553.

Christopher1993Dispersion
C. K. Tam and J. C. Webb, Dispersion-relation-preserving finite
	difference schemes for computational acoustics, Journal of Computational
Physics, 107 (1993), pp. 262–281.

Tao2006Nonlinear
T. Tao, Nonlinear Dispersive Equations: Local and Global
	Aalysis, American Mathematical Society, 2006.

Wang201New
J. Wang, A new error analysis of Crank–Nicolson Galerkin
	FEMs for a generalized nonlinear Schrödinger equation, Journal of
Scientific Computing, 60 (2014), pp. 390–407.

Xu2005Local
Y. Xu and C. W. Shu, Local discontinuous Galerkin methods for
	nonlinear Schrödinger equations, Journal of Computational Physics, 205
(2005), pp. 72–97.







Zabusky1965Interaction
N. J. Zabusky and M. D. Kruskal, Interaction of "solitons" in a
	collisionless plasma and the recurrence of initial states, Physical Review
Letters, 15 (1965), pp. 240–243.

Zwillinger1989Handbook
D. Zwillinger, Handbook of Differential Dquations, 3rd ed,
Boston, MA: Academic Press, 1997.


 






