Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Weak multi-symplectic reformulation and geometric numerical integration for the nonlinear Schrödinger equations with delta potentials

Weak multi-symplectic reformulation and geometric numerical integration for the nonlinear... Abstract We discuss multi-symplectic (MS) geometric numerical integration for a class of important models in condensed matter physics, i.e., nonlinear Schrödinger equations (NLSEs) with inclusions of finite delta potentials. Due to the inclusions, they cannot be reformulated as multi-symplectic Hamiltonian systems (MSHSs) such as normal NLSEs. With regard to integration, we propose a kind of weak MS reformulations and address some intrinsic local and global conservation laws. Based on theoretical suggestions, concatenating multi-symplectic Runge–Kutta (MSRK) methods are novelly constructed for these weak MSHSs, and some local and global conservation laws under MSRK discretizations are investigated analytically. Extensive numerical experimentation is then presented including cases of single- and double-delta potentials and cases of uniform and nonuniform gird points. Various numerical comparisons made between MS and non-MS schemes show us that the remarkable advantage of the former proposed here is the precise preservation of local and global normalization conservation laws. Other conservative properties, by using the two classes of schemes, can be preserved within some accuracies over long-time evolution. But, in general, the MS formulations behave more stably than the other in the numerics. 1. Introduction Nonlinear Schrödinger equation (NLSE) is one of the most important nonlinear models in modern science. It appears in quantum field theory, condensed matter, nonlinear optics, soliton waves, the theory of partial differential equations (PDEs), the theory of dynamical systems and so forth, in physics and mathematics (see, e.g., Berezin & Shubin, 1991; Haus & Wong, 1996; Dalfovo et al., 1999; Sulem & Sulem, 1999; Carretero-González et al., 2008). We focus on the model based on Bose–Einstein condensation, which was proposed theoretically in 1925 and realized experimentally in 1995 for rubidium gas. And it provides opportunities for exploring quantum phenomena on a macroscopic scale (Dalfovo et al., 1999; Pethick & Smith, 2002; Carretero-González et al., 2008). Let us briefly review the NLSE or Gross–Pitaevskii (GP) equation in condensed matter physics (Dalfovo et al., 1999; Pethick & Smith, 2002). We consider a sufficiently dilute ultracold atomic gas which consists of $$N$$ interacting bosons of mass $$m$$ and is confined by a time-independent external potential $$V_{{\rm ext}}({\bf r})$$. Then the many-body Hamiltonian, in terms of the second quantization, is given by   H^ =∫d3rΨ^†(r,t)[−ℏ22m∇2+Vext(r)]Ψ^(r,t) +12∫d3rd3r′Ψ^†(r,t)Ψ^†(r′,t)V(r−r′)Ψ^(r′,t)Ψ^(r,t), (1.1) where $$\hat{{\it{\Psi}}}({\bf r},t)$$ and $$\hat{{\it{\Psi}}}^\dagger({\bf r},t)$$ are the boson annihilation and creation field operators, respectively, $$V({\bf r}-{\bf r}')$$ is the two-body interatomic potential. The Heisenberg evolution equation for the field operator, i.e., $${\bf i}\hbar\, \partial_t\hat{{\it{\Psi}}} =[\hat{{\it{\Psi}}},\hat{H}]$$, gives rise to   iℏ∂tΨ^(r,t)=[−ℏ22m∇2+Vext(r)+∫d3r′Ψ^†(r′,t)V(r−r′)Ψ^(r′,t)]Ψ^(r,t). (1.2) Proceeding further, one writes   Ψ^(r,t)=Ψ(r,t)+Ψ^′(r,t), (1.3) where $${\it{\Psi}}({\bf r},t) \equiv \langle \hat{{\it{\Psi}}}({\bf r},t)\rangle$$ is the expectation value of the field operator and is commonly called the macroscopic wavefunction of the condensate. Below the critical temperature, the noncondensate part $$\hat{{\it{\Psi}}}'({\bf r},t)$$ is negligible and then the mean-field approximation, that is, $$\hat{{\it{\Psi}}}({\bf r},t) = {\it{\Psi}}({\bf r},t)$$ can be used. For a dilute ultracold gas at low energy, only binary collisions are relevant and can be characterized by the $$s$$-wave scattering length $$a$$, which is independent of the details of the two-body potential. Thus, the interatomic potential can be replaced by an effective delta-function interaction potential,   V(r−r′)=gδ(r−r′), (1.4) with the coupling constant $$g=4\pi \hbar^2 a/m$$. The substitutions of the effective delta-potential and the mean-field approximation into (1.1) lead to the following nonlinear Schrödinger equation or Gross–Pitaevskii equation   iℏ∂tΨ(r,t)=[−ℏ22m∇2+Vext(r)+g|Ψ(r,t)|2]Ψ(r,t). (1.5) As pointed out in Dalfovo et al. (1999), the validity of (1.5) is based on the following: (I) the $$s$$-wave scattering length $$a$$ should be much smaller than the average distance between atoms; hence the effective delta potential (1.4) can be used to replace the interatomic potential; (II) the number of atoms in the condensate should be much larger than 1; then the classical mean field given in (1.3) can be used to substitute the field operator approximately. For time-independent external potentials, the NLSE or GP model admits two important conserved quantities: (C1) The boson gas is only investigated nonrelativistically, and no real particles can be created or annihilated. Thus the total number of atoms (with a suitable normalization)   N=∫|Ψ(r,t)|2dr (1.6) should be conserved in the dynamics; (C2) The gas system is an isolated system, and thus the total energy   E=∫[ℏ22m|∇Ψ(r,t)|2+Vext|Ψ(r,t)|2+12g|Ψ(r,t)|4]dr (1.7) should also be conserved. We consider the $$1+1$$-dimensional NLSE for the macroscopic wave function   iℏ∂tΨ=−ℏ22m∂xxΨ+g|Ψ|2Ψ+VextΨ, (1.8) where $$ V_{{\rm ext}}=V_{{\rm ext}}(x)$$ is the external potential. For the BEC species (Dalfovo et al., 1999), $$a$$ may take either positive or negative values depending on repulsive or attractive interactions between the particles, respectively. With the help of the transformations,   t~=g2t,x~=|g|xandϕ(x~,t~)=Ψ(x,t)|g|=Ψ(x~|g|,t~g2)|g|, the parameter $$g$$ can be removed up to a sign. We will use atomic units $$m=\hbar=1$$ throughout the whole text and fix $$g=\pm 1$$ corresponding to the case of a repulsive and an attractive nonlinearity, respectively, except other specified. With regard to the case of $$V_{{\rm ext}}\equiv 0$$, there has been increasing interest not only on analytically solving this problem and its variants in the physical and mathematical areas (see, e.g., Haus & Wong, 1996; Kivshar & Luther-Davies, 1998; Witthaut et al., 2005, 2009; Carretero-González et al., 2008; Rapedius & Korsch, 2009) but also on various numerical treatments including symplectic or multi-symplectic methods (see, e.g., Reich, 2000; Islas et al., 2001; Bridges & Reich, 2006; Hong et al., 2007; Kong et al., 2015 and references therein). To reformulate this problem in the form of a multi-symplectic Hamiltonian system (MSHS), by separating the real and imaginary parts of the wave function $${\it{\Psi}}=q+{\bf i}p$$ and introducing a pair of conjugate momenta, $$v=\partial_x q$$ and $$w=\partial_x p$$, one obtains   M∂tz+K∂xz=∇zS(z), (1.9) where $$z=(q,p,v,w)^T$$ is the state variable; $$M$$ and $$K$$ are skew-symmetric matrices   M =(0−200200000000000),K=(00100001−10000−100), and $$S: \mathbf{R}^4\rightarrow \mathbf{R}$$ is the multi-symplectic Hamiltonian $$S(z)=\frac{1}{2}g(q^2+p^2)^2-\frac{1}{2}(v^2+w^2).$$ We call a system of the form (1.9) an MSHS, because it admits the following multi-symplectic conservation law (MSCL) (Bridges, 1997; Reich, 2000; Bridges & Reich, 2006):   ∂tω+∂xκ=0, (1.10) where $$\omega$$ and $$\kappa$$ are the pre-symplectic forms defined by   ω=12dz∧Mdzandκ=12dz∧Kdz. (1.11) Besides the MSCL, every MSHS has other two local conservation laws, i.e., (a) the energy conservation law (ECL)   ∂tE+∂xF=0, (1.12) with the energy density $$E=S(z)-\frac{1}{2}z^TK\partial_x z$$ and the energy flux density $$F=\frac{1}{2}z^TK\partial_t z$$ and (b) the momentum conservation law (MCL)   ∂tI+∂xG=0, (1.13) with the momentum density $$I=\frac{1}{2}z^TM\partial_x z$$ and the momentum flux density $$G=S(z)-\frac{1}{2}z^TM\partial_t z$$. The MSHS (1.9), originating from the nonlinear Schrödinger equation, has a more conservate law and it is referred to as the charge conservation law, or the particle-number conservation law, the mass conservation law in various physical problems:   ∂tNd+∂xNf=0, (1.14) where the density $$N_d=\frac{1}{4}(Mz)^TMz$$ and the flux density $$N_f=\frac{1}{2}(Mz)^TKz$$. No matter what conservation law is called, it commonly describes some physical quantity of finite value and can be normalized to unity. Hereafter we call it the local normalization conservation law (NCL). Under appropriate boundary conditions, e.g., zero or periodic boundary conditions, the local conservations laws (1.10), (1.12), (1.13) and (1.14) imply the corresponding global ones, which are obtained by integrating the local ones over the spatial domain and then eliminating the flux terms by suitably imposing boundary conditions (see, e.g., Bridges & Reich, 2006; Hong & Li, 2006): (I) The global symplecticity in time   ddt∫ω(z(x,t))dx=0. (1.15) (II) The global ECL   ddtE≜ddt∫E(z(x,t))dx=0. (1.16) (III) The global MCL   ddtI≜ddt∫I(z(x,t))dx=0. (1.17) (IV) The global NCL   ddtN≜ddt∫Nd(z(x,t))dx=0. (1.18) We should point out here that $$\mathscr{N}$$ appears in either (1.6) or (1.18), which is essentially the same. The expressions of $$\mathscr{E}$$ appearing in (1.7) and (1.16), as shown in Proposition 2.10, are equivalent under suitable boundary conditions. Concerning the case $$V_{{\rm ext}}\not\equiv 0$$ for (1.8), it does not seem to be very optimistic from an MS geometric viewpoint. To the best of our knowledge, published work about the MSHS reformulation of (1.8) and relevant MS geometric numerical methods applied to this model in the literature are rare especially for the case $$V_{{\rm ext}}$$ explicitly depending on $$x$$ or $$t$$. In this article, we focus on the equation (1.8) with finite external delta potentials $$V_{{\rm ext}}=\sum\limits_{\mu=1}^{J}\gamma_{\mu}\delta(x-x_{\mu}^{*})$$, where $$J$$ is a positive integer and $$x_{\mu}^{*}$$ ($$\mu=1,\cdots,J$$) are all given points with the assumption that $$x_1^{*}<x_2^{*}<\cdots<x_J^{*}$$. More precisely, the following NLSE   i∂tΨ=−12∂xxΨ+g|Ψ|2Ψ+[∑μ=1Jγμδ(x−xμ∗)]Ψ, (1.19) will be investigated theoretically and numerically from a viewpoint of the MSHS theory. There exist a lot of references investigated on stationary and dynamical delta-potential NLSEs, which have single, double, finite delta and delta-shell potentials, and even infinite ones like the Dirac–Comb potentials (Berezin & Shubin, 1991; Haus & Wong, 1996; Kivshar & Luther-Davies, 1998; Witthaut et al., 2005, 2009; Rapedius & Korsch, 2009; Cartarius et al., 2012). These are imposed by various kinds of boundary conditions, e.g., periodic or zero at infinite and finite distances, absorbing and reflecting, Siegert ones and many others arising in different physical problems. In this context, we just restrict to the first class, namely the infinitely distant periodic or zero ones, which are commonly related to bound-state problems, while others to the resonance state or scattering problems. Due to the presence of the delta potential, the state variables $$z$$ admit jumps and (1.19) cannot, be reformulated as an MSHS in the normal sense. Noting that the continuity equations (e.g., the conservation of mass, the conservation of momentum and the conservation of energy) in hydrodynamics (Landau & Lifshitz, 1987), are valid in the sense of integration and then rewritten in the form of differentiation we make an attempt to reformulate the above equation into an MSHS form in the sense of integration in Section 2. Based on the given weak MSHS form, relevant theoretical discussions including some local and global conservation laws are addressed as well. In Section 3, novel concatenating RK discretizations for (1.19) are constructed, the multi-symplecticity of them are revealed, and some local and global intrinsic conservation laws are investigated theoretically under multi-symplectic Runge–Kutta (MSRK) discretizations. Section 4 is devoted to numerical experiments and extensive relevant discussions. We conclude in Section 5. 2. The weak MSHS reformulation and conservation laws For a function $$f(x,t)$$ which is well defined in some deleted neighborhood of $$(\hat{x},t)$$, we define   fx=x^−=limx→x^−f(x,t)andfx=x^+=limx→x^+f(x,t). In what follows, we simplify $$f^-_{x=x_{\mu}^{*}}=f^-_\mu$$ and $$f^+_{x=x_{\mu}^{*}}=f^+_\mu$$ if there is no confusion. According to the distribution theory (see, Witthaut et al., 2005; Carretero-González et al., 2008 and references therein), the wave function $${\it{\Psi}}$$ in (1.19) is continuous, but its first spatial partial derivative admits the following discontinuities (or jumps) at $$x = x_{\mu}^{*}$$ that   ∂xΨμ+−∂xΨμ−=2γμΨ(xμ∗), because of the presence of the corresponding delta potentials, where $$\mu=1,2,\ldots,J$$. Hereafter, the points $$x= x_{\mu}^{*}$$, as proposed in Bai & Wang (2015) and Bai (2016), are referred to as interface points. Proposition 2.1 The equation (1.19), by introducing the state variable $$z=(q,p,v,w)^T$$ as specified in (1.9), is equivalent to the following system, in which the singular potentials $$\gamma_\mu\delta (x-x_{\mu}^{*})$$ are replaced by some extra jump conditions:   {M∂tz+K∂xz=∇zS(z),x≠xμ∗zμ+−zμ−=2γμ(0,0,qμ−,pμ−)T≜Ajump(μ)zμ−,  (2.1) for $$\mu=1,\ldots,J$$, where $$z^{\pm}_\mu=(q^{\pm}_\mu,p^{\pm}_\mu,v^{\pm}_\mu,w^{\pm}_\mu)^T$$ and the matrix   Ajump(μ)=(000000002γμ00002γμ00). Furthermore, it can be verified from (2.1) that the derivatives of $$v$$ and $$w$$, with respect to $$x$$, namely the corresponding second derivatives of $$q$$ and $$p$$, respectively, are continuous at the interface point $$x=x_{\mu}^{*}$$, i.e.,   ∂xvμ+=∂xvμ−and∂xwμ+=∂xwμ−. (2.2) It is apparent that the system (2.1), at any point $$x\neq x_{\mu}^{*}$$, is equivalent to the MSHS (1.9), and the local conservation laws (1.10), (1.12)–(1.14) are then automatically valid. Thus the local conservation laws are investigated only at the interface points subsequently. 2.1 Local conservation laws Noting the fact that the exterior differential operator $$d$$ can commute with the partial derivatives $$\partial_t$$ and $$\partial_x$$, we first give the variational equations in terms of (2.1)   {M∂tdz+K∂xdz=DzzS(z)dz,x≠xμ∗,dzμ+−dzμ−=Ajump(μ)dzμ−,  (2.3) for $$\mu=1,\ldots,J$$. Remark 2.2 The local linearizations, which are namely the differential 1-forms from the viewpoint of differential geometry (Bridges, 1997; Arnold, 1999), should have the same properties as the original state variables possess. Illustratively speaking, $$dq^+_\mu=dq^-_\mu$$ come from the fact that the local linearizations $$dq_\mu$$ should be continuous at the interface points $$x=x_{\mu}^{*}$$ as $$q$$ admits, and so $$dp^+_\mu=dp^-_\mu$$; $$dv^+_\mu=dv^-_\mu + 2\gamma dq^-_\mu$$ originate from the fact that $$dv$$ should have the same jumps at $$x=x_{\mu}^{*}$$ as $$v$$ has, and $$dw^+_\mu=dw^-_\mu +2\gamma dp^-_\mu$$. It is worth revisiting the deduction of the well-known Euler–Lagrange equations (Arnold, 1999; Hairer et al., 2006); the variation of the solution, in the action functional, is set to be vanished at endpoints due to the fixed boundary values. Thus it is necessary and natural to combine the discontinuities into variational demonstrations here. To be distinguished, the phrases ‘in a weak sense’ or ‘in the weak sense’ will be used in the following assertions. We say a conservation law $$\partial_t C_d+\partial_x C_f=0$$ holds in the weak sense at some point $$(\hat{x},\hat{t})$$, it means for any small enough rectangle space-time interval $$[x_1,x_2]\times[t_1,t_2]$$ with $$\hat{x}\in (x_1,x_2)$$ and $$\hat{t} \in (t_1,t_2)$$ such that the following integration   ∫x1x2(Cd(x,t2)−Cd(x,t1))dx+∫t1t2(Cf(x2,t)−Cf(x1,t))dt=0, (2.4) holds. We always check, in the weak sense, the validity of a conservation law through the above relation. Theorem 2.3 The system (2.1), in the weak sense we clarified above, admits the local MSCL (1.10) at each interface point $$x=x_{\mu}^{*}$$. Proof. Given any small enough rectangle interval $$[x_1,x_2]\times [t_1,t_2]$$ with $$x_\mu^{*}\in (x_1,x_2)$$ and excluding any other interface point in $$[x_1,x_2]$$, we will check the validity of the MSCL (1.10) in the weak sense. As suggested in (2.4), we check    ∫x1x2[ω(z(x,t2))−ω(z(x,t1))]dx+∫t1t2[κ(z(x2,t))−κ(z(x1,t))]dt =∫x1xμ∗−ϵ[ω(z(x,t2))−ω(z(x,t1))]dx+∫t1t2[κ(z(xμ∗−ϵ,t))−κ(z(x1,t))]dt +∫xμ∗−ϵxμ∗+η[ω(z(x,t2))−ω(z(x,t1))]dx+∫t1t2[κ(z(xμ∗+η,t))−κ(z(xμ∗−ϵ,t))]dt +∫xμ∗+ηx2[ω(z(x,t2))−ω(z(x,t1))]dx+∫t1t2[κ(z(x2,t))−κ(z(xμ∗+η,t))]dt, where $$\epsilon$$ and $$\eta$$ are small enough positive real numbers. Noting the fact that the MSCL (1.10) is valid in the normal sense at any point $$x\neq x_\mu^{*}$$, the first and third lines of the right-hand side (RHS) of the above are automatically vanishing, i.e.,    ∫x1x2[ω(z(x,t2))−ω(z(x,t1))]dx+∫t1t2[κ(z(x2,t))−κ(z(x1,t))]dt =2∫xμ∗−ϵxμ∗+η(dp∧dq)|t1t2dx+∫t1t2(dq∧dv+dp∧dw)|xμ∗−ϵxμ∗+ηdt. (2.5) Let us consider the limiting case as $$\epsilon$$ and $$\eta$$ tend to zero positively. With the aid of the jumps given in the second equation of (2.3), the first term in the RHS of (2.5) vanishes because of the continuity of $$dp$$ and $$dq$$ at the interface point $$x=x_\mu^{*}$$, and the second term becomes    ∫t1t2(dqμ+∧dvμ++dpμ+∧dwμ+−dqμ−∧dvμ−−dpμ−∧dwμ−)dt =∫t1t2(dqμ−∧(dvμ+−dvμ−)+dpμ−∧(dwμ+−dwμ−))dt =∫t1t2(dqμ−∧(2γμdqμ−)+dpμ−∧(2γμdpμ−))dt=0, (2.6) where the first equality uses the fact that $$dq^+_\mu =dq^-_\mu $$ and $$dp^{+}_\mu =dp^{-}_\mu$$. Thus (2.5) and (2.6) give rise to this theorem. □ Definition 2.4 For the purpose of being distinguished, a system of the abstract form (1.9) is called a normal MSHS, and a system of the form (2.1) is called a weak MSHS since it preserves the MSCL (1.10) in the weak sense at each interface point $$x=x_\mu^{*}$$. And conservation laws, which are proved valid in the weak sense, are also referred to with the term ‘weak’. Theorem 2.5 The system (2.1) admits the weak local NCL (1.14) at each interface point $$x=x_\mu^{*}$$. Proof. As shown in Theorem 2.3, it suffices to prove that   ∫xμ∗−ϵxμ∗+η[Nd(z(x,t2))−Nd(z(x,t1))]dx+∫t1t2[Nf(z(xμ∗+η,t))−Nf(z(xμ∗−ϵ,t))]dt=0, i.e.,   ∫xμ∗−ϵxμ∗+η(p2+q2)|t1t2dx+∫t1t2(qw−pv)|xμ∗−ϵxμ∗+ηdt=0 (2.7) with the positive parameters $$\epsilon$$ and $$\eta$$ tending to zero. The left-hand side (LHS) of 2.7, as $$\epsilon$$ and $$\eta$$ tend to zero positively, yields    ∫t1t2(qμ+wμ+−pμ+vμ+−qμ−wμ−+pμ−vμ−)dt=∫t1t2(qμ−(2γμpμ−)−pμ−(2γμqμ−))dt=0. (2.8) This completes the proof. □ Theorem 2.6 The system (2.1) admits the weak local ECL (1.12) at each interface point $$x=x_\mu^{*}$$. Proof. As discussed in Theorems 2.3 and 2.5, with the help of (2.1), we directly implement the following calculations    ∫xμ∗−ϵxμ∗+η[E(z(x,t2))−E(z(x,t1))]dx+∫t1t2[F(z(xμ∗+η,t))−F(z(xμ∗−ϵ,t))]dt =12∫xμ∗−ϵxμ∗+η(g(p2+q2)2−(q∂xv+p∂xw))|t1t2dx+12∫t1t2(q∂tv+p∂tw−v∂tq−w∂tp)|xμ∗−ϵxμ∗+ηdt. Let the positive parameters $$\epsilon$$ and $$\eta$$ tend to zero the first term of the right-hand side (RHS) of the above vanishes because of the fact that $$p$$ and $$q$$ are continuous at $$x=x_\mu^{*}$$ and also the same property for $$\partial_x\! v$$ and $$\partial_x\! w$$ as specified in (2.2). Also with $$\partial_t q^+_\mu=\partial_t q^-_\mu$$, $$\partial_t p^+_\mu=\partial_t p^-_\mu$$, $$\partial_t v^+_\mu-\partial_t v^-_\mu=2\gamma_{\mu}\partial_t q_\mu^-$$ and $$\partial_t w^+_\mu-\partial_t w^-_\mu=2\gamma_{\mu}\partial_t p_\mu^-$$, the second term reads   12∫t1t2(qμ−(∂tvμ+−∂tvμ−)+pμ−(∂twμ+−∂twμ−)+∂tqμ−(vμ−−vμ+)+∂tpμ−(wμ−−wμ+))dt =12∫t1t2(qμ−(2γμ∂tqμ−)+pμ−(2γμ∂tpμ−)−∂tqμ−(2γμqμ−)−∂tpμ−(2γμpμ−))dt=0. (2.9) Hence it completes the proof. □ Remark 2.7 The system (2.1), however, does not admit the weak local MCL (1.13) at each interface point. According to the Nöther theorem (Olver, 1993; Arnold, 1999), a continuous conservation law corresponds to a continuous variational symmetry. The time translation invariance implies the local ECL, and the space translation invariance implies the local MCL. Comparing (1.19) with (2.1), the delta potentials do not break the time translation invariance, and thus the local ECL is still valid. The space translation invariance, however, is broken at each interface point due to the extra jump conditions. 2.2 Global conservation laws It is evident that the local conservation laws are independent of whatever initial-boundary value conditions imposed. But it changes for the corresponding global ones. And the well posedness of the considered NLSEs need suitable initial and boundary conditions either. Throughout this context, we assume that the solution is smooth enough, except for the interface points $$x=x_\mu^{*}$$ for $$\mu=1,\ldots,J$$, the spatial interval we investigate is $$[x_L, x_R]$$ with all $$x_\mu^{*}\in (x_L, x_R)$$. Here $$x_L$$ and $$x_R$$ could be infinities in principle. To proceed, we supplement the system (2.1) by prescribing suitable conditions with:   {z|t=0=z0(x)≡(q0(x),p0(x),v0(x),w0(x))T,z(xL,t) = z(xR,t),∂xz(xL,t)=∂xz(xR,t),  (2.10) where $$(q_0(x),p_0(x), v_0(x), w_0(x))$$ satisfy the jump conditions as given in (2.1) and $$t\geq 0$$ for the second and third equalities. Theorem 2.8 Under the conditions prescribed in (2.10), the system (2.1) has the weak global symplecticity in time (1.15) with the integration interval $$[x_L, x_R]$$. Proof. As specified in Theorem 2.3, the local MSCL (1.10) is valid at each interface point $$x=x_\mu^{*}$$ in the weak sense, i.e., in the sense of integration. At any other point, the local MSCL is normally valid. Thus (1.10) can be taken as the integration over the spatial interval $$[x_L,x_R]$$ for both its left- and right-hand sides, i.e.,   ∫xLxR(∂ω∂t+∂κ∂x)dx=0. (2.11) As discussed for (2.4) and in Theorem 2.3, the second term in the LHS becomes   ∫xLxR∂κ∂xdx=κ(z(xR,t))−κ(z(xL,t))=0, (2.12) where the second equality is due to the periodic boundary conditions given in (2.10). Thus we have   ∫xLxR∂ω∂tdx=ddt∫xLxRωdx=0, where the first equality is due to the fact that the integration can commute with the temporal derivative. Hence the global symplecticity is conserved for all times,   ∫xLxRω(z(x,t))dx=∫xLxRω(z(x,0))dx=∫xLxRω(z0(x))dx,∀t≥0. (2.13) □ Similarly, we have the following statement about the the global ECL and NCL. Theorem 2.9 Under the assumptions of Theorem 2.8, the system (2.1) admits the weak global ECL (1.16) and the weak global NCL (1.18). As mentioned previously, the expressions of the total energy $$\mathscr{E}$$ given in (1.7) and (1.16) are not the same. To avoid the ambiguity, we give the following proposition. Proposition 2.10 Under the assumptions of Theorem 2.8, the total energies $$\mathscr{E}$$ appearing in (1.7) and (1.16), for the NLSE with the external delta potentials, are equivalent. Proof. By introducing the state variable $$z$$, the total energy appearing in (1.7), which is denoted by $$\mathscr{E}_1$$ to avoid confusion, reads   E1=∫xLxR[12(v2+w2)+∑μ=1Jγμδ(x−xμ∗)(q2+p2)+g2(q2+p2)2]dx. (2.14) With the aid of properties of the delta functions, we immediately get   E1 =[∫xLx1∗−ϵ1+∑μ=1J−1(∫xμ∗−ϵμxμ∗+ημ+∫xμ∗+ημxμ+1∗−ϵμ+1)+∫xJ∗−ϵJxJ∗+ηJ+∫xJ∗+ηJxR]12(v2+w2)dx +∑μ=1Jγμ[q2(xμ∗)+p2(xμ∗)]+g2∫xLxR(q2+p2)2dx, (2.15) where we suppress the time variable $$t$$, i.e., $$q(x_\mu^*)\equiv q(x_\mu^*,t)$$ and $$p(x_\mu^*)\equiv p(x_\mu^*,t)$$, $$\epsilon_\mu$$ and $$\eta_\mu$$ are small enough positive real numbers. Let $$\epsilon_\mu$$ and $$\eta_\mu$$ tend to zero, take the integration by parts and recall the continuities of $$q$$, $$p$$, $$\partial_x v$$ and $$\partial_x w$$ in the definite domain, we find the first line in the RHS of (2.15) yields    −12∫xLxR(q∂xv+p∂xw)dx−12∑μ=1J{qμ−[vμ+−vμ−]+pμ−[wμ+−wμ−]} =−12∫xLxR(q∂xv+p∂xw)dx−∑μ=1Jγμ[q2(xμ∗)+p2(xμ∗)]. (2.16) Combining (2.15) with (2.16) gives $$\mathscr{E}_1=\mathscr{E}$$. □ Hereafter, we only use the expression of the total energy as given in (1.16) while encountering. 3. Multi-symplectic Runge–Kutta methods In this section, we will introduce novel concatenating RK discretizations applied to the weak MSHS (2.1) and discuss the multi-symplecticity of the numerical methods. To simplify the presentation, we just limit our discussion to the case on a single external delta potential, i.e., $$J=1$$ in (1.19). We remark here that, it can be naturally extended to cases on finite ones. For convenience, we introduce some notions of MS schemes applied to a normal MSHS of the abstract form (1.9). Given a partition over the definite domain, one considers a certain discretization for an MSHS (1.9)   M∂ti,lzil+K∂xi,lzil=(∇zS(z))il, (3.1) at a mesh point $$(x_i,t_l)\in R^2$$, where $$\partial_t^{i,l}$$, $$\partial_x^{i,l}$$ are the corresponding discretizations of the partial derivatives $$\partial_t$$ and $$\partial_x$$, respectively, $$(\nabla_z S(z))_i^l=(\nabla_z S(z_i^l))$$, $$z_i^l$$ is the numerical value of $$z(x,t)$$. Definition 3.1 (Reich, 2000; Bridges & Reich, 2006) The numerical scheme (3.1) of (1.9) is said to be multi-symplectic if   ∂ti,lωil+∂xi,lκil=0, (3.2) still holds at any mesh point $$(x_i,t_l)$$, where   ωil=12(dz)il∧M(dz)ilandκil=12(dz)il∧K(dz)il. As far as we have known, the MS concatenating RK discretizations originated from Reich (2000), and who showed concatenating Gauss–Legendre schemes applied to a $$1+1$$-dimensional MSHS that leads to an MS integrator. Since then there are quite a few works about MSRK-type methods (see, e.g., Hong et al., 2005; Ryland & McLachlan, 2008; McLachlan et al., 2014 and references therein) and their applications to various special MSHSs (see, e.g., Bridges & Reich, 2006; Hong & Li, 2006; Hong et al., 2007 and references therein). Almost all the works, however, are restricted to normal MSHSs. Regarding the weak MSHS (2.1) reformulated here, we make an attempt to introduce the corresponding MS concatenating RK discretizations. In our numeric, the spatial and temporal intervals are taken to be $$[x_L,x_R]$$ and $$[T_{0},T_{1}]$$, respectively. Here and hereafter, we will suppress the index $$1$$, e.g., $$f^{\pm}=f^{\pm}_{1}$$, $$A_{\rm jump}=A_{\rm jump}^{(1)}$$ and so forth. As proposed for the normal MS numerical treatments, we begin with a uniform partition on the definite domain with the mesh points $$(x_i,t_l)$$ and the stepsizes $${\it{\Delta}} x=\frac{x_R-x_L}{N}$$ and $${\it{\Delta}} t=\frac{T_1-T_0}{\hat{N}}$$, where $$x_i=x_L+i {\it{\Delta}} x$$ for $$i=0,1,\ldots,N$$, $$t_l=T_0+l {\it{\Delta}} t$$ for $$l=0,1,\ldots,\hat{N}$$. Assume that the interface point $$x=x^{*}$$ lies in the subinterval $$[x_I,x_{I+1})$$ for some $$I$$ ($$0<I<N$$), we then supplement the interface point as a new grid point labeled by $$x_{I^{'}}$$. Thus the spatial grid points become $$x_L=x_0<x_1<\cdots<x_I\leq x_{I^{'}}=x^{*}<x_{I+1}<\cdots<x_N=x_R$$, and the mesh-lengths of the two new subintervals are denoted by $${\it{\Delta}} x_1$$ ($$\equiv x_{I^{'}}-x_I$$) and $${\it{\Delta}} x_2$$ ($$\equiv x_{I+1}-x_{I^{'}}$$), respectively. In case $${\it{\Delta}} x_1=0$$, which means that the interface point $$x=x^{*}$$ is a grid point after the first equidistant partition, i.e., the grid point $$x_{I}$$ coincides with $$x_{I^{'}}$$ and no further partition is needed. We are now in a position to perform numerical discretization, which is demonstrated by applying $$s$$-stage and $$r$$-stage RK methods to the system (2.1) to the $$x$$- and $$t$$-directions, respectively, on space-time grids. For a starting point $$(x_i,t_l)$$ with $$i\notin \{I,I^{'}\}$$, we have   {Zi,ml,k=zil,k+Δx∑n=1sa~mn∂xZi,nl,k,zi+1l,k=zil,k+Δx∑m=1sb~m∂xZi,ml,k,Zi,ml,k=zi,ml+Δt∑j=1rakj∂tZi,ml,j,zi,ml+1=zi,ml+Δt∑k=1rbk∂tZi,ml,k,M∂tZi,ml,k+K∂xZi,ml,k=∇zS(Zi,ml,k).  (3.3) Here we will specify some notations which first appeared: $$Z_{i,m}^{l,k}\approx z(x_i+\tilde{c}_m{\it{\Delta}} x,t_l+c_k{\it{\Delta}} t)$$, $$z_{i,m}^{l}\approx z(x_i+\tilde{c}_m{\it{\Delta}} x,t_{l})$$, $$z_{i,m}^{l+1}\approx z(x_i+\tilde{c}_m{\it{\Delta}} x,t_{l+1})$$, $$\partial_t Z_{i,m}^{l,k}\approx \partial_t z(x_i+\tilde{c}_m{\it{\Delta}} x,t_l+c_k{\it{\Delta}} t)$$, $$\partial_x Z_{i,m}^{l,k}\approx \partial_x z(x_i+\tilde{c}_m{\it{\Delta}} x,t_l+c_k{\it{\Delta}} t)$$, $$z_{i}^{l,k}\approx z(x_{i},t_l+c_k{\it{\Delta}} t)$$ and $$z_{i+1}^{l,k}\approx z(x_{i+1},t_l+c_k{\it{\Delta}} t)$$ with $$\tilde{c}_m=\sum\limits_{n=1}^s \tilde{a}_{mn}$$ and $$c_k=\sum\limits_{j=1}^r a_{kj}$$ as required (see, e.g., Hairer et al., 2006). And for the starting points $$(x_I,t_l)$$ and $$(x_{I^{'}},t_l)$$, the schemes read, respectively,   {ZI,ml,k=zIl,k+Δx1∑n=1sa~mn∂xZI,nl,k,zI ′,−l,k=zIl,k+Δx1∑m=1sb~m∂xZI,ml,k,ZI,ml,k=zI,ml+Δt∑j=1rakj∂tZI,ml,j,zI,ml+1=zI,ml+Δt∑k=1rbk∂tZI,ml,k,M∂tZI,ml,k+K∂xZI,ml,k=∇zS(ZI,ml,k)  (3.3B) and   {ZI ′,ml,k=zI ′,+l,k+Δx2∑n=1sa~mn∂xZI ′,nl,k,zI+1l,k=zI ′,+l,k+Δx2∑m=1sb~m∂xZI ′,ml,k,ZI ′,ml,k=zI ′,ml+Δt∑j=1rakj∂tZI ′,ml,j,zI ′,ml+1=zI ′,ml+Δt∑k=1rbk∂tZI ′,ml,k,M∂tZI ′,ml,k+K∂xZI ′,ml,k=∇zS(ZI ′,ml,k),  (3.3C) where $$Z_{I,m}^{l,k}\approx z(x_{I}+\tilde{c}_m {\it{\Delta}} x_1,t_l+c_k{\it{\Delta}} t)$$, $$Z_{I^{'}\!\!,m}^{l,k}\approx z(x_{I^{'}}+\tilde{c}_m {\it{\Delta}} x_2,t_l+c_k{\it{\Delta}} t)$$, $$z_{I^{'}\!\!,\pm}^{l,k}\approx z^{\pm}(x_{I^{'}},t_l+c_k{\it{\Delta}} t)$$, etc. And the jump conditions at $$x_{I^{'}}$$ are subjected to   zI ′,+l,k−zI ′,−l,k=AjumpzI ′,−l,k. (3.4) The variational equations corresponding to (3.3) are given by   {dZi,ml,k=dzil,k+Δx∑n=1sa~mnd(∂xZ)i,nl,k,dzi+1l,k=dzil,k+Δx∑m=1sb~md(∂xZ)i,ml,k,dZi,ml,k=dzi,ml+Δt∑j=1rakjd(∂tZ)i,ml,j,dzi,ml+1=dzi,ml+Δt∑k=1rbkd(∂tZ)i,ml,k,Md(∂tZ)i,ml,k+Kd(∂xZ)i,ml,k=DzzS(Zi,ml,k)dZi,ml,k,  (3.5) where $$i\notin \{I,I^{'}\}$$ and the abbreviations, i.e., $$dZ_{i,m}^{l,k}$$ for $$(dZ)_{i,m}^{l,k}$$, $$d(\partial_t Z)_{i,m}^{l,j}$$ for $$(d(\partial_t Z))_{i,m}^{l,j}$$ and so on, are utilized for the sake of simplicity. The variational equations for ($$3.3$$B) and ($$3.3$$C) are similar to the above, and we would not exhibit them here anymore. And the jump conditions for the variational equations corresponding to (3.4), yield   dzI ′,+l,k−dzI ′,−l,k=AjumpdzI ′,−l,k. (3.6) Subsequently, we mention the numerical algorithm (3.3) which means (3.3) together with ($$3.3$$B) and ($$3.3$$C). With the above settings, we have the following theorem on the multi-symplecticity of the concatenating scheme. Since the deductions are similar to that given in Reich (2000) and Hong et al. (2005), we just omit the proof here. Theorem 3.2 If the two RK methods we applied satisfy the following coefficient conditions   {bkbj−bkakj−bjajk=0,b~mb~n−b~ma~mn−b~na~nm=0,  (3.7) for all $$k,j=1,2,\ldots,r$$ and $$m,n=1,\ldots,s$$, then the numerical algorithm (3.3) is multi-symplectic with the following discrete MSCL: (a) For $$i\notin \{I,I^{'}\}$$,   1Δt∑m=1sb~m(ωi,ml+1−ωi,ml)+1Δx∑k=1rbk(κi+1l,k−κil,k)=0, (3.8) where $${\it{\omega}}_{i,m}^{l+1}=\frac{1}{2}{\,{\rm d}}z_{i,m}^{l+1}\wedge M{\,{\rm d}}z_{i,m}^{l+1}$$, $$\kappa_{i}^{l,k}=\frac{1}{2}{\,{\rm d}}z_{i}^{l,k}\wedge K{\,{\rm d}}z_{i}^{l,k}$$, and so forth. (b) For $$i= I$$,   1Δt∑m=1sb~m(ωI,ml+1−ωI,ml)+1Δx1∑k=1rbk(κI ′,−l,k−κIl,k)=0. (3.9) (c) For $$i= I^{'}$$,   1Δt∑m=1sb~m(ωI ′,ml+1−ωI ′,ml)+1Δx2∑k=1rbk(κI+1l,k−κI ′,+l,k)=0. (3.10) Taking a summation of (3.8)–(3.10) for all the desired indexes $$i$$, we have   [Δx∑i∉{I,I ′}∑m=1sb~m(ωi,ml+1−ωi,ml)+Δx1∑m=1sb~m(ωI,ml+1−ωI,ml)+Δx2∑m=1sb~m(ωI ′,ml+1−ωI ′,ml)]+[Δt∑i∉{I,I ′}∑k=1rbk(κi+1l,k−κil,k)+Δt∑k=1rbk(κI ′,−l,k−κIl,k)+Δt∑k=1rbk(κI+1l,k−κI ′,+l,k)]=0. Using the periodic conditions prescribed previously and taking the summation associatively, the terms in the above second square brackets can be simplified to   Δt∑k=1rbk(κI ′,−l,k−κI ′,+l,k) =Δt2∑k=1rbk(dqI ′,−l,k∧dvI ′,−l,k+dpI ′,−l,k∧dwI ′,−l,k−dqI ′,+l,k∧dvI ′,+l,k−dpI ′,+l,k∧dwI ′,+l,k)=0, where the last equality is similar to (2.6) with the jump condition (3.6) utilized. Thus, we have the following statement. Theorem 3.3 Under the assumptions of Theorem 3.2, the numerical scheme (3.3) preserves the discrete global symplecticity in time   Δx∑i∉{I,I ′}∑m=1sb~mωi,ml+1+Δx1∑m=1sb~mωI,ml+1+Δx2∑m=1sb~mωI ′,ml+1 =Δx∑i∉{I,I ′}∑m=1sb~mωi,ml+Δx1∑m=1sb~mωI,ml+Δx2∑m=1sb~mωI ′,ml, (3.11) with the periodic boundary conditions, i.e., the last two lines of (2.10) imposed. As discussed in Hong et al. (2005) and commented in Hong & Li (2006), two symplectic RK discretizations applied to a normal MSHS in time and space, respectively, leads to a multi-symplectic integrator. Theorem 3.2 tells us that a similar assertion holds for this weak MSHS (2.1). And (2.1) possesses not only this MSCL but also some local and global conservation laws. It is well known in numerical computation that one remarkably important standard in constructing algorithms is the ability of carrying as much of qualitative solution behavior as possible or as many of intrinsic characters as possible. The local and global NCLs, as we state next, are precisely carried under our novelly constructed MSRK discretizations. Theorem 3.4 Under the assumptions of Theorem 3.2, the numerical scheme (3.3) preserves the discrete NCL: (a) For $$i\notin \{I,I^{'}\}$$,   1Δt∑m=1sb~m(Nd(zi,ml+1)−Nd(zi,ml))+1Δx∑k=1rbk(Nf(zi+1l,k)−Nf(zil,k))=0. (3.12) (b) For $$i= I$$,   1Δt∑m=1sb~m(Nd(zI,ml+1)−Nd(zI,ml))+1Δx1∑k=1rbk(Nf(zI ′,−l,k)−Nf(zIl,k))=0. (3.13) (c) For $$i= I^{'}$$,   1Δt∑m=1sb~m(Nd(zI ′,ml+1)−Nd(zI ′,ml))+1Δx2∑k=1rbk(Nf(zI+1l,k)−Nf(zI ′,+l,k))=0. (3.14) Proof. We start with the following calculations   Nd(zi,ml+1)−Nd(zi,ml)=14[(Mzi,ml+1)TMzi,ml+1−(Mzi,ml)TMzi,ml] =14[Δt∑k=1rbk(Mzi,ml)TM∂tZi,ml,k+Δt∑k=1rbk(M∂tZi,ml,k)TMzi,ml  +Δt2∑k=1r∑j=1rbkbj(M∂tZi,ml,k)TM∂tZi,ml,j], (3.15) where the second equality is due to the second equation of (3.3). Making use of the first equation of (3.3), the first term in the right-hand side of (3.15) yields   Δt∑k=1rbk(Mzi,ml)TM∂tZi,ml,k =Δt∑k=1rbk(MZi,ml,k)TM∂tZi,ml,k−Δt2∑k=1r∑j=1rbkakj(M∂tZi,ml,j)TM∂tZi,ml,k. (3.16) Analogously, the second term reads   Δt∑k=1rbk(M∂tZi,ml,k)TMzi,ml =Δt∑k=1rbk(M∂tZi,ml,k)TMZi,ml,k−Δt2∑k=1r∑j=1rbkakj(M∂tZi,ml,k)TM∂tZi,ml,j =Δt∑k=1rbk(MZi,ml,k)TM∂tZi,ml,k−Δt2∑j=1r∑k=1rbjajk(M∂tZi,ml,j)TM∂tZi,ml,k, (3.17) where the first part in the right-hand side is due to the symmetry of the scalar product of vectors and the second the interchange of the summation indexes $$k$$ and $$j$$. By virtue of the first equation given in the MS conditions (3.7), (3.15)–(3.17), this tells us that   1Δt∑m=1sb~m(Nd(zi,ml+1)−Nd(zi,ml))=12∑m=1s∑k=1rb~mbk(MZi,ml,k)TM∂tZi,ml,k. (3.18) As deducted for (3.15)–(3.18), we have   1Δx∑k=1rbk(Nf(zi+1l,k)−Nf(zil,k))=12∑k=1r∑m=1sbkb~m(MZi,ml,k)TK∂xZi,ml,k. (3.19) Combining (3.18) and (3.19) with the last equation of (3.3) yields   1Δt∑m=1sb~m(Nd(zi,ml+1)−Nd(zi,ml))+1Δx∑k=1rbk(Nf(zi+1l,k)−Nf(zil,k))=−12∑k=1r∑m=1sbkb~m(MZi,ml,k)T∇zS(Zi,ml,k)=0, where the second equality is due to the fact that $$(Mz)^T\,\nabla_z S(z)\equiv 0$$. This completes the statement (a). The discussions for the last two statements are similar, and thus are omitted here. □ Theorem 3.5 Under the assumptions of Theorem 3.3, the numerical scheme (3.3) preserves the discrete global NCL   Δx∑i∉{I,I ′}∑m=1sb~mNd(zi,ml+1)+Δx1∑m=1sb~mNd(zI,ml+1)+Δx2∑m=1sb~mNd(zI ′,ml+1) =Δx∑i∉{I,I ′}∑m=1sb~mNd(zi,ml)+Δx1∑m=1sb~mNd(zI,ml)+Δx2∑m=1sb~mNd(zI ′,ml). (3.20) Proof. As discussed for Theorem 3.3, taking a product of $${\it{\Delta}} x {\it{\Delta}} t$$, $${\it{\Delta}} x_1 {\it{\Delta}} t$$ and $${\it{\Delta}} x_2 {\it{\Delta}} t$$ with both sides of (3.12)–(3.14), respectively, and then summing them up, we have   [Δx∑i∉{I,I′}∑m=1sb~m(Nd(zi,ml+1)−Nd(zi,ml))  +Δx1∑m=1sb~m(Nd(zI,ml+1)−Nd(zI,ml))+Δx2∑m=1sb~m(Nd(zI′,ml+1+1)−Nd(zI′,ml))] +[Δt∑i∉{I,I′}∑m=1sbk(Nf(zi+1l,k)−Nf(zil,k))  +Δt∑k=1rbk(Nf(zI′,−l,k)−Nf(zIl,k))+Δt∑k=1rbk(Nf(zI+1l,k)−Nf(zI′,+l,k))=0.] By applying periodic boundary conditions, one reduces the part in the above second square brackets to   Δt∑k=1rbk(Nf(zI′,−l,k)−Nf(zI′,+l,k)) =Δt∑k=1rbk{(qI′,−l,k,wI′,−l,k−pI′,−l,k,vI′,−l,k)−(qI′,+l,k,wI′,+l,k,pI′,+l,k,vI′,+l,k,)}=0, where the last equality can be readily verified by the jump conditions (3.4). Thus it leads to the proof. □ The normalization conservation is an exceptionally important, maybe the most important, conserved quantity in quantum physics (Berezin & Shubin, 1991; Dalfovo et al., 1999; Sulem & Sulem, 1999; Carretero-Gonzáalez et al., 2008). Together with its probability interpretation, it provides the fundamental basis for understanding NLSEs physically and studying NLSEs and other intrinsic conservation laws mathematically. It is well-known that symplectic integrators for Hamiltonian ODEs may not preserve the energy exactly (Hairer et al., 2006). Concerning MS integrators for normal MSHSs, similar phenomena (see, e.g., Reich, 2000; Bridges & Reich, 2006; Hong & Li, 2006; Hong et al., 2007) are frequently observed from numerical experiments for nonlinear RHS $$\nabla_z S(z)$$, but a rigorous proof remains to be given. The existing discussions are commonly restricted to the fact that MSRK methods can preserve the energy of linear normal MSHSs (Reich, 2000; Islas et al., 2001; Hong et al., 2005; Bridges & Reich, 2006; Hong & Li, 2006), i.e., $$S(z)$$ is quadratic in $$z$$. As for a weak MSHS of (2.1), we have the following similar statement. Theorem 3.6 For a system of the abstract form (2.1) (restrict to $$J=1$$ for simplicity), provided that (3.7) is satisfied and $$S(z)$$ is quadratic in $$z$$, i.e., $$S(z)=\frac{1}{2}z^THz+\vec{h}^Tz$$ with $$H$$ a symmetric matrix and $$\vec{h}$$ a vector, then the numerical scheme (3.3) preserves the following discrete local ECL: (a) For $$i\notin \{I,I^{'}\}$$,   1Δt∑m=1sb~m(E(zi,ml+1)−E(zi,ml))+1Δx∑k=1rbk(F(zi+1l,k)−F(zil,k))=0. (3.21) (b) For $$i= I$$,   1Δt∑m=1sb~m(E(zI,ml+1)−E(zI,ml))+1Δx1∑k=1rbk(F(zI ′,−l,k)−F(zIl,k))=0. (3.22) (c) For $$i= I^{'}$$,   1Δt∑m=1sb~m(E(zI ′,ml+1)−E(zI ′,ml))+1Δx2∑k=1rbk(F(zI+1l,k)−F(zI ′,+l,k))=0. (3.23) With periodic boundary conditions prescribed in Theorem 3.3, the algorithm has the following discrete global ECL:   Δx∑i∉{I,I ′}∑m=1sb~mE(zi,ml+1)+Δx1∑m=1sb~mE(zI,ml+1)+Δx2∑m=1sb~mE(zI ′,ml+1) =Δx∑i∉{I,I ′}∑m=1sb~mE(zi,ml)+Δx1∑m=1sb~mE(zI,ml)+Δx2∑m=1sb~mE(zI ′,ml). (3.24) 4. Numerical experiments Based on the theoretical results in the previous section, numerical experiments presented in this section reveal numerical characters of the MSRK methods applied to the system (2.1). 4.1 Numerical experiment A: a single-delta potential Here we concern the single-delta potential case that $$V_{\rm ext}(x)=\gamma \delta(x)$$ in (1.8), which can be reformulated as the weak MSHS (2.1) as investigated previously. We discretize the system (2.1) by the MSRK method with $$r=1$$ and $$s=1$$, i.e., $$a_{11}=\tilde{a}_{11}=1/2$$, $$c_1=\tilde{c}_1=1/2$$ and $$b_1=\tilde{b}_1=1$$, which is an implicit second-order MSRK method (in short, denoted by MM2) and is equivalent to the central box scheme (see, e.g., Reich, 2000; Bridges & Reich, 2006; Hong & Li, 2006). For $$i\notin \{I,I^{'}\}$$, it reads   M(zi+12l+1−zi+12l)Δt+K(zi+1l+12−zil+12)Δx=∇zS(zi+12l+12) (4.1) with that   zi+12l=(zi+1l+zil)/2,  zi+12l+1=(zi+1l+1+zil+1)/2,  zil+12=(zil+1+zil)/2,zi+1l+12=(zi+1l+1+zi+1l)/2,  zi+12l+12=(zi+1l+1+zi+1l+zil+1+zil)/4; for $$i=I$$, according to ($$3.3$$B), one can use $$z^{l+\frac{1}{2}}_{I^{'}\!\!,-}$$ and $${\it{\Delta}} x_1$$ to substitute $$z^{l+\frac{1}{2}}_{i+1}$$ and $${\it{\Delta}} x$$ appearing in (4.1), respectively; for $$i=I^{'}$$, according to ($$3.3$$C), one can use $$z^{l+\frac{1}{2}}_{I^{'}\!\!,+}$$, $$z^{l+\frac{1}{2}}_{I+1}$$ and $${\it{\Delta}} x_2$$ to substitute $$z^{l+\frac{1}{2}}_{i}$$, $$z^{l+\frac{1}{2}}_{i+1}$$ and $${\it{\Delta}} x$$, respectively. Here   zI+12l=(zI ′,−l+zIl)/2,zI+12l+12=(zI ′,−l+1+zI ′,−l+zIl+1+zIl)/4,zI ′+12l=(zI ′,+l+zI+1l)/2,zI ′+12l+12=(zI ′,+l+1+zI ′,+l+zI+1l+1+zI+1l)/4 and so on. Relevant to (3.4), the jump conditions for MM2 are specified as   zI ′,+l−zI ′,−l=AjumpzI ′,−landzI ′,+l+1−zI ′,−l+1=AjumpzI ′,−l+1. (4.2) For the purpose of numerical comparisons, we will consider two more implicit second-order nonmulti-symplectic methods here: one is about the symplectic RK method with $$r=1$$ in time and an immersed interface method (see, e.g., Bai & Wang, 2015; Bai, 2016 and references therein) in space (in short, denoted by M2IIM), the other is about the trapezoidal scheme in time and IIM in space (in short, denoted by T2IIM). For the spatial index $$i\notin \{I,I+1\}$$, the scheme M2IIM yields   {2(qil+1−qil)Δt+pi+1l+12−2pil+12+pi−1l+12Δx2=2g[(qil+12)2+(pil+12)2]pil+12,2(pil+1−pil)Δt−qi+1l+12−2qil+12+qi−1l+12Δx2=−2g[(qil+12)2+(pil+12)2]qil+12  (4.3) with $$p^{l+\frac{1}{2}}_{i}=\big(p^{l+1}_{i}+p^{l}_{i}\big)/2$$ processed previously; and for $$i=I$$ and $$i=I+1$$, the scheme reads, respectively,   {2(qIl+1−qIl)Δt+ϕ1(pIl+12)=2g[(qIl+12)2+(pIl+12)2]pIl+12,2(pIl+1−pIl)Δt−ϕ1(qIl+12)=−2g[(qIl+12)2+(pIl+12)2]qIl+12  (4.3B) and   {2(qI+1l+1−qI+1l)Δt+ϕ2(pI+1l+12)=2g[(qI+1l+12)2+(pI+1l+12)2]pI+1l+12,2(pI+1l+1−pI+1l)Δt−ϕ2(qI+1l+12)=−2g[(qI+1l+12)2+(pI+1l+12)2]qI+1l+12  (4.3C) with   {ϕ1(yI)≡βI,1yI−1+βI,2yI+βI,3yI+1ϕ2(yI+1)≡βI+1,1yI+βI+1,2yI+1+βI+1,3yI+2,  (4.4) where the linear coefficients $$\beta_{I,1}$$, $$\beta_{I,2}$$ and $$\beta_{I,3}$$ solve the following linear equations   {βI,1+βI,2+(1+2γxI+1)βI,3=0,xI−1βI,1+xIβI,2+xI+1βI,3=0,xI−12βI,1+xI2βI,2+xI+12βI,3=2,  (4.5) and $$\beta_{I+1,1}$$, $$\beta_{I+1,2}$$ and $$\beta_{I+1,3}$$ solve   { (1−2γxI)βI+1,1+βI+1,2+βI+1,3=0,xIβI+1,1+xI+1βI+1,2+xI+2βI+1,3=0,xI2βI+1,1+xI+12βI+1,2+xI+22βI+1,3=2.  (4.6) Let $$\tilde{z}=(q,p)^T$$ and $$\varphi(\tilde{z})=\big(2g(q^2+p^2)p,-2g(q^2+p^2)q\big)^T$$ here; then the scheme T2IIM is obtained by means of making the substitutions of   φ(z~il+1)+φ(z~il)2,φ(z~Il+1)+φ(z~Il)2andφ(z~I+1l+1)+φ(z~I+1l)2 for the RHSs of (4.3), ($$4.3$$B) and ($$4.3$$C), respectively. It is well known that conservation laws are not only important in physics are but also powerful tools to test and judge numerical algorithms for differential equations, especially for those which could not be solved analytically. To investigate numerically the intrinsic conservation laws addressed before, we need to define, in terms of different numerical schemes, the corresponding discrete analogies of them. Concerning the MS scheme MM2, the discrete local ECL $$(E_{{\rm loc}}^*)_i^l$$, for $$i\notin \{I,I^{'}\}$$ is given by (see Islas et al., 2001; Bridges & Reich, 2006; Hong & Li, 2006)   (Eloc∗)il=1Δt[(S(zi+12l+1)−12(zi+12l+1)TKzi+1l+1−zil+1Δx)−(S(zi+12l)−12(zi+12l)TKzi+1l−zilΔx)]  +1Δx[(zi+1l+12)TKzi+1l+1−zi+1lΔt−(zil+12)TKzil+1−zilΔt]; for $$i=I$$, the expression can be obtained by replacing $$z_{i+1}$$ and $${\it{\Delta}} x$$ with $$z_{I^{'}\!\!,-}$$ and $${\it{\Delta}} x_1$$, respectively; for $$i=I^{'}$$, by replacing $$z_{i}$$, $$z_{i+1}$$ and $${\it{\Delta}} x$$ with $$z_{I^{'}\!\!,+}$$, $$z_{I+1}$$ and $${\it{\Delta}} x_2$$, respectively. Actually, $$(E_{{\rm loc}}^*)_i^l$$ is the numerical error for the ECL (1.12) at the point $$(x_{i+\frac{1}{2}},t_{l+\frac{1}{2}})$$. For simplicity, we hence define the following maximum errors for ECL   (Emax∗)l=maxi{(Eloc∗)il}, (4.7) where $$i$$ runs over all spatial index values. Similarly the discrete local NCL $$(N_{{\rm loc}}^*)_i^l$$ corresponding to (1.14), for $$i\notin \{I,I^{'}\}$$, is given by   (Nloc∗)il=1Δt[Nd(zi+12l+1)−Nd(zi+12l)]+1Δx[Nf(zi+1l+12)−Nf(zil+12)]. And the expressions of $$(N_{{\rm loc}}^*)_{I}^{l}$$ and $$(N_{{\rm loc}}^*)_{I^{'}}^{l}$$ can be similarly obtained as for the local ECL. In response to (4.7), we define the following maximum errors for local NCL   (Nmax∗)l=maxi{(Nloc∗)il}. (4.8) Moreover, the discrete total energy and normalization for MM2 are given by   (EMS)l=Δx∑i∉{I,I′}(S(zi+12l)−12(zi+12l)TKzi+1l−zilΔx)+Δx1(S(zI+12l)−12(zI+12l)TKzI′,−l−zIlΔx1)+Δx2(S(zI′+12l)−12(zI′+12l)TKzI+1l−zI′,+lΔx2) and   (NMS)l=Δx∑i∉{I,I ′}Nd(zi+12l)+Δx1Nd(zI+12l)+Δx2Nd(zI ′+12l), respectively. Thus the relative errors, also referred to as the propagating errors, for the discrete global ECL and NCL are, respectively, defined by   (EMSrel)l=(EMS)l−(EMS)0|(EMS)0|and(NMSrel)l=(NMS)l−(NMS)0|(NMS)0|. (4.9) We remark here that the two non-MS schemes are both discretized for the second-order spatial derivative directly and have no corresponding approximations for the first-order ones, and thus they are not quite suitable to be used in investigating the local conservation laws at the present stage. The corresponding discrete total ones, however, are given by   (EnMS)l=Δx2∑ig[(qil)2+(pil)2]2−Δx2∑i∉{I,I+1}[qi−1l−2qil+qi+1lΔx2qil+pi−1l−2pil+pi+1lΔx2pil]−Δx2[ϕ1(qIl)qIl+ϕ2(qI+1l)qI+1l+ϕ1(pIl)pIl+ϕ2(pI+1l)pI+1l] and   (NnMS)l=Δx∑i[(qil)2+(pil)2] with $$\phi_1$$ and $$\phi_2$$ specified in (4.4). Similarly, the corresponding propagating errors $$(\mathscr{E}_{{\rm nMS}}^{{\rm rel}})^{l}$$ and $$(\mathscr{N}_{{\rm nMS}}^{{\rm rel}})^{l}$$ can be defined as addressed in (4.9). The numerical experiment we will investigate is the case of an attractive nonlinearity, i.e., $$g = -1$$ for the nonlinear Schrödinger equation (1.19). For the sake of better testing our methods, the experiment we choose to implement here possesses the analytic solution (Witthaut et al., 2005)   Ψ(x,t)=ψ(x)exp⁡(i(2γ−1)2t/8), (4.10) where   ψ(x)={λsech(λ(x−x0)),x>0λsech(λ(x+x0)),x<0  with $$\lambda$$ and $$x_0$$ is being subjected to   λ=12−γandtanh(λx0)=γλ. (4.11) For the corresponding stationary nonlinear Schrödinger equations, bound state solutions exist provided $$\gamma<1/4$$ (Witthaut et al., 2005). The above analytic solution concerns the evolution of this bound state solution with the energy $$\mathscr{E}=-(2\gamma-1)^2/8$$. Due to the bound-state-existence condition, in our numerical experimentation here, we adopt a delta-potential barrier with $$\gamma=0.2$$, $$\lambda$$ and $$x_0$$ are then determined by (4.11). Since the exact solution $${\it{\Psi}}(t,x)$$ is exponentially small away from $$x = 0$$ for any fixing $$t$$, the periodic boundary conditions are applied to both the MS and non-MS schemes in numerics. In practical calculations, we take $$x_L = -30$$ and $$x_R = 30$$. We consider two circumstances, which admit, after the first equidistant partition, whether the interface point $$x=0$$ is a uniform grid point or not. If it is, we choose $$N=600$$ and $$N=599$$ for the other. The former means that the spatial stepsize $${\it{\Delta}} x=(x_R-x_L)/N=0.1$$, the latter that $${\it{\Delta}} x=(x_R-x_L)/N\approx 0.1002$$. The temporal stepsize is fixed to be $${\it{\Delta}} t=0.1$$ and the temporal interval $$[0,200]$$. The fixed point iteration method is utilized to solve the nonlinear algebraic systems generated by the numerical schemes, and each iteration will be terminated when the maximum absolute error of the adjacent two iterative values is less than $$10^{-15}$$. The two circumstances are implemented independently. The numerical results plotted here put two legends on them, ‘S1’ for the first circumstance $$N=600$$ and ‘S2’ for the second one $$N=599$$. Figure 1 exhibits the maximum global errors of the numerical wave function obtained by the three schemes. Here we use   max0≤i≤N(max(|qil−q(xi,tl)|,|pil−p(xi,tl)|)) to denote the maximum global error of the wave function at time $$t_l$$, where $$p_i^l$$ is the numerical solution of $$q(x,t)$$ at $$(x_i,t_l)$$ and $$q(x_l,t_i)$$ is the corresponding analytic solution due to (4.10). It follows from the three plots that, all the errors obtained are all in reasonable oscillation and normal linear growth. Considering the fact that the schemes we implemented here are all of second-order accuracy, the errors seem comparatively nice for the temporal step size $${\it{\Delta}} t=0.1$$ we used. The errors for the MS scheme, whether $$x=0$$ is a uniform grid point or not, are both of $$\mathscr{O}(10^{-3})$$ and have almost the same evolutions with time. But those, which are obtained by the other two non-MS schemes, behave very differently for the two circumstances: when $$x=0$$ is a uniform grid point, the errors are of $$\mathscr{O}(10^{-4})$$; and they increase by an order of magnitude, i.e., are of $$\mathscr{O}(10^{-3})$$ for the other. Fig. 1. View largeDownload slide The maximum errors of the wave function obtained by MM2, M2IIM and T2IIM. Fig. 1. View largeDownload slide The maximum errors of the wave function obtained by MM2, M2IIM and T2IIM. Figure 2 displays the propagating errors of the total normalization and energy obtained by MM2. They are given in (4.9). The errors of the normalization, whether $$x=0$$ is a grid point or not, are both in the scale of $$\mathscr{O}(10^{-15})$$, almost the roundoff error of computer. Simultaneously they have similar reasonable linear growths. The errors of the total energy are both in the scale of $$\mathscr{O}(10^{-4})$$ and have stable oscillations something like sinusoidal waves. Fig. 2. View largeDownload slide The propagating errors of the total normalization and energy obtained by MM2, respectively. Fig. 2. View largeDownload slide The propagating errors of the total normalization and energy obtained by MM2, respectively. The propagating errors of the total normalization, which are obtained by M2IIM and T2IIM, are exhibited in Fig. 3. We can see for the two non-MS schemes that the errors for the same conserved quantity look very different whether the interface point is a uniform grid point or not. When it is, the error of the normalization obtained by M2IIM is of $$\mathscr{O}(10^{-14})$$, which is almost as good as that by the MS scheme MM2. In the other circumstance, however, it increases to $$\mathscr{O}(10^{-8})$$, about $$6$$ orders of magnitude lost. For T2IIM, the errors are only of $$\mathscr{O}(10^{-9})$$ and $$\mathscr{O}(10^{-8})$$ for the two circumstances, respectively. Numerical results reveal that the precise normalization conservation is one remarkable advantage of the MS scheme compared with the non-MS schemes. Fig. 3. View largeDownload slide The propagating errors of the total normalization, the left two plots are obtained by M2IIM and the right by T2IIM. Fig. 3. View largeDownload slide The propagating errors of the total normalization, the left two plots are obtained by M2IIM and the right by T2IIM. Figure 4 shows the errors of the total energy obtained by M2IIM and T2IIM. For the first circumstance, the errors obtained by the two non-MS schemes are both of orders $$\mathscr{O}(10^{-9})$$; for the second one, however, they both grow rapidly from $$\mathscr{O}(10^{-9})$$ to $$\mathscr{O}(10^{-5})$$ and then stay in the latter scale. The corresponding errors obtained by MM2, which are displayed in the right plot of Fig. 2, are only of orders $$\mathscr{O}(10^{-4})$$ for both the two circumstances and several orders worse than that given here. We think the reason might be due to different approximations of $$\partial_x v$$ and $$\partial_x w$$ arising in the energy density. Regarding MM2, $$\partial_x v$$ and $$\partial_x w$$ are regarded as the first-order spatial derivatives and thus approximated by the numerical values of $$v$$ and $$w$$, which admit discontinuities at interface points. For the other two non-MS schemes, however, $$\partial_x v$$ and $$\partial_x w$$ are regarded as the second-order ones and approximated by the numerical values of $$p$$ and $$q$$. We remind here that $$\partial_x v$$ and $$\partial_x w$$ are continuous. The discontinuities of $$v$$ and $$w$$ we highlighted in our proposed MS disretizations may be the drawbacks in such approximations. Taking this into account, Bai (2016) proposed MS Runge–Kutta–Nyström methods applied to the delta-potential NLSEs and numerically reveal that the total energy obtained by the implemented MS scheme is preserved much better than by the non-MS one. Fig. 4. View largeDownload slide The propagating errors of the total energy, the left two plots are obtained by M2IIM and the right by T2IIM. Fig. 4. View largeDownload slide The propagating errors of the total energy, the left two plots are obtained by M2IIM and the right by T2IIM. The parity symmetry is one of the three important discrete symmetries, i.e., charge conjugate, parity and time reversal (CPT) symmetries in quantum physics (see, e.g., Cartarius et al., 2012). The exact solution (4.10) we investigate here is of even parity, i.e., the even symmetry of the wave function $${\it{\Psi}}(x,t)={\it{\Psi}}(-x,t)$$. We use   max0≤i≤N(max(|pil−pN−il|,|qil−qN−il|)) to denote the maximum error of the parity conservation at time $$t_l$$. The errors for the two circumstances are presented in Fig. 5. For the MS scheme MM2, the errors are both of orders $$\mathscr{O}(10^{-12})$$; for either M2IIM or T2IIM, the errors are in the scale of $$\mathscr{O}(10^{-10})$$ and $$\mathscr{O}(10^{-11})$$ for the two circumstances, respectively. If the six errors are amplified, we find that they have very similar evolution behaviors. Fig. 5. View largeDownload slide The maximum errors in the parity conservation obtained by MM2, M2IIM and T2IIM. Fig. 5. View largeDownload slide The maximum errors in the parity conservation obtained by MM2, M2IIM and T2IIM. Figure 6 shows the errors of the local conservation laws obtained by the MS scheme MM2. The maximum errors of the local NCL are both of $$\mathscr{O}(10^{-16})$$, almost the roundoff error of the computer, for the two circumstances. The numerical results match our theoretical analysis in Theorem 2.5 very well. The maximum errors of the local ECL, given by (4.7), are in the scale of $$\mathscr{O}(10^{-6})$$ and $$\mathscr{O}(10^{-3})$$ for the two circumstances, respectively. Fig. 6. View largeDownload slide The maximum errors of the local conservation laws by MM2, the left two plots for the local NCL and the right for the local ECL. Fig. 6. View largeDownload slide The maximum errors of the local conservation laws by MM2, the left two plots for the local NCL and the right for the local ECL. 4.2 Numerical experiment B: a double-delta potential For comparison, the numerical experiment we will consider here is the case of the repulsive nonlinearity and the symmetric double-delta potential well, i.e., $$g = 1$$ and $$V_{\rm ext}(x)=\tilde{\gamma}\big[\delta(x-\alpha)+\delta(x+\alpha)\big]$$ ($$\tilde{\gamma}<0,\, \alpha>0$$), respectively, for the NLSE (1.19). The bound state problems for the corresponding stationary linear case have been investigated using the explicit jump IIM and the Peskin’s IB Bai & Wang (2015). The nonlinear case, however, is far distinguished. By virtue of the strategies proposed in Kivshar & Luther-Davies (1998) and Witthaut et al. (2005, 2009), we can construct a solitary wave solution   Ψ~(x,t)=ψ~(x)exp⁡(iλ~2t), (4.12) where   ψ~(x)={2λ~csch(2λ~(x−x~0)),x>α,λ~tan(λ~x),−α≤x≤α,2λ~csch(2λ~(x+x~0)),x<−α.  We stress here that the above is a bound state solution to this nonlinear Schrödinger equation and evolves with the energy $$\tilde{\mathscr{E}}=-\tilde{\lambda}^2$$. One could verify that the above solution is of odd parity, i.e., $$\tilde{{\it{\Psi}}}(-x,t)=-\tilde{{\it{\Psi}}}(x,t)$$, which is in contrast to the even parity proposed in Experiment A. When $$\alpha$$ is fixed, the parameters $$\tilde{\gamma}$$, $$\tilde{\lambda}$$ and $$\tilde{x}_0$$ can be determined by   {2csch(2λ~(α−x~0))=tan⁡(λ~α),2λ~csch(2λ~(α−x~0))coth(2λ~(α−x~0))+λ~sec2⁡(λ~α)=−2γ~tan⁡(λ~α),2λ~tan⁡(λ~α)−2λ~2α+22λ~(coth(2λ~(α−x~0))−1)=1.  (4.13) Here the first equation of (4.13) comes from the continuity of the wave function at the interface points, the second from the jumps and the third from the normalization requirement that   ∫−∞+∞|Ψ~(x,t)|2dx=1. The existence of bound state solutions needs some more requirements, e.g., $$\tilde{x}_0<\alpha$$, $$\tilde{\lambda}\,\alpha<{\pi}/2$$ and so on. In numerics, we adopt $$\alpha=3$$ as investigated for the stationary linear case in Bai & Wang (2015). And the other three parameters we choose are $$\tilde{\lambda}=0.345692582829775$$, $$\tilde{x}_0=1.444829457402318$$ and $$\tilde{\gamma}=-0.775835479368852$$, which solve the transcendental equations (4.13) numerically with the terminating criterion $$10^{-15}$$. Since the exact solution $$\tilde{{\it{\Psi}}}(t,x)$$ is also exponentially small away from $$x = 0$$ for any fixing $$t$$, the periodic boundary conditions are still used in numerical implementation. The spatial domain, the temporal interval, the stop criterion for the nonlinear iteration and the temporal stepsize are all taken to the same as utilized in Experiment A. With regard to the choices of spatial stepsize, we take two circumstances into account. The first concerns that both the two interface points $$x=\pm \alpha$$ are uniform partition grid points and the second is that both are not. For comparisons, we adopt the same two spatial stepsizes as used in Experiment A, which also correspond to the two circumstances in this experiment. For simplicity, we only consider two schemes, the MS one MM2 and the non-MS one T2IIM proposed before. And all the expressions could be obtained similarly and thus would not be presented here. The two circumstances of different spatial stepsizes are also implemented independently and put the two legends ‘D1’ and ‘D2’ on the plots, respectively. Figure 7 presents the maximum errors of numerical solutions obtained by using the two schemes. The two errors are both in the scale of $$\mathscr{O}(10^{-3})$$, just the same magnitude to that for the single-delta case as exhibited in Fig. 1. But the ones by T2IIM given here behave very different from that for the single-delta case: (1) the two errors are of orders $$\mathscr{O}(10^{-1})$$ and $$\mathscr{O}(10^{-2})$$, respectively, which are several orders of magnitude higher than that for the single-delta case as shown in Fig. 1; (2) the first error is one order of magnitude higher than the second one in this double-delta case, but it just reverses in the preceding single-delta case. Fig. 7. View largeDownload slide The maximum errors of numerical solutions, the left is obtained by MM2 and the right by T2IIM. Fig. 7. View largeDownload slide The maximum errors of numerical solutions, the left is obtained by MM2 and the right by T2IIM. Figure 8 exhibits the errors of the total normalization. Both the errors obtained by MM2 here reach almost the roundoff error of computer, just like that in the single-delta case given in the left plot of Fig. 2. The ones by T2IIM are both of $$\mathscr{O}(10^{-7})$$, about two and one orders of magnitude higher than that in the single-delta case, respectively. Fig. 8. View largeDownload slide The propagating errors of the total normalization obtained by MM2 and T2IIM. Fig. 8. View largeDownload slide The propagating errors of the total normalization obtained by MM2 and T2IIM. The propagating errors of the total energy are given in Fig. 9. Those two obtained by MM2 are in the scale of $$\mathscr{O}(10^{-3})$$, about one order of magnitude higher than that in the single-delta case as displayed in Fig. 2. The two errors by T2IIM are of orders $$\mathscr{O}(10^{-7})$$ and $$\mathscr{O}(10^{-5}),$$ respectively, while the corresponding ones are $$\mathscr{O}(10^{-9})$$ and $$\mathscr{O}(10^{-5})$$ in the single-delta case as shown in Fig. 4. Fig. 9. View largeDownload slide The propagating errors of the total energy, the left one obtained by MM2 and the other by T2IIM. Fig. 9. View largeDownload slide The propagating errors of the total energy, the left one obtained by MM2 and the other by T2IIM. Figure 10 gives the errors in the parity conservation. The errors are obtained by using MM2 and are both in the scale of $$\mathscr{O}(10^{-7})$$, almost five orders of magnitude higher than that in the single-delta case as exhibited in Fig. 5. The two errors by T2IIM are of orders of $$\mathscr{O}(10^{-8})$$ and $$\mathscr{O}(10^{-10})$$, respectively, about two and one orders of magnitude higher than that in the single-delta case. Fig. 10. View largeDownload slide The maximum errors in the parity conservation, the left and the middle are obtained by MM2, the right by T2IIM. Fig. 10. View largeDownload slide The maximum errors in the parity conservation, the left and the middle are obtained by MM2, the right by T2IIM. Figure 11 shows the maximum errors of the two local conservation laws. The errors of the local NCL arrive at the roundoff error, just like that in the single-delta case exhibited in Fig. 6. The two errors of the local ECL are both in the scale of $$\mathscr{O}(10^{-3})$$, while the corresponding ones are $$\mathscr{O}(10^{-6})$$ and $$\mathscr{O}(10^{-3})$$, respectively, in the single-delta case. Fig. 11. View largeDownload slide The maximum errors of local conservation laws obtained by MM2, the left and the middle for the local NCL, and the right for the local ECL. Fig. 11. View largeDownload slide The maximum errors of local conservation laws obtained by MM2, the left and the middle for the local NCL, and the right for the local ECL. 5. Conclusions We make an attempt to propose a weak MSHS reformulation for the NLSE with delta potentials. Based on this weak reformulation, some intrinsic local and global conservation laws are newly addressed, a kind of novel concatenating RK discretizations are constructed and the multi-symplecticity of the concatenating algorithms are revealed. Under the MSRK discretizations, some discrete properties are investigated theoretically and, among them, what is important is the precise preservation of the local and global NCLs. Extensive numerical experimentation validates our method. Comparisons with non-MS schemes reveal some priorities of our MS formulations, and the first is of course the roundoff-error preservation of the NCLs. Furthermore, the longitudinal comparisons made between two circumstances show that for either the single- or double-delta case the errors by our MS schemes vary much less than that by the non-MS ones whether the interface points are uniform grid points or not; the transverse comparisons made between two cases reveal that, the errors by our MS schemes change much less than that by non-MS ones when the number of delta potentials change from one to two. It might be concluded that the MS schemes have more stable dynamical numerical behaviors than the non-MS ones. In physics, there exist a large number of important PDEs which admit various jumps or discontinuities, e.g., the nonlinear Dirac equations, the Klein–Gordon equations and many other coupled nonlinear equations with several delta potentials, some hydrodynamic equations with discontinuous coefficients or singular sources (Landau & Lifshitz, 1987; Haus & Wong, 1996; Kivshar & Luther-Davies, 1998; Bridges & Reich, 2006; Cartarius et al., 2012 and references therein) and so on. They may be explored in a similar way as proposed in this paper. this context, although further research is needed. Acknowledgements We would like to thank the anonymous referees for their valuable comments and suggestions. Funding National Natural Science Foundation of China (Nos.11001125 and 91130003) and A Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions. References Arnold, V. I. ( 1999) Mathematical Methods of Classical Mechanics , 2nd edn. New York: Springer. Bai, J. ( 2016) Multi-symplectic Runge-Kutta-Nyström methods for nonsmooth nonlinear Schrödinger equations. J. Math. Anal. Appl. , 444, 721– 736. Google Scholar CrossRef Search ADS   Bai, J. & Wang, L. ( 2015) EJIIM for the stationary Schrödinger equations with delta potential wells. Appl. Math. Comput. , 254, 113– 124. Berezin, F. A. & Shubin, M. A. ( 1991) The Schrödinger Equation . Dordrecht: Kluwer Academic Publishers. Google Scholar CrossRef Search ADS   Bridges, T. J. ( 1997) Multi-symplectic structures and wave propagation. Math. Proc. Camb. Phil. Soc. , 121, 147– 190. Google Scholar CrossRef Search ADS   Bridges, T. J. & Reich, S. ( 2006) Numerical methods for Hamiltonian PDEs. J. Phys. A Math. Gen. , 39, 5287– 5320. Google Scholar CrossRef Search ADS   Carretero-González, R., Frantzeskakis, D. J. & Kevrekidis, P. G. ( 2008) Nonlinear waves in Bose-Einstein condensates: physical relevance and mathematical techniques. Nonlinearity , 21, R139– R202. Google Scholar CrossRef Search ADS   Cartarius, H., Haag, D., Dast, D. & Wunner, G. ( 2012) Nonlinear Schrödinger equation for a $$\mathscr{PT}$$-symmetric delta-function double well. J. Phys. A Math. Theor. , 45, 444008 ( 15 pp). Google Scholar CrossRef Search ADS   Dalfovo, F., Giorgini, S., Pitaevskii, L. P. & Stringari, S. ( 1999) Theory of Bose-Einstein condensation in trapped gases. Rev. Mod. Phys. , 71, 463– 512. Google Scholar CrossRef Search ADS   Hairer, E., Lubich, C. & Wanner, G. ( 2006) Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations , 2nd edn. Berlin: Springer. Haus, H. A. & Wong, W. S. ( 1996) Solitons in optical communications. Rev. Mod. Phys. , 68, 423– 444. Google Scholar CrossRef Search ADS   Hong, J. & Li, C. ( 2006) Multi-symplectic Runge-Kutta methods for nonlinear Dirac equations. J. Comput. Phys. , 211, 448– 472. Google Scholar CrossRef Search ADS   Hong, J., Liu, H. & Sun, G. ( 2005) The multi-symplecticity of partitioned Runge-Kutta methods for Hamiltonian PDEs. Math. Comput. , 75, 167– 181. Google Scholar CrossRef Search ADS   Hong, J., Liu, X. & Li, C. ( 2007) Multi-symplectic Runge-Kutta-Nyström methods for nonlinear Schrödinger equations with variable coefficients. J. Comput. Phys. , 226, 1968– 1984. Google Scholar CrossRef Search ADS   Islas, A. L., Karpeev, D. A. & Schober, C. M. ( 2001) Geometric integrators for the nonlinear Schrödinger equation. J. Comput. Phys. , 173, 116– 148. Google Scholar CrossRef Search ADS   Kivshar, Y. S. & Luther-Davies, B. ( 1998) Dark optical solitons: physics and applications. Phys. Rep. , 298, 81– 197. Google Scholar CrossRef Search ADS   Kong, L., Hong, J., Ji, L. & Zhu, P. ( 2015) Compact and efficient conservative schemes for coupled nonlinear Schrödinger equations. Numer. Meth. PDE , 31, 1814– 1843. Google Scholar CrossRef Search ADS   Landau, L. D. & Lifshitz, E. M. ( 2008) Fluid Mechanics , 2nd edn. Beijing: Beijing World Publishing Corporation. McLachlan, R. I., Ryland, B. N. & Sun, Y. ( 2014) High order multisymplectic Runge-Kutta methods. SIAM J. Sci. Comput. , 36, A2199– A2226. Google Scholar CrossRef Search ADS   Olver, P. J. ( 1993) Applications of Lie Groups to Differential Equations , 2nd edn. New-York: Springer. Google Scholar CrossRef Search ADS   Pethick, C. J. & Smith, H. ( 2002) Bose-Einstein Condensation in Dilute Gases . New-York: Cambridge University Press. Rapedius, K. & Korsch, H. J. ( 2009) Resonance solutions of the nonlinear Schrödinger equation in an open double-well potential. J. Phys. B , 42, 044005 ( 12 pp). Google Scholar CrossRef Search ADS   Reich, S. ( 2000) Multi-symplectic Runge-Kutta collocation methods for Hamiltonian wave equation. J. Comput. Phys. , 157, 473– 499. Google Scholar CrossRef Search ADS   Ryland, B. N. & McLachlan, R. I. ( 2008) On multi-symplecticity of partitioned Runge-Kutta methods. SIAM J. Sci. Comput. , 30, 1318– 1340. Google Scholar CrossRef Search ADS   Sulem, C. & Sulem, P.-L. ( 1999) The Nonlinear Schrödinger Equation: Self-Focusing and Wave Collapse . New York: Springer. Witthaut, D., Mossmann, S. & Korsch, H. J. ( 2005) Bound and resonance states of the nonlinear Schrödinger equation in simple model systems. J. Phys. A Math. Gen. , 38, 1777– 1792. Google Scholar CrossRef Search ADS   Witthaut, D., Rapedius, K. & Korsch, H. J. ( 2009) The nonlinear Schrodinger equation for the delta-comb potential: quasi-classical chaos and bifurcations of periodic stationary solutions. J. Nonlinear Math. Phys. , 16, 207– 233. Google Scholar CrossRef Search ADS   © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

Weak multi-symplectic reformulation and geometric numerical integration for the nonlinear Schrödinger equations with delta potentials

Download PDF
31 pages

Loading next page...
 
/lp/ou_press/weak-multi-symplectic-reformulation-and-geometric-numerical-9MCrrW09VW

References (26)

Publisher
Oxford University Press
Copyright
© The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.
ISSN
0272-4979
eISSN
1464-3642
DOI
10.1093/imanum/drw062
Publisher site
See Article on Publisher Site

Abstract

Abstract We discuss multi-symplectic (MS) geometric numerical integration for a class of important models in condensed matter physics, i.e., nonlinear Schrödinger equations (NLSEs) with inclusions of finite delta potentials. Due to the inclusions, they cannot be reformulated as multi-symplectic Hamiltonian systems (MSHSs) such as normal NLSEs. With regard to integration, we propose a kind of weak MS reformulations and address some intrinsic local and global conservation laws. Based on theoretical suggestions, concatenating multi-symplectic Runge–Kutta (MSRK) methods are novelly constructed for these weak MSHSs, and some local and global conservation laws under MSRK discretizations are investigated analytically. Extensive numerical experimentation is then presented including cases of single- and double-delta potentials and cases of uniform and nonuniform gird points. Various numerical comparisons made between MS and non-MS schemes show us that the remarkable advantage of the former proposed here is the precise preservation of local and global normalization conservation laws. Other conservative properties, by using the two classes of schemes, can be preserved within some accuracies over long-time evolution. But, in general, the MS formulations behave more stably than the other in the numerics. 1. Introduction Nonlinear Schrödinger equation (NLSE) is one of the most important nonlinear models in modern science. It appears in quantum field theory, condensed matter, nonlinear optics, soliton waves, the theory of partial differential equations (PDEs), the theory of dynamical systems and so forth, in physics and mathematics (see, e.g., Berezin & Shubin, 1991; Haus & Wong, 1996; Dalfovo et al., 1999; Sulem & Sulem, 1999; Carretero-González et al., 2008). We focus on the model based on Bose–Einstein condensation, which was proposed theoretically in 1925 and realized experimentally in 1995 for rubidium gas. And it provides opportunities for exploring quantum phenomena on a macroscopic scale (Dalfovo et al., 1999; Pethick & Smith, 2002; Carretero-González et al., 2008). Let us briefly review the NLSE or Gross–Pitaevskii (GP) equation in condensed matter physics (Dalfovo et al., 1999; Pethick & Smith, 2002). We consider a sufficiently dilute ultracold atomic gas which consists of $$N$$ interacting bosons of mass $$m$$ and is confined by a time-independent external potential $$V_{{\rm ext}}({\bf r})$$. Then the many-body Hamiltonian, in terms of the second quantization, is given by   H^ =∫d3rΨ^†(r,t)[−ℏ22m∇2+Vext(r)]Ψ^(r,t) +12∫d3rd3r′Ψ^†(r,t)Ψ^†(r′,t)V(r−r′)Ψ^(r′,t)Ψ^(r,t), (1.1) where $$\hat{{\it{\Psi}}}({\bf r},t)$$ and $$\hat{{\it{\Psi}}}^\dagger({\bf r},t)$$ are the boson annihilation and creation field operators, respectively, $$V({\bf r}-{\bf r}')$$ is the two-body interatomic potential. The Heisenberg evolution equation for the field operator, i.e., $${\bf i}\hbar\, \partial_t\hat{{\it{\Psi}}} =[\hat{{\it{\Psi}}},\hat{H}]$$, gives rise to   iℏ∂tΨ^(r,t)=[−ℏ22m∇2+Vext(r)+∫d3r′Ψ^†(r′,t)V(r−r′)Ψ^(r′,t)]Ψ^(r,t). (1.2) Proceeding further, one writes   Ψ^(r,t)=Ψ(r,t)+Ψ^′(r,t), (1.3) where $${\it{\Psi}}({\bf r},t) \equiv \langle \hat{{\it{\Psi}}}({\bf r},t)\rangle$$ is the expectation value of the field operator and is commonly called the macroscopic wavefunction of the condensate. Below the critical temperature, the noncondensate part $$\hat{{\it{\Psi}}}'({\bf r},t)$$ is negligible and then the mean-field approximation, that is, $$\hat{{\it{\Psi}}}({\bf r},t) = {\it{\Psi}}({\bf r},t)$$ can be used. For a dilute ultracold gas at low energy, only binary collisions are relevant and can be characterized by the $$s$$-wave scattering length $$a$$, which is independent of the details of the two-body potential. Thus, the interatomic potential can be replaced by an effective delta-function interaction potential,   V(r−r′)=gδ(r−r′), (1.4) with the coupling constant $$g=4\pi \hbar^2 a/m$$. The substitutions of the effective delta-potential and the mean-field approximation into (1.1) lead to the following nonlinear Schrödinger equation or Gross–Pitaevskii equation   iℏ∂tΨ(r,t)=[−ℏ22m∇2+Vext(r)+g|Ψ(r,t)|2]Ψ(r,t). (1.5) As pointed out in Dalfovo et al. (1999), the validity of (1.5) is based on the following: (I) the $$s$$-wave scattering length $$a$$ should be much smaller than the average distance between atoms; hence the effective delta potential (1.4) can be used to replace the interatomic potential; (II) the number of atoms in the condensate should be much larger than 1; then the classical mean field given in (1.3) can be used to substitute the field operator approximately. For time-independent external potentials, the NLSE or GP model admits two important conserved quantities: (C1) The boson gas is only investigated nonrelativistically, and no real particles can be created or annihilated. Thus the total number of atoms (with a suitable normalization)   N=∫|Ψ(r,t)|2dr (1.6) should be conserved in the dynamics; (C2) The gas system is an isolated system, and thus the total energy   E=∫[ℏ22m|∇Ψ(r,t)|2+Vext|Ψ(r,t)|2+12g|Ψ(r,t)|4]dr (1.7) should also be conserved. We consider the $$1+1$$-dimensional NLSE for the macroscopic wave function   iℏ∂tΨ=−ℏ22m∂xxΨ+g|Ψ|2Ψ+VextΨ, (1.8) where $$ V_{{\rm ext}}=V_{{\rm ext}}(x)$$ is the external potential. For the BEC species (Dalfovo et al., 1999), $$a$$ may take either positive or negative values depending on repulsive or attractive interactions between the particles, respectively. With the help of the transformations,   t~=g2t,x~=|g|xandϕ(x~,t~)=Ψ(x,t)|g|=Ψ(x~|g|,t~g2)|g|, the parameter $$g$$ can be removed up to a sign. We will use atomic units $$m=\hbar=1$$ throughout the whole text and fix $$g=\pm 1$$ corresponding to the case of a repulsive and an attractive nonlinearity, respectively, except other specified. With regard to the case of $$V_{{\rm ext}}\equiv 0$$, there has been increasing interest not only on analytically solving this problem and its variants in the physical and mathematical areas (see, e.g., Haus & Wong, 1996; Kivshar & Luther-Davies, 1998; Witthaut et al., 2005, 2009; Carretero-González et al., 2008; Rapedius & Korsch, 2009) but also on various numerical treatments including symplectic or multi-symplectic methods (see, e.g., Reich, 2000; Islas et al., 2001; Bridges & Reich, 2006; Hong et al., 2007; Kong et al., 2015 and references therein). To reformulate this problem in the form of a multi-symplectic Hamiltonian system (MSHS), by separating the real and imaginary parts of the wave function $${\it{\Psi}}=q+{\bf i}p$$ and introducing a pair of conjugate momenta, $$v=\partial_x q$$ and $$w=\partial_x p$$, one obtains   M∂tz+K∂xz=∇zS(z), (1.9) where $$z=(q,p,v,w)^T$$ is the state variable; $$M$$ and $$K$$ are skew-symmetric matrices   M =(0−200200000000000),K=(00100001−10000−100), and $$S: \mathbf{R}^4\rightarrow \mathbf{R}$$ is the multi-symplectic Hamiltonian $$S(z)=\frac{1}{2}g(q^2+p^2)^2-\frac{1}{2}(v^2+w^2).$$ We call a system of the form (1.9) an MSHS, because it admits the following multi-symplectic conservation law (MSCL) (Bridges, 1997; Reich, 2000; Bridges & Reich, 2006):   ∂tω+∂xκ=0, (1.10) where $$\omega$$ and $$\kappa$$ are the pre-symplectic forms defined by   ω=12dz∧Mdzandκ=12dz∧Kdz. (1.11) Besides the MSCL, every MSHS has other two local conservation laws, i.e., (a) the energy conservation law (ECL)   ∂tE+∂xF=0, (1.12) with the energy density $$E=S(z)-\frac{1}{2}z^TK\partial_x z$$ and the energy flux density $$F=\frac{1}{2}z^TK\partial_t z$$ and (b) the momentum conservation law (MCL)   ∂tI+∂xG=0, (1.13) with the momentum density $$I=\frac{1}{2}z^TM\partial_x z$$ and the momentum flux density $$G=S(z)-\frac{1}{2}z^TM\partial_t z$$. The MSHS (1.9), originating from the nonlinear Schrödinger equation, has a more conservate law and it is referred to as the charge conservation law, or the particle-number conservation law, the mass conservation law in various physical problems:   ∂tNd+∂xNf=0, (1.14) where the density $$N_d=\frac{1}{4}(Mz)^TMz$$ and the flux density $$N_f=\frac{1}{2}(Mz)^TKz$$. No matter what conservation law is called, it commonly describes some physical quantity of finite value and can be normalized to unity. Hereafter we call it the local normalization conservation law (NCL). Under appropriate boundary conditions, e.g., zero or periodic boundary conditions, the local conservations laws (1.10), (1.12), (1.13) and (1.14) imply the corresponding global ones, which are obtained by integrating the local ones over the spatial domain and then eliminating the flux terms by suitably imposing boundary conditions (see, e.g., Bridges & Reich, 2006; Hong & Li, 2006): (I) The global symplecticity in time   ddt∫ω(z(x,t))dx=0. (1.15) (II) The global ECL   ddtE≜ddt∫E(z(x,t))dx=0. (1.16) (III) The global MCL   ddtI≜ddt∫I(z(x,t))dx=0. (1.17) (IV) The global NCL   ddtN≜ddt∫Nd(z(x,t))dx=0. (1.18) We should point out here that $$\mathscr{N}$$ appears in either (1.6) or (1.18), which is essentially the same. The expressions of $$\mathscr{E}$$ appearing in (1.7) and (1.16), as shown in Proposition 2.10, are equivalent under suitable boundary conditions. Concerning the case $$V_{{\rm ext}}\not\equiv 0$$ for (1.8), it does not seem to be very optimistic from an MS geometric viewpoint. To the best of our knowledge, published work about the MSHS reformulation of (1.8) and relevant MS geometric numerical methods applied to this model in the literature are rare especially for the case $$V_{{\rm ext}}$$ explicitly depending on $$x$$ or $$t$$. In this article, we focus on the equation (1.8) with finite external delta potentials $$V_{{\rm ext}}=\sum\limits_{\mu=1}^{J}\gamma_{\mu}\delta(x-x_{\mu}^{*})$$, where $$J$$ is a positive integer and $$x_{\mu}^{*}$$ ($$\mu=1,\cdots,J$$) are all given points with the assumption that $$x_1^{*}<x_2^{*}<\cdots<x_J^{*}$$. More precisely, the following NLSE   i∂tΨ=−12∂xxΨ+g|Ψ|2Ψ+[∑μ=1Jγμδ(x−xμ∗)]Ψ, (1.19) will be investigated theoretically and numerically from a viewpoint of the MSHS theory. There exist a lot of references investigated on stationary and dynamical delta-potential NLSEs, which have single, double, finite delta and delta-shell potentials, and even infinite ones like the Dirac–Comb potentials (Berezin & Shubin, 1991; Haus & Wong, 1996; Kivshar & Luther-Davies, 1998; Witthaut et al., 2005, 2009; Rapedius & Korsch, 2009; Cartarius et al., 2012). These are imposed by various kinds of boundary conditions, e.g., periodic or zero at infinite and finite distances, absorbing and reflecting, Siegert ones and many others arising in different physical problems. In this context, we just restrict to the first class, namely the infinitely distant periodic or zero ones, which are commonly related to bound-state problems, while others to the resonance state or scattering problems. Due to the presence of the delta potential, the state variables $$z$$ admit jumps and (1.19) cannot, be reformulated as an MSHS in the normal sense. Noting that the continuity equations (e.g., the conservation of mass, the conservation of momentum and the conservation of energy) in hydrodynamics (Landau & Lifshitz, 1987), are valid in the sense of integration and then rewritten in the form of differentiation we make an attempt to reformulate the above equation into an MSHS form in the sense of integration in Section 2. Based on the given weak MSHS form, relevant theoretical discussions including some local and global conservation laws are addressed as well. In Section 3, novel concatenating RK discretizations for (1.19) are constructed, the multi-symplecticity of them are revealed, and some local and global intrinsic conservation laws are investigated theoretically under multi-symplectic Runge–Kutta (MSRK) discretizations. Section 4 is devoted to numerical experiments and extensive relevant discussions. We conclude in Section 5. 2. The weak MSHS reformulation and conservation laws For a function $$f(x,t)$$ which is well defined in some deleted neighborhood of $$(\hat{x},t)$$, we define   fx=x^−=limx→x^−f(x,t)andfx=x^+=limx→x^+f(x,t). In what follows, we simplify $$f^-_{x=x_{\mu}^{*}}=f^-_\mu$$ and $$f^+_{x=x_{\mu}^{*}}=f^+_\mu$$ if there is no confusion. According to the distribution theory (see, Witthaut et al., 2005; Carretero-González et al., 2008 and references therein), the wave function $${\it{\Psi}}$$ in (1.19) is continuous, but its first spatial partial derivative admits the following discontinuities (or jumps) at $$x = x_{\mu}^{*}$$ that   ∂xΨμ+−∂xΨμ−=2γμΨ(xμ∗), because of the presence of the corresponding delta potentials, where $$\mu=1,2,\ldots,J$$. Hereafter, the points $$x= x_{\mu}^{*}$$, as proposed in Bai & Wang (2015) and Bai (2016), are referred to as interface points. Proposition 2.1 The equation (1.19), by introducing the state variable $$z=(q,p,v,w)^T$$ as specified in (1.9), is equivalent to the following system, in which the singular potentials $$\gamma_\mu\delta (x-x_{\mu}^{*})$$ are replaced by some extra jump conditions:   {M∂tz+K∂xz=∇zS(z),x≠xμ∗zμ+−zμ−=2γμ(0,0,qμ−,pμ−)T≜Ajump(μ)zμ−,  (2.1) for $$\mu=1,\ldots,J$$, where $$z^{\pm}_\mu=(q^{\pm}_\mu,p^{\pm}_\mu,v^{\pm}_\mu,w^{\pm}_\mu)^T$$ and the matrix   Ajump(μ)=(000000002γμ00002γμ00). Furthermore, it can be verified from (2.1) that the derivatives of $$v$$ and $$w$$, with respect to $$x$$, namely the corresponding second derivatives of $$q$$ and $$p$$, respectively, are continuous at the interface point $$x=x_{\mu}^{*}$$, i.e.,   ∂xvμ+=∂xvμ−and∂xwμ+=∂xwμ−. (2.2) It is apparent that the system (2.1), at any point $$x\neq x_{\mu}^{*}$$, is equivalent to the MSHS (1.9), and the local conservation laws (1.10), (1.12)–(1.14) are then automatically valid. Thus the local conservation laws are investigated only at the interface points subsequently. 2.1 Local conservation laws Noting the fact that the exterior differential operator $$d$$ can commute with the partial derivatives $$\partial_t$$ and $$\partial_x$$, we first give the variational equations in terms of (2.1)   {M∂tdz+K∂xdz=DzzS(z)dz,x≠xμ∗,dzμ+−dzμ−=Ajump(μ)dzμ−,  (2.3) for $$\mu=1,\ldots,J$$. Remark 2.2 The local linearizations, which are namely the differential 1-forms from the viewpoint of differential geometry (Bridges, 1997; Arnold, 1999), should have the same properties as the original state variables possess. Illustratively speaking, $$dq^+_\mu=dq^-_\mu$$ come from the fact that the local linearizations $$dq_\mu$$ should be continuous at the interface points $$x=x_{\mu}^{*}$$ as $$q$$ admits, and so $$dp^+_\mu=dp^-_\mu$$; $$dv^+_\mu=dv^-_\mu + 2\gamma dq^-_\mu$$ originate from the fact that $$dv$$ should have the same jumps at $$x=x_{\mu}^{*}$$ as $$v$$ has, and $$dw^+_\mu=dw^-_\mu +2\gamma dp^-_\mu$$. It is worth revisiting the deduction of the well-known Euler–Lagrange equations (Arnold, 1999; Hairer et al., 2006); the variation of the solution, in the action functional, is set to be vanished at endpoints due to the fixed boundary values. Thus it is necessary and natural to combine the discontinuities into variational demonstrations here. To be distinguished, the phrases ‘in a weak sense’ or ‘in the weak sense’ will be used in the following assertions. We say a conservation law $$\partial_t C_d+\partial_x C_f=0$$ holds in the weak sense at some point $$(\hat{x},\hat{t})$$, it means for any small enough rectangle space-time interval $$[x_1,x_2]\times[t_1,t_2]$$ with $$\hat{x}\in (x_1,x_2)$$ and $$\hat{t} \in (t_1,t_2)$$ such that the following integration   ∫x1x2(Cd(x,t2)−Cd(x,t1))dx+∫t1t2(Cf(x2,t)−Cf(x1,t))dt=0, (2.4) holds. We always check, in the weak sense, the validity of a conservation law through the above relation. Theorem 2.3 The system (2.1), in the weak sense we clarified above, admits the local MSCL (1.10) at each interface point $$x=x_{\mu}^{*}$$. Proof. Given any small enough rectangle interval $$[x_1,x_2]\times [t_1,t_2]$$ with $$x_\mu^{*}\in (x_1,x_2)$$ and excluding any other interface point in $$[x_1,x_2]$$, we will check the validity of the MSCL (1.10) in the weak sense. As suggested in (2.4), we check    ∫x1x2[ω(z(x,t2))−ω(z(x,t1))]dx+∫t1t2[κ(z(x2,t))−κ(z(x1,t))]dt =∫x1xμ∗−ϵ[ω(z(x,t2))−ω(z(x,t1))]dx+∫t1t2[κ(z(xμ∗−ϵ,t))−κ(z(x1,t))]dt +∫xμ∗−ϵxμ∗+η[ω(z(x,t2))−ω(z(x,t1))]dx+∫t1t2[κ(z(xμ∗+η,t))−κ(z(xμ∗−ϵ,t))]dt +∫xμ∗+ηx2[ω(z(x,t2))−ω(z(x,t1))]dx+∫t1t2[κ(z(x2,t))−κ(z(xμ∗+η,t))]dt, where $$\epsilon$$ and $$\eta$$ are small enough positive real numbers. Noting the fact that the MSCL (1.10) is valid in the normal sense at any point $$x\neq x_\mu^{*}$$, the first and third lines of the right-hand side (RHS) of the above are automatically vanishing, i.e.,    ∫x1x2[ω(z(x,t2))−ω(z(x,t1))]dx+∫t1t2[κ(z(x2,t))−κ(z(x1,t))]dt =2∫xμ∗−ϵxμ∗+η(dp∧dq)|t1t2dx+∫t1t2(dq∧dv+dp∧dw)|xμ∗−ϵxμ∗+ηdt. (2.5) Let us consider the limiting case as $$\epsilon$$ and $$\eta$$ tend to zero positively. With the aid of the jumps given in the second equation of (2.3), the first term in the RHS of (2.5) vanishes because of the continuity of $$dp$$ and $$dq$$ at the interface point $$x=x_\mu^{*}$$, and the second term becomes    ∫t1t2(dqμ+∧dvμ++dpμ+∧dwμ+−dqμ−∧dvμ−−dpμ−∧dwμ−)dt =∫t1t2(dqμ−∧(dvμ+−dvμ−)+dpμ−∧(dwμ+−dwμ−))dt =∫t1t2(dqμ−∧(2γμdqμ−)+dpμ−∧(2γμdpμ−))dt=0, (2.6) where the first equality uses the fact that $$dq^+_\mu =dq^-_\mu $$ and $$dp^{+}_\mu =dp^{-}_\mu$$. Thus (2.5) and (2.6) give rise to this theorem. □ Definition 2.4 For the purpose of being distinguished, a system of the abstract form (1.9) is called a normal MSHS, and a system of the form (2.1) is called a weak MSHS since it preserves the MSCL (1.10) in the weak sense at each interface point $$x=x_\mu^{*}$$. And conservation laws, which are proved valid in the weak sense, are also referred to with the term ‘weak’. Theorem 2.5 The system (2.1) admits the weak local NCL (1.14) at each interface point $$x=x_\mu^{*}$$. Proof. As shown in Theorem 2.3, it suffices to prove that   ∫xμ∗−ϵxμ∗+η[Nd(z(x,t2))−Nd(z(x,t1))]dx+∫t1t2[Nf(z(xμ∗+η,t))−Nf(z(xμ∗−ϵ,t))]dt=0, i.e.,   ∫xμ∗−ϵxμ∗+η(p2+q2)|t1t2dx+∫t1t2(qw−pv)|xμ∗−ϵxμ∗+ηdt=0 (2.7) with the positive parameters $$\epsilon$$ and $$\eta$$ tending to zero. The left-hand side (LHS) of 2.7, as $$\epsilon$$ and $$\eta$$ tend to zero positively, yields    ∫t1t2(qμ+wμ+−pμ+vμ+−qμ−wμ−+pμ−vμ−)dt=∫t1t2(qμ−(2γμpμ−)−pμ−(2γμqμ−))dt=0. (2.8) This completes the proof. □ Theorem 2.6 The system (2.1) admits the weak local ECL (1.12) at each interface point $$x=x_\mu^{*}$$. Proof. As discussed in Theorems 2.3 and 2.5, with the help of (2.1), we directly implement the following calculations    ∫xμ∗−ϵxμ∗+η[E(z(x,t2))−E(z(x,t1))]dx+∫t1t2[F(z(xμ∗+η,t))−F(z(xμ∗−ϵ,t))]dt =12∫xμ∗−ϵxμ∗+η(g(p2+q2)2−(q∂xv+p∂xw))|t1t2dx+12∫t1t2(q∂tv+p∂tw−v∂tq−w∂tp)|xμ∗−ϵxμ∗+ηdt. Let the positive parameters $$\epsilon$$ and $$\eta$$ tend to zero the first term of the right-hand side (RHS) of the above vanishes because of the fact that $$p$$ and $$q$$ are continuous at $$x=x_\mu^{*}$$ and also the same property for $$\partial_x\! v$$ and $$\partial_x\! w$$ as specified in (2.2). Also with $$\partial_t q^+_\mu=\partial_t q^-_\mu$$, $$\partial_t p^+_\mu=\partial_t p^-_\mu$$, $$\partial_t v^+_\mu-\partial_t v^-_\mu=2\gamma_{\mu}\partial_t q_\mu^-$$ and $$\partial_t w^+_\mu-\partial_t w^-_\mu=2\gamma_{\mu}\partial_t p_\mu^-$$, the second term reads   12∫t1t2(qμ−(∂tvμ+−∂tvμ−)+pμ−(∂twμ+−∂twμ−)+∂tqμ−(vμ−−vμ+)+∂tpμ−(wμ−−wμ+))dt =12∫t1t2(qμ−(2γμ∂tqμ−)+pμ−(2γμ∂tpμ−)−∂tqμ−(2γμqμ−)−∂tpμ−(2γμpμ−))dt=0. (2.9) Hence it completes the proof. □ Remark 2.7 The system (2.1), however, does not admit the weak local MCL (1.13) at each interface point. According to the Nöther theorem (Olver, 1993; Arnold, 1999), a continuous conservation law corresponds to a continuous variational symmetry. The time translation invariance implies the local ECL, and the space translation invariance implies the local MCL. Comparing (1.19) with (2.1), the delta potentials do not break the time translation invariance, and thus the local ECL is still valid. The space translation invariance, however, is broken at each interface point due to the extra jump conditions. 2.2 Global conservation laws It is evident that the local conservation laws are independent of whatever initial-boundary value conditions imposed. But it changes for the corresponding global ones. And the well posedness of the considered NLSEs need suitable initial and boundary conditions either. Throughout this context, we assume that the solution is smooth enough, except for the interface points $$x=x_\mu^{*}$$ for $$\mu=1,\ldots,J$$, the spatial interval we investigate is $$[x_L, x_R]$$ with all $$x_\mu^{*}\in (x_L, x_R)$$. Here $$x_L$$ and $$x_R$$ could be infinities in principle. To proceed, we supplement the system (2.1) by prescribing suitable conditions with:   {z|t=0=z0(x)≡(q0(x),p0(x),v0(x),w0(x))T,z(xL,t) = z(xR,t),∂xz(xL,t)=∂xz(xR,t),  (2.10) where $$(q_0(x),p_0(x), v_0(x), w_0(x))$$ satisfy the jump conditions as given in (2.1) and $$t\geq 0$$ for the second and third equalities. Theorem 2.8 Under the conditions prescribed in (2.10), the system (2.1) has the weak global symplecticity in time (1.15) with the integration interval $$[x_L, x_R]$$. Proof. As specified in Theorem 2.3, the local MSCL (1.10) is valid at each interface point $$x=x_\mu^{*}$$ in the weak sense, i.e., in the sense of integration. At any other point, the local MSCL is normally valid. Thus (1.10) can be taken as the integration over the spatial interval $$[x_L,x_R]$$ for both its left- and right-hand sides, i.e.,   ∫xLxR(∂ω∂t+∂κ∂x)dx=0. (2.11) As discussed for (2.4) and in Theorem 2.3, the second term in the LHS becomes   ∫xLxR∂κ∂xdx=κ(z(xR,t))−κ(z(xL,t))=0, (2.12) where the second equality is due to the periodic boundary conditions given in (2.10). Thus we have   ∫xLxR∂ω∂tdx=ddt∫xLxRωdx=0, where the first equality is due to the fact that the integration can commute with the temporal derivative. Hence the global symplecticity is conserved for all times,   ∫xLxRω(z(x,t))dx=∫xLxRω(z(x,0))dx=∫xLxRω(z0(x))dx,∀t≥0. (2.13) □ Similarly, we have the following statement about the the global ECL and NCL. Theorem 2.9 Under the assumptions of Theorem 2.8, the system (2.1) admits the weak global ECL (1.16) and the weak global NCL (1.18). As mentioned previously, the expressions of the total energy $$\mathscr{E}$$ given in (1.7) and (1.16) are not the same. To avoid the ambiguity, we give the following proposition. Proposition 2.10 Under the assumptions of Theorem 2.8, the total energies $$\mathscr{E}$$ appearing in (1.7) and (1.16), for the NLSE with the external delta potentials, are equivalent. Proof. By introducing the state variable $$z$$, the total energy appearing in (1.7), which is denoted by $$\mathscr{E}_1$$ to avoid confusion, reads   E1=∫xLxR[12(v2+w2)+∑μ=1Jγμδ(x−xμ∗)(q2+p2)+g2(q2+p2)2]dx. (2.14) With the aid of properties of the delta functions, we immediately get   E1 =[∫xLx1∗−ϵ1+∑μ=1J−1(∫xμ∗−ϵμxμ∗+ημ+∫xμ∗+ημxμ+1∗−ϵμ+1)+∫xJ∗−ϵJxJ∗+ηJ+∫xJ∗+ηJxR]12(v2+w2)dx +∑μ=1Jγμ[q2(xμ∗)+p2(xμ∗)]+g2∫xLxR(q2+p2)2dx, (2.15) where we suppress the time variable $$t$$, i.e., $$q(x_\mu^*)\equiv q(x_\mu^*,t)$$ and $$p(x_\mu^*)\equiv p(x_\mu^*,t)$$, $$\epsilon_\mu$$ and $$\eta_\mu$$ are small enough positive real numbers. Let $$\epsilon_\mu$$ and $$\eta_\mu$$ tend to zero, take the integration by parts and recall the continuities of $$q$$, $$p$$, $$\partial_x v$$ and $$\partial_x w$$ in the definite domain, we find the first line in the RHS of (2.15) yields    −12∫xLxR(q∂xv+p∂xw)dx−12∑μ=1J{qμ−[vμ+−vμ−]+pμ−[wμ+−wμ−]} =−12∫xLxR(q∂xv+p∂xw)dx−∑μ=1Jγμ[q2(xμ∗)+p2(xμ∗)]. (2.16) Combining (2.15) with (2.16) gives $$\mathscr{E}_1=\mathscr{E}$$. □ Hereafter, we only use the expression of the total energy as given in (1.16) while encountering. 3. Multi-symplectic Runge–Kutta methods In this section, we will introduce novel concatenating RK discretizations applied to the weak MSHS (2.1) and discuss the multi-symplecticity of the numerical methods. To simplify the presentation, we just limit our discussion to the case on a single external delta potential, i.e., $$J=1$$ in (1.19). We remark here that, it can be naturally extended to cases on finite ones. For convenience, we introduce some notions of MS schemes applied to a normal MSHS of the abstract form (1.9). Given a partition over the definite domain, one considers a certain discretization for an MSHS (1.9)   M∂ti,lzil+K∂xi,lzil=(∇zS(z))il, (3.1) at a mesh point $$(x_i,t_l)\in R^2$$, where $$\partial_t^{i,l}$$, $$\partial_x^{i,l}$$ are the corresponding discretizations of the partial derivatives $$\partial_t$$ and $$\partial_x$$, respectively, $$(\nabla_z S(z))_i^l=(\nabla_z S(z_i^l))$$, $$z_i^l$$ is the numerical value of $$z(x,t)$$. Definition 3.1 (Reich, 2000; Bridges & Reich, 2006) The numerical scheme (3.1) of (1.9) is said to be multi-symplectic if   ∂ti,lωil+∂xi,lκil=0, (3.2) still holds at any mesh point $$(x_i,t_l)$$, where   ωil=12(dz)il∧M(dz)ilandκil=12(dz)il∧K(dz)il. As far as we have known, the MS concatenating RK discretizations originated from Reich (2000), and who showed concatenating Gauss–Legendre schemes applied to a $$1+1$$-dimensional MSHS that leads to an MS integrator. Since then there are quite a few works about MSRK-type methods (see, e.g., Hong et al., 2005; Ryland & McLachlan, 2008; McLachlan et al., 2014 and references therein) and their applications to various special MSHSs (see, e.g., Bridges & Reich, 2006; Hong & Li, 2006; Hong et al., 2007 and references therein). Almost all the works, however, are restricted to normal MSHSs. Regarding the weak MSHS (2.1) reformulated here, we make an attempt to introduce the corresponding MS concatenating RK discretizations. In our numeric, the spatial and temporal intervals are taken to be $$[x_L,x_R]$$ and $$[T_{0},T_{1}]$$, respectively. Here and hereafter, we will suppress the index $$1$$, e.g., $$f^{\pm}=f^{\pm}_{1}$$, $$A_{\rm jump}=A_{\rm jump}^{(1)}$$ and so forth. As proposed for the normal MS numerical treatments, we begin with a uniform partition on the definite domain with the mesh points $$(x_i,t_l)$$ and the stepsizes $${\it{\Delta}} x=\frac{x_R-x_L}{N}$$ and $${\it{\Delta}} t=\frac{T_1-T_0}{\hat{N}}$$, where $$x_i=x_L+i {\it{\Delta}} x$$ for $$i=0,1,\ldots,N$$, $$t_l=T_0+l {\it{\Delta}} t$$ for $$l=0,1,\ldots,\hat{N}$$. Assume that the interface point $$x=x^{*}$$ lies in the subinterval $$[x_I,x_{I+1})$$ for some $$I$$ ($$0<I<N$$), we then supplement the interface point as a new grid point labeled by $$x_{I^{'}}$$. Thus the spatial grid points become $$x_L=x_0<x_1<\cdots<x_I\leq x_{I^{'}}=x^{*}<x_{I+1}<\cdots<x_N=x_R$$, and the mesh-lengths of the two new subintervals are denoted by $${\it{\Delta}} x_1$$ ($$\equiv x_{I^{'}}-x_I$$) and $${\it{\Delta}} x_2$$ ($$\equiv x_{I+1}-x_{I^{'}}$$), respectively. In case $${\it{\Delta}} x_1=0$$, which means that the interface point $$x=x^{*}$$ is a grid point after the first equidistant partition, i.e., the grid point $$x_{I}$$ coincides with $$x_{I^{'}}$$ and no further partition is needed. We are now in a position to perform numerical discretization, which is demonstrated by applying $$s$$-stage and $$r$$-stage RK methods to the system (2.1) to the $$x$$- and $$t$$-directions, respectively, on space-time grids. For a starting point $$(x_i,t_l)$$ with $$i\notin \{I,I^{'}\}$$, we have   {Zi,ml,k=zil,k+Δx∑n=1sa~mn∂xZi,nl,k,zi+1l,k=zil,k+Δx∑m=1sb~m∂xZi,ml,k,Zi,ml,k=zi,ml+Δt∑j=1rakj∂tZi,ml,j,zi,ml+1=zi,ml+Δt∑k=1rbk∂tZi,ml,k,M∂tZi,ml,k+K∂xZi,ml,k=∇zS(Zi,ml,k).  (3.3) Here we will specify some notations which first appeared: $$Z_{i,m}^{l,k}\approx z(x_i+\tilde{c}_m{\it{\Delta}} x,t_l+c_k{\it{\Delta}} t)$$, $$z_{i,m}^{l}\approx z(x_i+\tilde{c}_m{\it{\Delta}} x,t_{l})$$, $$z_{i,m}^{l+1}\approx z(x_i+\tilde{c}_m{\it{\Delta}} x,t_{l+1})$$, $$\partial_t Z_{i,m}^{l,k}\approx \partial_t z(x_i+\tilde{c}_m{\it{\Delta}} x,t_l+c_k{\it{\Delta}} t)$$, $$\partial_x Z_{i,m}^{l,k}\approx \partial_x z(x_i+\tilde{c}_m{\it{\Delta}} x,t_l+c_k{\it{\Delta}} t)$$, $$z_{i}^{l,k}\approx z(x_{i},t_l+c_k{\it{\Delta}} t)$$ and $$z_{i+1}^{l,k}\approx z(x_{i+1},t_l+c_k{\it{\Delta}} t)$$ with $$\tilde{c}_m=\sum\limits_{n=1}^s \tilde{a}_{mn}$$ and $$c_k=\sum\limits_{j=1}^r a_{kj}$$ as required (see, e.g., Hairer et al., 2006). And for the starting points $$(x_I,t_l)$$ and $$(x_{I^{'}},t_l)$$, the schemes read, respectively,   {ZI,ml,k=zIl,k+Δx1∑n=1sa~mn∂xZI,nl,k,zI ′,−l,k=zIl,k+Δx1∑m=1sb~m∂xZI,ml,k,ZI,ml,k=zI,ml+Δt∑j=1rakj∂tZI,ml,j,zI,ml+1=zI,ml+Δt∑k=1rbk∂tZI,ml,k,M∂tZI,ml,k+K∂xZI,ml,k=∇zS(ZI,ml,k)  (3.3B) and   {ZI ′,ml,k=zI ′,+l,k+Δx2∑n=1sa~mn∂xZI ′,nl,k,zI+1l,k=zI ′,+l,k+Δx2∑m=1sb~m∂xZI ′,ml,k,ZI ′,ml,k=zI ′,ml+Δt∑j=1rakj∂tZI ′,ml,j,zI ′,ml+1=zI ′,ml+Δt∑k=1rbk∂tZI ′,ml,k,M∂tZI ′,ml,k+K∂xZI ′,ml,k=∇zS(ZI ′,ml,k),  (3.3C) where $$Z_{I,m}^{l,k}\approx z(x_{I}+\tilde{c}_m {\it{\Delta}} x_1,t_l+c_k{\it{\Delta}} t)$$, $$Z_{I^{'}\!\!,m}^{l,k}\approx z(x_{I^{'}}+\tilde{c}_m {\it{\Delta}} x_2,t_l+c_k{\it{\Delta}} t)$$, $$z_{I^{'}\!\!,\pm}^{l,k}\approx z^{\pm}(x_{I^{'}},t_l+c_k{\it{\Delta}} t)$$, etc. And the jump conditions at $$x_{I^{'}}$$ are subjected to   zI ′,+l,k−zI ′,−l,k=AjumpzI ′,−l,k. (3.4) The variational equations corresponding to (3.3) are given by   {dZi,ml,k=dzil,k+Δx∑n=1sa~mnd(∂xZ)i,nl,k,dzi+1l,k=dzil,k+Δx∑m=1sb~md(∂xZ)i,ml,k,dZi,ml,k=dzi,ml+Δt∑j=1rakjd(∂tZ)i,ml,j,dzi,ml+1=dzi,ml+Δt∑k=1rbkd(∂tZ)i,ml,k,Md(∂tZ)i,ml,k+Kd(∂xZ)i,ml,k=DzzS(Zi,ml,k)dZi,ml,k,  (3.5) where $$i\notin \{I,I^{'}\}$$ and the abbreviations, i.e., $$dZ_{i,m}^{l,k}$$ for $$(dZ)_{i,m}^{l,k}$$, $$d(\partial_t Z)_{i,m}^{l,j}$$ for $$(d(\partial_t Z))_{i,m}^{l,j}$$ and so on, are utilized for the sake of simplicity. The variational equations for ($$3.3$$B) and ($$3.3$$C) are similar to the above, and we would not exhibit them here anymore. And the jump conditions for the variational equations corresponding to (3.4), yield   dzI ′,+l,k−dzI ′,−l,k=AjumpdzI ′,−l,k. (3.6) Subsequently, we mention the numerical algorithm (3.3) which means (3.3) together with ($$3.3$$B) and ($$3.3$$C). With the above settings, we have the following theorem on the multi-symplecticity of the concatenating scheme. Since the deductions are similar to that given in Reich (2000) and Hong et al. (2005), we just omit the proof here. Theorem 3.2 If the two RK methods we applied satisfy the following coefficient conditions   {bkbj−bkakj−bjajk=0,b~mb~n−b~ma~mn−b~na~nm=0,  (3.7) for all $$k,j=1,2,\ldots,r$$ and $$m,n=1,\ldots,s$$, then the numerical algorithm (3.3) is multi-symplectic with the following discrete MSCL: (a) For $$i\notin \{I,I^{'}\}$$,   1Δt∑m=1sb~m(ωi,ml+1−ωi,ml)+1Δx∑k=1rbk(κi+1l,k−κil,k)=0, (3.8) where $${\it{\omega}}_{i,m}^{l+1}=\frac{1}{2}{\,{\rm d}}z_{i,m}^{l+1}\wedge M{\,{\rm d}}z_{i,m}^{l+1}$$, $$\kappa_{i}^{l,k}=\frac{1}{2}{\,{\rm d}}z_{i}^{l,k}\wedge K{\,{\rm d}}z_{i}^{l,k}$$, and so forth. (b) For $$i= I$$,   1Δt∑m=1sb~m(ωI,ml+1−ωI,ml)+1Δx1∑k=1rbk(κI ′,−l,k−κIl,k)=0. (3.9) (c) For $$i= I^{'}$$,   1Δt∑m=1sb~m(ωI ′,ml+1−ωI ′,ml)+1Δx2∑k=1rbk(κI+1l,k−κI ′,+l,k)=0. (3.10) Taking a summation of (3.8)–(3.10) for all the desired indexes $$i$$, we have   [Δx∑i∉{I,I ′}∑m=1sb~m(ωi,ml+1−ωi,ml)+Δx1∑m=1sb~m(ωI,ml+1−ωI,ml)+Δx2∑m=1sb~m(ωI ′,ml+1−ωI ′,ml)]+[Δt∑i∉{I,I ′}∑k=1rbk(κi+1l,k−κil,k)+Δt∑k=1rbk(κI ′,−l,k−κIl,k)+Δt∑k=1rbk(κI+1l,k−κI ′,+l,k)]=0. Using the periodic conditions prescribed previously and taking the summation associatively, the terms in the above second square brackets can be simplified to   Δt∑k=1rbk(κI ′,−l,k−κI ′,+l,k) =Δt2∑k=1rbk(dqI ′,−l,k∧dvI ′,−l,k+dpI ′,−l,k∧dwI ′,−l,k−dqI ′,+l,k∧dvI ′,+l,k−dpI ′,+l,k∧dwI ′,+l,k)=0, where the last equality is similar to (2.6) with the jump condition (3.6) utilized. Thus, we have the following statement. Theorem 3.3 Under the assumptions of Theorem 3.2, the numerical scheme (3.3) preserves the discrete global symplecticity in time   Δx∑i∉{I,I ′}∑m=1sb~mωi,ml+1+Δx1∑m=1sb~mωI,ml+1+Δx2∑m=1sb~mωI ′,ml+1 =Δx∑i∉{I,I ′}∑m=1sb~mωi,ml+Δx1∑m=1sb~mωI,ml+Δx2∑m=1sb~mωI ′,ml, (3.11) with the periodic boundary conditions, i.e., the last two lines of (2.10) imposed. As discussed in Hong et al. (2005) and commented in Hong & Li (2006), two symplectic RK discretizations applied to a normal MSHS in time and space, respectively, leads to a multi-symplectic integrator. Theorem 3.2 tells us that a similar assertion holds for this weak MSHS (2.1). And (2.1) possesses not only this MSCL but also some local and global conservation laws. It is well known in numerical computation that one remarkably important standard in constructing algorithms is the ability of carrying as much of qualitative solution behavior as possible or as many of intrinsic characters as possible. The local and global NCLs, as we state next, are precisely carried under our novelly constructed MSRK discretizations. Theorem 3.4 Under the assumptions of Theorem 3.2, the numerical scheme (3.3) preserves the discrete NCL: (a) For $$i\notin \{I,I^{'}\}$$,   1Δt∑m=1sb~m(Nd(zi,ml+1)−Nd(zi,ml))+1Δx∑k=1rbk(Nf(zi+1l,k)−Nf(zil,k))=0. (3.12) (b) For $$i= I$$,   1Δt∑m=1sb~m(Nd(zI,ml+1)−Nd(zI,ml))+1Δx1∑k=1rbk(Nf(zI ′,−l,k)−Nf(zIl,k))=0. (3.13) (c) For $$i= I^{'}$$,   1Δt∑m=1sb~m(Nd(zI ′,ml+1)−Nd(zI ′,ml))+1Δx2∑k=1rbk(Nf(zI+1l,k)−Nf(zI ′,+l,k))=0. (3.14) Proof. We start with the following calculations   Nd(zi,ml+1)−Nd(zi,ml)=14[(Mzi,ml+1)TMzi,ml+1−(Mzi,ml)TMzi,ml] =14[Δt∑k=1rbk(Mzi,ml)TM∂tZi,ml,k+Δt∑k=1rbk(M∂tZi,ml,k)TMzi,ml  +Δt2∑k=1r∑j=1rbkbj(M∂tZi,ml,k)TM∂tZi,ml,j], (3.15) where the second equality is due to the second equation of (3.3). Making use of the first equation of (3.3), the first term in the right-hand side of (3.15) yields   Δt∑k=1rbk(Mzi,ml)TM∂tZi,ml,k =Δt∑k=1rbk(MZi,ml,k)TM∂tZi,ml,k−Δt2∑k=1r∑j=1rbkakj(M∂tZi,ml,j)TM∂tZi,ml,k. (3.16) Analogously, the second term reads   Δt∑k=1rbk(M∂tZi,ml,k)TMzi,ml =Δt∑k=1rbk(M∂tZi,ml,k)TMZi,ml,k−Δt2∑k=1r∑j=1rbkakj(M∂tZi,ml,k)TM∂tZi,ml,j =Δt∑k=1rbk(MZi,ml,k)TM∂tZi,ml,k−Δt2∑j=1r∑k=1rbjajk(M∂tZi,ml,j)TM∂tZi,ml,k, (3.17) where the first part in the right-hand side is due to the symmetry of the scalar product of vectors and the second the interchange of the summation indexes $$k$$ and $$j$$. By virtue of the first equation given in the MS conditions (3.7), (3.15)–(3.17), this tells us that   1Δt∑m=1sb~m(Nd(zi,ml+1)−Nd(zi,ml))=12∑m=1s∑k=1rb~mbk(MZi,ml,k)TM∂tZi,ml,k. (3.18) As deducted for (3.15)–(3.18), we have   1Δx∑k=1rbk(Nf(zi+1l,k)−Nf(zil,k))=12∑k=1r∑m=1sbkb~m(MZi,ml,k)TK∂xZi,ml,k. (3.19) Combining (3.18) and (3.19) with the last equation of (3.3) yields   1Δt∑m=1sb~m(Nd(zi,ml+1)−Nd(zi,ml))+1Δx∑k=1rbk(Nf(zi+1l,k)−Nf(zil,k))=−12∑k=1r∑m=1sbkb~m(MZi,ml,k)T∇zS(Zi,ml,k)=0, where the second equality is due to the fact that $$(Mz)^T\,\nabla_z S(z)\equiv 0$$. This completes the statement (a). The discussions for the last two statements are similar, and thus are omitted here. □ Theorem 3.5 Under the assumptions of Theorem 3.3, the numerical scheme (3.3) preserves the discrete global NCL   Δx∑i∉{I,I ′}∑m=1sb~mNd(zi,ml+1)+Δx1∑m=1sb~mNd(zI,ml+1)+Δx2∑m=1sb~mNd(zI ′,ml+1) =Δx∑i∉{I,I ′}∑m=1sb~mNd(zi,ml)+Δx1∑m=1sb~mNd(zI,ml)+Δx2∑m=1sb~mNd(zI ′,ml). (3.20) Proof. As discussed for Theorem 3.3, taking a product of $${\it{\Delta}} x {\it{\Delta}} t$$, $${\it{\Delta}} x_1 {\it{\Delta}} t$$ and $${\it{\Delta}} x_2 {\it{\Delta}} t$$ with both sides of (3.12)–(3.14), respectively, and then summing them up, we have   [Δx∑i∉{I,I′}∑m=1sb~m(Nd(zi,ml+1)−Nd(zi,ml))  +Δx1∑m=1sb~m(Nd(zI,ml+1)−Nd(zI,ml))+Δx2∑m=1sb~m(Nd(zI′,ml+1+1)−Nd(zI′,ml))] +[Δt∑i∉{I,I′}∑m=1sbk(Nf(zi+1l,k)−Nf(zil,k))  +Δt∑k=1rbk(Nf(zI′,−l,k)−Nf(zIl,k))+Δt∑k=1rbk(Nf(zI+1l,k)−Nf(zI′,+l,k))=0.] By applying periodic boundary conditions, one reduces the part in the above second square brackets to   Δt∑k=1rbk(Nf(zI′,−l,k)−Nf(zI′,+l,k)) =Δt∑k=1rbk{(qI′,−l,k,wI′,−l,k−pI′,−l,k,vI′,−l,k)−(qI′,+l,k,wI′,+l,k,pI′,+l,k,vI′,+l,k,)}=0, where the last equality can be readily verified by the jump conditions (3.4). Thus it leads to the proof. □ The normalization conservation is an exceptionally important, maybe the most important, conserved quantity in quantum physics (Berezin & Shubin, 1991; Dalfovo et al., 1999; Sulem & Sulem, 1999; Carretero-Gonzáalez et al., 2008). Together with its probability interpretation, it provides the fundamental basis for understanding NLSEs physically and studying NLSEs and other intrinsic conservation laws mathematically. It is well-known that symplectic integrators for Hamiltonian ODEs may not preserve the energy exactly (Hairer et al., 2006). Concerning MS integrators for normal MSHSs, similar phenomena (see, e.g., Reich, 2000; Bridges & Reich, 2006; Hong & Li, 2006; Hong et al., 2007) are frequently observed from numerical experiments for nonlinear RHS $$\nabla_z S(z)$$, but a rigorous proof remains to be given. The existing discussions are commonly restricted to the fact that MSRK methods can preserve the energy of linear normal MSHSs (Reich, 2000; Islas et al., 2001; Hong et al., 2005; Bridges & Reich, 2006; Hong & Li, 2006), i.e., $$S(z)$$ is quadratic in $$z$$. As for a weak MSHS of (2.1), we have the following similar statement. Theorem 3.6 For a system of the abstract form (2.1) (restrict to $$J=1$$ for simplicity), provided that (3.7) is satisfied and $$S(z)$$ is quadratic in $$z$$, i.e., $$S(z)=\frac{1}{2}z^THz+\vec{h}^Tz$$ with $$H$$ a symmetric matrix and $$\vec{h}$$ a vector, then the numerical scheme (3.3) preserves the following discrete local ECL: (a) For $$i\notin \{I,I^{'}\}$$,   1Δt∑m=1sb~m(E(zi,ml+1)−E(zi,ml))+1Δx∑k=1rbk(F(zi+1l,k)−F(zil,k))=0. (3.21) (b) For $$i= I$$,   1Δt∑m=1sb~m(E(zI,ml+1)−E(zI,ml))+1Δx1∑k=1rbk(F(zI ′,−l,k)−F(zIl,k))=0. (3.22) (c) For $$i= I^{'}$$,   1Δt∑m=1sb~m(E(zI ′,ml+1)−E(zI ′,ml))+1Δx2∑k=1rbk(F(zI+1l,k)−F(zI ′,+l,k))=0. (3.23) With periodic boundary conditions prescribed in Theorem 3.3, the algorithm has the following discrete global ECL:   Δx∑i∉{I,I ′}∑m=1sb~mE(zi,ml+1)+Δx1∑m=1sb~mE(zI,ml+1)+Δx2∑m=1sb~mE(zI ′,ml+1) =Δx∑i∉{I,I ′}∑m=1sb~mE(zi,ml)+Δx1∑m=1sb~mE(zI,ml)+Δx2∑m=1sb~mE(zI ′,ml). (3.24) 4. Numerical experiments Based on the theoretical results in the previous section, numerical experiments presented in this section reveal numerical characters of the MSRK methods applied to the system (2.1). 4.1 Numerical experiment A: a single-delta potential Here we concern the single-delta potential case that $$V_{\rm ext}(x)=\gamma \delta(x)$$ in (1.8), which can be reformulated as the weak MSHS (2.1) as investigated previously. We discretize the system (2.1) by the MSRK method with $$r=1$$ and $$s=1$$, i.e., $$a_{11}=\tilde{a}_{11}=1/2$$, $$c_1=\tilde{c}_1=1/2$$ and $$b_1=\tilde{b}_1=1$$, which is an implicit second-order MSRK method (in short, denoted by MM2) and is equivalent to the central box scheme (see, e.g., Reich, 2000; Bridges & Reich, 2006; Hong & Li, 2006). For $$i\notin \{I,I^{'}\}$$, it reads   M(zi+12l+1−zi+12l)Δt+K(zi+1l+12−zil+12)Δx=∇zS(zi+12l+12) (4.1) with that   zi+12l=(zi+1l+zil)/2,  zi+12l+1=(zi+1l+1+zil+1)/2,  zil+12=(zil+1+zil)/2,zi+1l+12=(zi+1l+1+zi+1l)/2,  zi+12l+12=(zi+1l+1+zi+1l+zil+1+zil)/4; for $$i=I$$, according to ($$3.3$$B), one can use $$z^{l+\frac{1}{2}}_{I^{'}\!\!,-}$$ and $${\it{\Delta}} x_1$$ to substitute $$z^{l+\frac{1}{2}}_{i+1}$$ and $${\it{\Delta}} x$$ appearing in (4.1), respectively; for $$i=I^{'}$$, according to ($$3.3$$C), one can use $$z^{l+\frac{1}{2}}_{I^{'}\!\!,+}$$, $$z^{l+\frac{1}{2}}_{I+1}$$ and $${\it{\Delta}} x_2$$ to substitute $$z^{l+\frac{1}{2}}_{i}$$, $$z^{l+\frac{1}{2}}_{i+1}$$ and $${\it{\Delta}} x$$, respectively. Here   zI+12l=(zI ′,−l+zIl)/2,zI+12l+12=(zI ′,−l+1+zI ′,−l+zIl+1+zIl)/4,zI ′+12l=(zI ′,+l+zI+1l)/2,zI ′+12l+12=(zI ′,+l+1+zI ′,+l+zI+1l+1+zI+1l)/4 and so on. Relevant to (3.4), the jump conditions for MM2 are specified as   zI ′,+l−zI ′,−l=AjumpzI ′,−landzI ′,+l+1−zI ′,−l+1=AjumpzI ′,−l+1. (4.2) For the purpose of numerical comparisons, we will consider two more implicit second-order nonmulti-symplectic methods here: one is about the symplectic RK method with $$r=1$$ in time and an immersed interface method (see, e.g., Bai & Wang, 2015; Bai, 2016 and references therein) in space (in short, denoted by M2IIM), the other is about the trapezoidal scheme in time and IIM in space (in short, denoted by T2IIM). For the spatial index $$i\notin \{I,I+1\}$$, the scheme M2IIM yields   {2(qil+1−qil)Δt+pi+1l+12−2pil+12+pi−1l+12Δx2=2g[(qil+12)2+(pil+12)2]pil+12,2(pil+1−pil)Δt−qi+1l+12−2qil+12+qi−1l+12Δx2=−2g[(qil+12)2+(pil+12)2]qil+12  (4.3) with $$p^{l+\frac{1}{2}}_{i}=\big(p^{l+1}_{i}+p^{l}_{i}\big)/2$$ processed previously; and for $$i=I$$ and $$i=I+1$$, the scheme reads, respectively,   {2(qIl+1−qIl)Δt+ϕ1(pIl+12)=2g[(qIl+12)2+(pIl+12)2]pIl+12,2(pIl+1−pIl)Δt−ϕ1(qIl+12)=−2g[(qIl+12)2+(pIl+12)2]qIl+12  (4.3B) and   {2(qI+1l+1−qI+1l)Δt+ϕ2(pI+1l+12)=2g[(qI+1l+12)2+(pI+1l+12)2]pI+1l+12,2(pI+1l+1−pI+1l)Δt−ϕ2(qI+1l+12)=−2g[(qI+1l+12)2+(pI+1l+12)2]qI+1l+12  (4.3C) with   {ϕ1(yI)≡βI,1yI−1+βI,2yI+βI,3yI+1ϕ2(yI+1)≡βI+1,1yI+βI+1,2yI+1+βI+1,3yI+2,  (4.4) where the linear coefficients $$\beta_{I,1}$$, $$\beta_{I,2}$$ and $$\beta_{I,3}$$ solve the following linear equations   {βI,1+βI,2+(1+2γxI+1)βI,3=0,xI−1βI,1+xIβI,2+xI+1βI,3=0,xI−12βI,1+xI2βI,2+xI+12βI,3=2,  (4.5) and $$\beta_{I+1,1}$$, $$\beta_{I+1,2}$$ and $$\beta_{I+1,3}$$ solve   { (1−2γxI)βI+1,1+βI+1,2+βI+1,3=0,xIβI+1,1+xI+1βI+1,2+xI+2βI+1,3=0,xI2βI+1,1+xI+12βI+1,2+xI+22βI+1,3=2.  (4.6) Let $$\tilde{z}=(q,p)^T$$ and $$\varphi(\tilde{z})=\big(2g(q^2+p^2)p,-2g(q^2+p^2)q\big)^T$$ here; then the scheme T2IIM is obtained by means of making the substitutions of   φ(z~il+1)+φ(z~il)2,φ(z~Il+1)+φ(z~Il)2andφ(z~I+1l+1)+φ(z~I+1l)2 for the RHSs of (4.3), ($$4.3$$B) and ($$4.3$$C), respectively. It is well known that conservation laws are not only important in physics are but also powerful tools to test and judge numerical algorithms for differential equations, especially for those which could not be solved analytically. To investigate numerically the intrinsic conservation laws addressed before, we need to define, in terms of different numerical schemes, the corresponding discrete analogies of them. Concerning the MS scheme MM2, the discrete local ECL $$(E_{{\rm loc}}^*)_i^l$$, for $$i\notin \{I,I^{'}\}$$ is given by (see Islas et al., 2001; Bridges & Reich, 2006; Hong & Li, 2006)   (Eloc∗)il=1Δt[(S(zi+12l+1)−12(zi+12l+1)TKzi+1l+1−zil+1Δx)−(S(zi+12l)−12(zi+12l)TKzi+1l−zilΔx)]  +1Δx[(zi+1l+12)TKzi+1l+1−zi+1lΔt−(zil+12)TKzil+1−zilΔt]; for $$i=I$$, the expression can be obtained by replacing $$z_{i+1}$$ and $${\it{\Delta}} x$$ with $$z_{I^{'}\!\!,-}$$ and $${\it{\Delta}} x_1$$, respectively; for $$i=I^{'}$$, by replacing $$z_{i}$$, $$z_{i+1}$$ and $${\it{\Delta}} x$$ with $$z_{I^{'}\!\!,+}$$, $$z_{I+1}$$ and $${\it{\Delta}} x_2$$, respectively. Actually, $$(E_{{\rm loc}}^*)_i^l$$ is the numerical error for the ECL (1.12) at the point $$(x_{i+\frac{1}{2}},t_{l+\frac{1}{2}})$$. For simplicity, we hence define the following maximum errors for ECL   (Emax∗)l=maxi{(Eloc∗)il}, (4.7) where $$i$$ runs over all spatial index values. Similarly the discrete local NCL $$(N_{{\rm loc}}^*)_i^l$$ corresponding to (1.14), for $$i\notin \{I,I^{'}\}$$, is given by   (Nloc∗)il=1Δt[Nd(zi+12l+1)−Nd(zi+12l)]+1Δx[Nf(zi+1l+12)−Nf(zil+12)]. And the expressions of $$(N_{{\rm loc}}^*)_{I}^{l}$$ and $$(N_{{\rm loc}}^*)_{I^{'}}^{l}$$ can be similarly obtained as for the local ECL. In response to (4.7), we define the following maximum errors for local NCL   (Nmax∗)l=maxi{(Nloc∗)il}. (4.8) Moreover, the discrete total energy and normalization for MM2 are given by   (EMS)l=Δx∑i∉{I,I′}(S(zi+12l)−12(zi+12l)TKzi+1l−zilΔx)+Δx1(S(zI+12l)−12(zI+12l)TKzI′,−l−zIlΔx1)+Δx2(S(zI′+12l)−12(zI′+12l)TKzI+1l−zI′,+lΔx2) and   (NMS)l=Δx∑i∉{I,I ′}Nd(zi+12l)+Δx1Nd(zI+12l)+Δx2Nd(zI ′+12l), respectively. Thus the relative errors, also referred to as the propagating errors, for the discrete global ECL and NCL are, respectively, defined by   (EMSrel)l=(EMS)l−(EMS)0|(EMS)0|and(NMSrel)l=(NMS)l−(NMS)0|(NMS)0|. (4.9) We remark here that the two non-MS schemes are both discretized for the second-order spatial derivative directly and have no corresponding approximations for the first-order ones, and thus they are not quite suitable to be used in investigating the local conservation laws at the present stage. The corresponding discrete total ones, however, are given by   (EnMS)l=Δx2∑ig[(qil)2+(pil)2]2−Δx2∑i∉{I,I+1}[qi−1l−2qil+qi+1lΔx2qil+pi−1l−2pil+pi+1lΔx2pil]−Δx2[ϕ1(qIl)qIl+ϕ2(qI+1l)qI+1l+ϕ1(pIl)pIl+ϕ2(pI+1l)pI+1l] and   (NnMS)l=Δx∑i[(qil)2+(pil)2] with $$\phi_1$$ and $$\phi_2$$ specified in (4.4). Similarly, the corresponding propagating errors $$(\mathscr{E}_{{\rm nMS}}^{{\rm rel}})^{l}$$ and $$(\mathscr{N}_{{\rm nMS}}^{{\rm rel}})^{l}$$ can be defined as addressed in (4.9). The numerical experiment we will investigate is the case of an attractive nonlinearity, i.e., $$g = -1$$ for the nonlinear Schrödinger equation (1.19). For the sake of better testing our methods, the experiment we choose to implement here possesses the analytic solution (Witthaut et al., 2005)   Ψ(x,t)=ψ(x)exp⁡(i(2γ−1)2t/8), (4.10) where   ψ(x)={λsech(λ(x−x0)),x>0λsech(λ(x+x0)),x<0  with $$\lambda$$ and $$x_0$$ is being subjected to   λ=12−γandtanh(λx0)=γλ. (4.11) For the corresponding stationary nonlinear Schrödinger equations, bound state solutions exist provided $$\gamma<1/4$$ (Witthaut et al., 2005). The above analytic solution concerns the evolution of this bound state solution with the energy $$\mathscr{E}=-(2\gamma-1)^2/8$$. Due to the bound-state-existence condition, in our numerical experimentation here, we adopt a delta-potential barrier with $$\gamma=0.2$$, $$\lambda$$ and $$x_0$$ are then determined by (4.11). Since the exact solution $${\it{\Psi}}(t,x)$$ is exponentially small away from $$x = 0$$ for any fixing $$t$$, the periodic boundary conditions are applied to both the MS and non-MS schemes in numerics. In practical calculations, we take $$x_L = -30$$ and $$x_R = 30$$. We consider two circumstances, which admit, after the first equidistant partition, whether the interface point $$x=0$$ is a uniform grid point or not. If it is, we choose $$N=600$$ and $$N=599$$ for the other. The former means that the spatial stepsize $${\it{\Delta}} x=(x_R-x_L)/N=0.1$$, the latter that $${\it{\Delta}} x=(x_R-x_L)/N\approx 0.1002$$. The temporal stepsize is fixed to be $${\it{\Delta}} t=0.1$$ and the temporal interval $$[0,200]$$. The fixed point iteration method is utilized to solve the nonlinear algebraic systems generated by the numerical schemes, and each iteration will be terminated when the maximum absolute error of the adjacent two iterative values is less than $$10^{-15}$$. The two circumstances are implemented independently. The numerical results plotted here put two legends on them, ‘S1’ for the first circumstance $$N=600$$ and ‘S2’ for the second one $$N=599$$. Figure 1 exhibits the maximum global errors of the numerical wave function obtained by the three schemes. Here we use   max0≤i≤N(max(|qil−q(xi,tl)|,|pil−p(xi,tl)|)) to denote the maximum global error of the wave function at time $$t_l$$, where $$p_i^l$$ is the numerical solution of $$q(x,t)$$ at $$(x_i,t_l)$$ and $$q(x_l,t_i)$$ is the corresponding analytic solution due to (4.10). It follows from the three plots that, all the errors obtained are all in reasonable oscillation and normal linear growth. Considering the fact that the schemes we implemented here are all of second-order accuracy, the errors seem comparatively nice for the temporal step size $${\it{\Delta}} t=0.1$$ we used. The errors for the MS scheme, whether $$x=0$$ is a uniform grid point or not, are both of $$\mathscr{O}(10^{-3})$$ and have almost the same evolutions with time. But those, which are obtained by the other two non-MS schemes, behave very differently for the two circumstances: when $$x=0$$ is a uniform grid point, the errors are of $$\mathscr{O}(10^{-4})$$; and they increase by an order of magnitude, i.e., are of $$\mathscr{O}(10^{-3})$$ for the other. Fig. 1. View largeDownload slide The maximum errors of the wave function obtained by MM2, M2IIM and T2IIM. Fig. 1. View largeDownload slide The maximum errors of the wave function obtained by MM2, M2IIM and T2IIM. Figure 2 displays the propagating errors of the total normalization and energy obtained by MM2. They are given in (4.9). The errors of the normalization, whether $$x=0$$ is a grid point or not, are both in the scale of $$\mathscr{O}(10^{-15})$$, almost the roundoff error of computer. Simultaneously they have similar reasonable linear growths. The errors of the total energy are both in the scale of $$\mathscr{O}(10^{-4})$$ and have stable oscillations something like sinusoidal waves. Fig. 2. View largeDownload slide The propagating errors of the total normalization and energy obtained by MM2, respectively. Fig. 2. View largeDownload slide The propagating errors of the total normalization and energy obtained by MM2, respectively. The propagating errors of the total normalization, which are obtained by M2IIM and T2IIM, are exhibited in Fig. 3. We can see for the two non-MS schemes that the errors for the same conserved quantity look very different whether the interface point is a uniform grid point or not. When it is, the error of the normalization obtained by M2IIM is of $$\mathscr{O}(10^{-14})$$, which is almost as good as that by the MS scheme MM2. In the other circumstance, however, it increases to $$\mathscr{O}(10^{-8})$$, about $$6$$ orders of magnitude lost. For T2IIM, the errors are only of $$\mathscr{O}(10^{-9})$$ and $$\mathscr{O}(10^{-8})$$ for the two circumstances, respectively. Numerical results reveal that the precise normalization conservation is one remarkable advantage of the MS scheme compared with the non-MS schemes. Fig. 3. View largeDownload slide The propagating errors of the total normalization, the left two plots are obtained by M2IIM and the right by T2IIM. Fig. 3. View largeDownload slide The propagating errors of the total normalization, the left two plots are obtained by M2IIM and the right by T2IIM. Figure 4 shows the errors of the total energy obtained by M2IIM and T2IIM. For the first circumstance, the errors obtained by the two non-MS schemes are both of orders $$\mathscr{O}(10^{-9})$$; for the second one, however, they both grow rapidly from $$\mathscr{O}(10^{-9})$$ to $$\mathscr{O}(10^{-5})$$ and then stay in the latter scale. The corresponding errors obtained by MM2, which are displayed in the right plot of Fig. 2, are only of orders $$\mathscr{O}(10^{-4})$$ for both the two circumstances and several orders worse than that given here. We think the reason might be due to different approximations of $$\partial_x v$$ and $$\partial_x w$$ arising in the energy density. Regarding MM2, $$\partial_x v$$ and $$\partial_x w$$ are regarded as the first-order spatial derivatives and thus approximated by the numerical values of $$v$$ and $$w$$, which admit discontinuities at interface points. For the other two non-MS schemes, however, $$\partial_x v$$ and $$\partial_x w$$ are regarded as the second-order ones and approximated by the numerical values of $$p$$ and $$q$$. We remind here that $$\partial_x v$$ and $$\partial_x w$$ are continuous. The discontinuities of $$v$$ and $$w$$ we highlighted in our proposed MS disretizations may be the drawbacks in such approximations. Taking this into account, Bai (2016) proposed MS Runge–Kutta–Nyström methods applied to the delta-potential NLSEs and numerically reveal that the total energy obtained by the implemented MS scheme is preserved much better than by the non-MS one. Fig. 4. View largeDownload slide The propagating errors of the total energy, the left two plots are obtained by M2IIM and the right by T2IIM. Fig. 4. View largeDownload slide The propagating errors of the total energy, the left two plots are obtained by M2IIM and the right by T2IIM. The parity symmetry is one of the three important discrete symmetries, i.e., charge conjugate, parity and time reversal (CPT) symmetries in quantum physics (see, e.g., Cartarius et al., 2012). The exact solution (4.10) we investigate here is of even parity, i.e., the even symmetry of the wave function $${\it{\Psi}}(x,t)={\it{\Psi}}(-x,t)$$. We use   max0≤i≤N(max(|pil−pN−il|,|qil−qN−il|)) to denote the maximum error of the parity conservation at time $$t_l$$. The errors for the two circumstances are presented in Fig. 5. For the MS scheme MM2, the errors are both of orders $$\mathscr{O}(10^{-12})$$; for either M2IIM or T2IIM, the errors are in the scale of $$\mathscr{O}(10^{-10})$$ and $$\mathscr{O}(10^{-11})$$ for the two circumstances, respectively. If the six errors are amplified, we find that they have very similar evolution behaviors. Fig. 5. View largeDownload slide The maximum errors in the parity conservation obtained by MM2, M2IIM and T2IIM. Fig. 5. View largeDownload slide The maximum errors in the parity conservation obtained by MM2, M2IIM and T2IIM. Figure 6 shows the errors of the local conservation laws obtained by the MS scheme MM2. The maximum errors of the local NCL are both of $$\mathscr{O}(10^{-16})$$, almost the roundoff error of the computer, for the two circumstances. The numerical results match our theoretical analysis in Theorem 2.5 very well. The maximum errors of the local ECL, given by (4.7), are in the scale of $$\mathscr{O}(10^{-6})$$ and $$\mathscr{O}(10^{-3})$$ for the two circumstances, respectively. Fig. 6. View largeDownload slide The maximum errors of the local conservation laws by MM2, the left two plots for the local NCL and the right for the local ECL. Fig. 6. View largeDownload slide The maximum errors of the local conservation laws by MM2, the left two plots for the local NCL and the right for the local ECL. 4.2 Numerical experiment B: a double-delta potential For comparison, the numerical experiment we will consider here is the case of the repulsive nonlinearity and the symmetric double-delta potential well, i.e., $$g = 1$$ and $$V_{\rm ext}(x)=\tilde{\gamma}\big[\delta(x-\alpha)+\delta(x+\alpha)\big]$$ ($$\tilde{\gamma}<0,\, \alpha>0$$), respectively, for the NLSE (1.19). The bound state problems for the corresponding stationary linear case have been investigated using the explicit jump IIM and the Peskin’s IB Bai & Wang (2015). The nonlinear case, however, is far distinguished. By virtue of the strategies proposed in Kivshar & Luther-Davies (1998) and Witthaut et al. (2005, 2009), we can construct a solitary wave solution   Ψ~(x,t)=ψ~(x)exp⁡(iλ~2t), (4.12) where   ψ~(x)={2λ~csch(2λ~(x−x~0)),x>α,λ~tan(λ~x),−α≤x≤α,2λ~csch(2λ~(x+x~0)),x<−α.  We stress here that the above is a bound state solution to this nonlinear Schrödinger equation and evolves with the energy $$\tilde{\mathscr{E}}=-\tilde{\lambda}^2$$. One could verify that the above solution is of odd parity, i.e., $$\tilde{{\it{\Psi}}}(-x,t)=-\tilde{{\it{\Psi}}}(x,t)$$, which is in contrast to the even parity proposed in Experiment A. When $$\alpha$$ is fixed, the parameters $$\tilde{\gamma}$$, $$\tilde{\lambda}$$ and $$\tilde{x}_0$$ can be determined by   {2csch(2λ~(α−x~0))=tan⁡(λ~α),2λ~csch(2λ~(α−x~0))coth(2λ~(α−x~0))+λ~sec2⁡(λ~α)=−2γ~tan⁡(λ~α),2λ~tan⁡(λ~α)−2λ~2α+22λ~(coth(2λ~(α−x~0))−1)=1.  (4.13) Here the first equation of (4.13) comes from the continuity of the wave function at the interface points, the second from the jumps and the third from the normalization requirement that   ∫−∞+∞|Ψ~(x,t)|2dx=1. The existence of bound state solutions needs some more requirements, e.g., $$\tilde{x}_0<\alpha$$, $$\tilde{\lambda}\,\alpha<{\pi}/2$$ and so on. In numerics, we adopt $$\alpha=3$$ as investigated for the stationary linear case in Bai & Wang (2015). And the other three parameters we choose are $$\tilde{\lambda}=0.345692582829775$$, $$\tilde{x}_0=1.444829457402318$$ and $$\tilde{\gamma}=-0.775835479368852$$, which solve the transcendental equations (4.13) numerically with the terminating criterion $$10^{-15}$$. Since the exact solution $$\tilde{{\it{\Psi}}}(t,x)$$ is also exponentially small away from $$x = 0$$ for any fixing $$t$$, the periodic boundary conditions are still used in numerical implementation. The spatial domain, the temporal interval, the stop criterion for the nonlinear iteration and the temporal stepsize are all taken to the same as utilized in Experiment A. With regard to the choices of spatial stepsize, we take two circumstances into account. The first concerns that both the two interface points $$x=\pm \alpha$$ are uniform partition grid points and the second is that both are not. For comparisons, we adopt the same two spatial stepsizes as used in Experiment A, which also correspond to the two circumstances in this experiment. For simplicity, we only consider two schemes, the MS one MM2 and the non-MS one T2IIM proposed before. And all the expressions could be obtained similarly and thus would not be presented here. The two circumstances of different spatial stepsizes are also implemented independently and put the two legends ‘D1’ and ‘D2’ on the plots, respectively. Figure 7 presents the maximum errors of numerical solutions obtained by using the two schemes. The two errors are both in the scale of $$\mathscr{O}(10^{-3})$$, just the same magnitude to that for the single-delta case as exhibited in Fig. 1. But the ones by T2IIM given here behave very different from that for the single-delta case: (1) the two errors are of orders $$\mathscr{O}(10^{-1})$$ and $$\mathscr{O}(10^{-2})$$, respectively, which are several orders of magnitude higher than that for the single-delta case as shown in Fig. 1; (2) the first error is one order of magnitude higher than the second one in this double-delta case, but it just reverses in the preceding single-delta case. Fig. 7. View largeDownload slide The maximum errors of numerical solutions, the left is obtained by MM2 and the right by T2IIM. Fig. 7. View largeDownload slide The maximum errors of numerical solutions, the left is obtained by MM2 and the right by T2IIM. Figure 8 exhibits the errors of the total normalization. Both the errors obtained by MM2 here reach almost the roundoff error of computer, just like that in the single-delta case given in the left plot of Fig. 2. The ones by T2IIM are both of $$\mathscr{O}(10^{-7})$$, about two and one orders of magnitude higher than that in the single-delta case, respectively. Fig. 8. View largeDownload slide The propagating errors of the total normalization obtained by MM2 and T2IIM. Fig. 8. View largeDownload slide The propagating errors of the total normalization obtained by MM2 and T2IIM. The propagating errors of the total energy are given in Fig. 9. Those two obtained by MM2 are in the scale of $$\mathscr{O}(10^{-3})$$, about one order of magnitude higher than that in the single-delta case as displayed in Fig. 2. The two errors by T2IIM are of orders $$\mathscr{O}(10^{-7})$$ and $$\mathscr{O}(10^{-5}),$$ respectively, while the corresponding ones are $$\mathscr{O}(10^{-9})$$ and $$\mathscr{O}(10^{-5})$$ in the single-delta case as shown in Fig. 4. Fig. 9. View largeDownload slide The propagating errors of the total energy, the left one obtained by MM2 and the other by T2IIM. Fig. 9. View largeDownload slide The propagating errors of the total energy, the left one obtained by MM2 and the other by T2IIM. Figure 10 gives the errors in the parity conservation. The errors are obtained by using MM2 and are both in the scale of $$\mathscr{O}(10^{-7})$$, almost five orders of magnitude higher than that in the single-delta case as exhibited in Fig. 5. The two errors by T2IIM are of orders of $$\mathscr{O}(10^{-8})$$ and $$\mathscr{O}(10^{-10})$$, respectively, about two and one orders of magnitude higher than that in the single-delta case. Fig. 10. View largeDownload slide The maximum errors in the parity conservation, the left and the middle are obtained by MM2, the right by T2IIM. Fig. 10. View largeDownload slide The maximum errors in the parity conservation, the left and the middle are obtained by MM2, the right by T2IIM. Figure 11 shows the maximum errors of the two local conservation laws. The errors of the local NCL arrive at the roundoff error, just like that in the single-delta case exhibited in Fig. 6. The two errors of the local ECL are both in the scale of $$\mathscr{O}(10^{-3})$$, while the corresponding ones are $$\mathscr{O}(10^{-6})$$ and $$\mathscr{O}(10^{-3})$$, respectively, in the single-delta case. Fig. 11. View largeDownload slide The maximum errors of local conservation laws obtained by MM2, the left and the middle for the local NCL, and the right for the local ECL. Fig. 11. View largeDownload slide The maximum errors of local conservation laws obtained by MM2, the left and the middle for the local NCL, and the right for the local ECL. 5. Conclusions We make an attempt to propose a weak MSHS reformulation for the NLSE with delta potentials. Based on this weak reformulation, some intrinsic local and global conservation laws are newly addressed, a kind of novel concatenating RK discretizations are constructed and the multi-symplecticity of the concatenating algorithms are revealed. Under the MSRK discretizations, some discrete properties are investigated theoretically and, among them, what is important is the precise preservation of the local and global NCLs. Extensive numerical experimentation validates our method. Comparisons with non-MS schemes reveal some priorities of our MS formulations, and the first is of course the roundoff-error preservation of the NCLs. Furthermore, the longitudinal comparisons made between two circumstances show that for either the single- or double-delta case the errors by our MS schemes vary much less than that by the non-MS ones whether the interface points are uniform grid points or not; the transverse comparisons made between two cases reveal that, the errors by our MS schemes change much less than that by non-MS ones when the number of delta potentials change from one to two. It might be concluded that the MS schemes have more stable dynamical numerical behaviors than the non-MS ones. In physics, there exist a large number of important PDEs which admit various jumps or discontinuities, e.g., the nonlinear Dirac equations, the Klein–Gordon equations and many other coupled nonlinear equations with several delta potentials, some hydrodynamic equations with discontinuous coefficients or singular sources (Landau & Lifshitz, 1987; Haus & Wong, 1996; Kivshar & Luther-Davies, 1998; Bridges & Reich, 2006; Cartarius et al., 2012 and references therein) and so on. They may be explored in a similar way as proposed in this paper. this context, although further research is needed. Acknowledgements We would like to thank the anonymous referees for their valuable comments and suggestions. Funding National Natural Science Foundation of China (Nos.11001125 and 91130003) and A Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions. References Arnold, V. I. ( 1999) Mathematical Methods of Classical Mechanics , 2nd edn. New York: Springer. Bai, J. ( 2016) Multi-symplectic Runge-Kutta-Nyström methods for nonsmooth nonlinear Schrödinger equations. J. Math. Anal. Appl. , 444, 721– 736. Google Scholar CrossRef Search ADS   Bai, J. & Wang, L. ( 2015) EJIIM for the stationary Schrödinger equations with delta potential wells. Appl. Math. Comput. , 254, 113– 124. Berezin, F. A. & Shubin, M. A. ( 1991) The Schrödinger Equation . Dordrecht: Kluwer Academic Publishers. Google Scholar CrossRef Search ADS   Bridges, T. J. ( 1997) Multi-symplectic structures and wave propagation. Math. Proc. Camb. Phil. Soc. , 121, 147– 190. Google Scholar CrossRef Search ADS   Bridges, T. J. & Reich, S. ( 2006) Numerical methods for Hamiltonian PDEs. J. Phys. A Math. Gen. , 39, 5287– 5320. Google Scholar CrossRef Search ADS   Carretero-González, R., Frantzeskakis, D. J. & Kevrekidis, P. G. ( 2008) Nonlinear waves in Bose-Einstein condensates: physical relevance and mathematical techniques. Nonlinearity , 21, R139– R202. Google Scholar CrossRef Search ADS   Cartarius, H., Haag, D., Dast, D. & Wunner, G. ( 2012) Nonlinear Schrödinger equation for a $$\mathscr{PT}$$-symmetric delta-function double well. J. Phys. A Math. Theor. , 45, 444008 ( 15 pp). Google Scholar CrossRef Search ADS   Dalfovo, F., Giorgini, S., Pitaevskii, L. P. & Stringari, S. ( 1999) Theory of Bose-Einstein condensation in trapped gases. Rev. Mod. Phys. , 71, 463– 512. Google Scholar CrossRef Search ADS   Hairer, E., Lubich, C. & Wanner, G. ( 2006) Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations , 2nd edn. Berlin: Springer. Haus, H. A. & Wong, W. S. ( 1996) Solitons in optical communications. Rev. Mod. Phys. , 68, 423– 444. Google Scholar CrossRef Search ADS   Hong, J. & Li, C. ( 2006) Multi-symplectic Runge-Kutta methods for nonlinear Dirac equations. J. Comput. Phys. , 211, 448– 472. Google Scholar CrossRef Search ADS   Hong, J., Liu, H. & Sun, G. ( 2005) The multi-symplecticity of partitioned Runge-Kutta methods for Hamiltonian PDEs. Math. Comput. , 75, 167– 181. Google Scholar CrossRef Search ADS   Hong, J., Liu, X. & Li, C. ( 2007) Multi-symplectic Runge-Kutta-Nyström methods for nonlinear Schrödinger equations with variable coefficients. J. Comput. Phys. , 226, 1968– 1984. Google Scholar CrossRef Search ADS   Islas, A. L., Karpeev, D. A. & Schober, C. M. ( 2001) Geometric integrators for the nonlinear Schrödinger equation. J. Comput. Phys. , 173, 116– 148. Google Scholar CrossRef Search ADS   Kivshar, Y. S. & Luther-Davies, B. ( 1998) Dark optical solitons: physics and applications. Phys. Rep. , 298, 81– 197. Google Scholar CrossRef Search ADS   Kong, L., Hong, J., Ji, L. & Zhu, P. ( 2015) Compact and efficient conservative schemes for coupled nonlinear Schrödinger equations. Numer. Meth. PDE , 31, 1814– 1843. Google Scholar CrossRef Search ADS   Landau, L. D. & Lifshitz, E. M. ( 2008) Fluid Mechanics , 2nd edn. Beijing: Beijing World Publishing Corporation. McLachlan, R. I., Ryland, B. N. & Sun, Y. ( 2014) High order multisymplectic Runge-Kutta methods. SIAM J. Sci. Comput. , 36, A2199– A2226. Google Scholar CrossRef Search ADS   Olver, P. J. ( 1993) Applications of Lie Groups to Differential Equations , 2nd edn. New-York: Springer. Google Scholar CrossRef Search ADS   Pethick, C. J. & Smith, H. ( 2002) Bose-Einstein Condensation in Dilute Gases . New-York: Cambridge University Press. Rapedius, K. & Korsch, H. J. ( 2009) Resonance solutions of the nonlinear Schrödinger equation in an open double-well potential. J. Phys. B , 42, 044005 ( 12 pp). Google Scholar CrossRef Search ADS   Reich, S. ( 2000) Multi-symplectic Runge-Kutta collocation methods for Hamiltonian wave equation. J. Comput. Phys. , 157, 473– 499. Google Scholar CrossRef Search ADS   Ryland, B. N. & McLachlan, R. I. ( 2008) On multi-symplecticity of partitioned Runge-Kutta methods. SIAM J. Sci. Comput. , 30, 1318– 1340. Google Scholar CrossRef Search ADS   Sulem, C. & Sulem, P.-L. ( 1999) The Nonlinear Schrödinger Equation: Self-Focusing and Wave Collapse . New York: Springer. Witthaut, D., Mossmann, S. & Korsch, H. J. ( 2005) Bound and resonance states of the nonlinear Schrödinger equation in simple model systems. J. Phys. A Math. Gen. , 38, 1777– 1792. Google Scholar CrossRef Search ADS   Witthaut, D., Rapedius, K. & Korsch, H. J. ( 2009) The nonlinear Schrodinger equation for the delta-comb potential: quasi-classical chaos and bifurcations of periodic stationary solutions. J. Nonlinear Math. Phys. , 16, 207– 233. Google Scholar CrossRef Search ADS   © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

Journal

IMA Journal of Numerical Analysis – Oxford University Press

Published: Jan 1, 2018

There are no references for this article.