Cite this Article
Fu Hao, Zhou Wei-En, Qian Xu, Song Song-He, Zhang Li-Ying. Conformal structure-preserving method for damped nonlinear Schrödinger equation. Chinese Physics B, 2016, 25(11): 110201  
Permissions
Conformal structure-preserving method for damped nonlinear Schrödinger equation
Fu Hao1, †, , Zhou Wei-En1, Qian Xu1, Song Song-He1, Zhang Li-Ying2
College of Science, National University of Defense Technology, Changsha 410073, China
School of Mathematical Science, China University of Mining and Technology, Beijing 100083, China

 

† Corresponding author. E-mail: fuhaosnsn@126.com

Project supported by the National Natural Science Foundation of China (Grant Nos. 11571366, 11501570, and 11601514) and the Open Foundation of State Key Laboratory of High Performance Computing of China (Grant No. JC15-02-02).

Abstract
Abstract

In this paper, we propose a conformal momentum-preserving method to solve a damped nonlinear Schrödinger (DNLS) equation. Based on its damped multi-symplectic formulation, the DNLS system can be split into a Hamiltonian part and a dissipative part. For the Hamiltonian part, the average vector field (AVF) method and implicit midpoint method are employed in spatial and temporal discretizations, respectively. For the dissipative part, we can solve it exactly. The proposed method conserves the conformal momentum conservation law in any local time–space region. With periodic boundary conditions, this method also preserves the total conformal momentum and the dissipation rate of momentum exactly. Numerical experiments are presented to demonstrate the conservative properties of the proposed method.

PACS: 02.60.Lj;02.60.–x;02.60.Cb;02.70.Bf
Keyword:conformal conservation law;splitting method;conformal momentum-preserving method;damped multi-symplectic formulation
1. Introduction

We consider the following damped nonlinear Schrödinger (DNLS) equation:

where γ > 0 is a constant parameter and a > 0 is a damping coefficient. If a = 0, equation (1) is the general nonlinear Schrödinger equation (NLS), and it is one of the most important models of mathematical physics, with application to different fields such as plasma physics, nonlinear optics, water waves, bimolecular dynamics, and many other fields. It plays an essential role in mathematical and physical contexts. The general NLS equation can be cast into a multi-symplectic form; thus many geometric integrators, which can preserve the physical properties of the original problem under appropriate discretizations,[1–8] are used to solve it (see Refs. [9]–[12] and references therein). Bridges and Reich[13] proposed a multi-symplectic Preissman scheme by using an implicit midpoint rule, then introduced the idea of a multi-symplectic Fourier transform for multi-symplectic PDEs with periodic boundary conditions in Ref. [14]. Next, Chen and Qin proposed a multi-symplectic Fourier pseudospectral method on real space.[15] In Ref. [16], it has been shown that concatenating symplectic Nyström methods in the spatial direction and symplectic Runge–Kutta methods in the temporal direction leads to multi-symplectic integrators for the NLS equation by Hong et al. Except the multi-symplectic conservation law, multi-symplectic Hamiltonian partial differential equations (PDE) also have the energy and momentum conservation laws. It is well known that conservation law plays an important role in conservative PDEs, especially in the theory of solitons. In Ref. [17], Gong et al. gave several systematic methods for discretizing general multi-symplectic formulations of Hamiltonian PDEs, including a local energy-preserving algorithm, a class of global energy-preserving methods and a local momentum-preserving algorithm by using the AVF method and implicit midpoint method.

For the damped equation (1), it describes several resonant phenomena in nonlinear media, in particular, the nonlinear Faraday resonance in a vertically oscillating water trough[18] and the effect of phase-sensitive amplifiers on solitons in optical fibers.[19] In the context of symmetric Hamiltonian systems with linear damping, which are known as conformal Hamiltonian systems,[20] McLachlan and Perlmutter developed a reduction conformal theory and showed that standard symplectic methods do not generally preserve the conformal symplectic structure, but conformal symplectic methods do. In Ref. [21], McLachlan and Quispel have constructed numerical methods that preserve conformal properties, such as symplecticity and volume preservation for ODEs. Moore generalized these results to multi-symplectic PDEs with added dissipation, and proposed a conformal multi-symplectic scheme for the force-damped semi-linear wave equation in Ref. [22]. Subsequently, Moore et al.[23] derived conformal energy, momentum, and other quantities that arise from linear symmetry conservation laws and developed conformal Preissman box scheme and conformal discrete gradient methods. These methods are proven to preserve the conservation laws exactly. Comparing with standard structure-preserving algorithms, conformal structure-preserving algorithms preserve the dissipation rate exactly. Our aim in this article is to develop a numerical method that not only preserves a conformal property, but also conserves the dissipation rate. In order to achieve this goal, a conformal momentum-preserving method is constructed by the standard conservative method and splitting technique. In other words, after using the splitting method in the temporal direction, we discrete the split system using the implicit midpoint scheme in the temporal direction and AVF method in the spatial direction respectively.

The rest of the paper is organized as follows. In Section 2, the damped multi-symplectic formulation of the DNLS equation and some of its conformal conservation laws are presented. In Section 3, some necessary operators together with their properties are given. Then, we derive a conformal momentum-preserving method for the DNLS equation, and also prove some conservative properties of the method. Numerical experiments for the solitary wave solution are presented in Section 4, which show the effectiveness and advantages of this method. Finally, we make some conclusions.

2. Damped multi-symplectic formulation and conformal conservation laws

The form of the damped multi-symplectic formulation[23] is

where z ∈ ℝd, M, and K are skew-symmetric matrices, S : ℝd → ℝ is a smooth function, ∇ is a gradient operator, and a and b are non-negative real numbers. Many damped PDEs can be rewritten in this form, such as a semi-linear wave equation with a simple linear damping term, a damped Klein–Gordon equation, etc. The corresponding conformal multi-symplectic conservation law is

where and . If a = 0, the corresponding conformal energy conservation law is

where and . If b = 0, the corresponding conformal momentum conservation law is

where and . Notice that integrating Eq. (5) in space, with the periodic boundary condition, we obtain the total conformal momentum as

where . We say that the total conformal momentum is preserved if a numerical method satisfies

where is the temporal step size, and j is the spatial index.

Consider the DNLS equation (1). Let ψ = p + iq, equation (1) can be written as

We introduce a pair of conjugate momenta v = px, w = qx and obtain

where

and Hamiltonian . Notice that, S does not depend on damping coefficient a in this case.

Proposition 1 The DNLS (1) admits the conformal norm conservation law (CNCL)

and the conformal momentum conservation law (CMCL)

Proof Multiplying Eq. (8) by 2p and Eq. (9) by −2q gives

Adding Eq. (13) to Eq. (14) and noting

we obtain CNCL (11).

Multiplying Eq. (8) by qx and Eq. (9) by px obtains

Adding Eq. (15) to Eq. (16) and noting

so we obtain CMCL (12).

Integrating Eqs. (11) and (12) with respect to x with the periodic boundary condition, we obtain the total conformal norm and total conformal momentum

3. Conformal momentum-preserving method for DNLS equation
3.1. Operator definitions and properties

In order to derive the algorithms conveniently, we give some notations and operator definitions. We introduce a uniform grid (xj, tn) ∈ ℝ2 with mesh-length Δt in the temporal direction and mesh-length Δx in the spatial direction, and denote the value of the function z(x, t) at the mesh point (xj, tn) by , where xj = −L + jΔx, j = 0,…,J, and tn = nΔt, n = 0,…,N. Define the finite difference operators as

and averaging operators

then we have the following properties:[23]

Commutative law:

Generalized discrete Leibnitz rule:

For some special cases, we have

Similarly, we can obtain a series of analogous Leibnitz rules in the temporal direction. These discrete Leibnitz rules are essential for proving a method that preserves any of the conservation laws mentioned in this article.

As we know, the standard geometric integrators generally cannot preserve the conformal conservation law. Mclachlan and Quispel[21] have shown that splitting methods may be used to numerically preserve the conformal symplecticity for Hamiltonian ordinary differential equations. In this article, we use a splitting technique to construct the numerical method that preserves conformal conservation laws exactly from the standard conservative method. To explain the process, we consider the dynamical system

where f1(z) is the Hamiltonian vector field and f2(z) = Dz is a dissipative term. This dynamical system can be split into the Hamiltonian part zt = f1(z) and the dissipative part zt = f2(z). For the Hamiltonian part, we apply the standard geometric integrator to obtain zn+1 = ΨΔt (zn), and for the dissipative part, we solve it exactly, i.e., zn+1 = ΦΔt(zn) = expDΔt (zn). Thus, the system (19) may be solved by decomposing the flow maps ΦΔt ○ ΨΔt, which motivate the definitions of non-standard difference operators

along with the averaging operators

Obviously, if we take a = 0, then and . We check that the operators and also satisfy the following properties:[23]

Commutative law:

Discrete Leibnitz rules:

3.2. Full-discrete scheme of DNLS equation

In this subsection, we would like to construct a conformal momentum-preserving (CMP) method and study the long-time behavior of such a method. This conformal method is derived from the standard conservative method and splitting technique. This idea is to use the implicit midpoint scheme and the AVF method. At first, we use a splitting technique in the temporal direction, then discretize time derivatives by using the implicit midpoint scheme and space derivatives by using the AVF method, finally obtaining a full-discrete system

Theorem 1 The full-discrete system (20) exactly conserves the following N-discrete conformal momentum conservation law:

for j = 0,1,…,N − 1, where

Proof Using the communicative law and the discrete Leibnitz rule shows

Taking the inner product of Eq. (20) by and considering that

Thus, we have

Corollary 1 With the periodic or homogeneous boundary conditions, the conformal momentum-preserving method preserves the total conformal momentum exactly, namely,

where the total conformal momentum . That is to say, the conformal momentum-preserving method preserves the dissipation rate of the total momentum exactly.

In order to prove Corollary 1, we only need to sum the discrete conformal conservation law (21) over all space index j and use the periodic or homogeneous boundary conditions. Because of the exponential dissipation of the total conformal momentum (18), we can use

to measure the dissipation rate of the momentum. If r = 0, we say the numerical method preserves the dissipation rate of the momentum. Thus, the method we propose conserves the dissipation rate of the momentum. However, the standard structure-preserving algorithms do not, in general, preserve the conformal conservation law and dissipation rate; this is explained and demonstrated in Ref. [23].

Applying the conformal momentum-preserving method for the damped multi-symplectic formulation of the DNLS equation (10), we obtain the scheme as follows:

Eliminating the auxiliary variables, we obtain the following scheme:

where

According to Theorem 1, the scheme (25) or (26) preserves the discrete conformal momentum conservation law

for j = 0,…,N − 1, where

The numerically induced residual RM of the conformal momentum conservation law is given by

The corresponding discrete total conformal momentum is

and the discrete total conformal norm is

4. Numerical experiments

In this section, we conduct several typical experiments to illustrate the performance of the proposed method. All the solutions are computed with a time step of Δt = 0.01 and periodic boundary conditions in [−64,64]. We choose grid number N = 513 and the time interval is taken as [0,60] for all examples. To test the conservative properties, we conduct long time simulation for one soliton propagation and double solitons collision propagation in Examples 1 and 2, respectively. For the purpose of numerical comparisons, we apply a standard local momentum-preserving (LMP) algorithm,[17] implicit midpoint scheme in time and the AVF method in space to the DNLS equation. In addition, we solve the implicit nonlinear equations by using the fixed-point iteration method with tolerance ɛ = 10−13.

Example 1 One-soliton solution. In optics, a soliton wave arises from a balance between nonlinear and dispersive effects. Many physical systems can be modeled quite successfully using equations that admit soliton solutions. Indeed, solitons and solitary waves have been observed in numerous situations and often dominate long-time behavior. For the DNLS equation (1), if γ ≪ 1, the soliton will be affected by the linear effect (diffraction) much earlier than the nonlinear effect, it will just diffract without any nonlinear behavior. If γ ≫ 1, the nonlinear effect will be more evident than the linear effect, and because of self phase modulation, the soliton will change its shape periodically. If γ ≈ 1, the two effects balance each other. In this example, we choose γ = 1, that is to say, the shape of soliton does not change except the damping property during propagation. Consider the cubic DNLS equation with initial condition

where p = 20 is the initial phase. We choose different damping coefficients a = 0.01 and a = 0.03 to observe the evolution of the periodic solution, respectively.

Figures 1 and 2 show the results by using the CMP method with a = 0.01 and a = 0.03, respectively. As can be seen from Figs. 1(a) and 2(a), the propagation of solitary wave over the time interval [0,60] is traveling from left and right as required, presenting the good preservation of the phase space structure. Figures 1(b) and 2(b) show the shape of propagation of solitary at t = 0 and t = 60. We can easily observe that the case of a = 0.03 is damping faster than that of a = 0.01. The evolution of discrete conformal norm and momentum are exhibited in Figs. 1(d), 1(e) and 2(d), 2(e), respectively. Figures 1(c) and 2(c) display the residuals of CMCL, and the error of CMCL is O(10−12). Figures (f) and (f) show the errors of the dissipation rate of momentum is about O(10−12). The CMCL and dissipation rate are conserved up to the accuracy of the CMP method, which support the results of Theorem 1 and Corollary 1.

Fig. 1. The one-soliton solution with a = 0.01. (a) Propagation of solitary wave |ψ|; (b) the shape of solution at t = 0 and t = 60; (c) the residual of ; (d) the evolution of total conformal momentum GM; (e) the evolution of total conformal norm GN; (f) error of the dissipation rate of momentum r.
Fig. 2. The one-soliton solution with a = 0.03. (a) Propagation of solitary wave |ψ|; (b) the shape of solution at t = 0 and t = 60; (c) the residual of ; (d) the evolution of total conformal momentum GM; (e) the evolution of total conformal norm GN; (f) error of the dissipation rate of momentum r.

Example 2 Two-soliton solution. As the fibre technology advances, the interest in optical solitons grows rapidly. Various soliton collision scenarios such as transmission, reflection, and creation of a new soliton have been reported. If a system is integrable, solitary waves collide elastically; that is, they preserve their shape after collision. However, if the system is non-integrable, the collision may be highly non-trivial and inelastic; that is, the shapes of the solitons change after collision. Additionally, the soliton interactions can lead to large and rapidly decaying oscillating radiative tails. In this example, we choose γ = 1, that is to say, the DNLS system (1) is integrable. Considering a multiple solitary wave interaction, the initial condition is given by

where the initial phase is p = 20. The initial condition represents two waves separated by a distance equal to 40 units: the faster one initially located at x = −20 and the slower one at x = 20. In this case, the initial phases of two solitary waves are symmetric with equal amplitude and different velocities, which ensures that the initial momentum is not equal to zero. We also choose different damping coefficients a = 0.01 and a = 0.03 to observe the evolutions of the periodic solution, respectively.

Figures 3 and 4 provide the results by using the CMP method with a = 0.01 and a = 0.03, respectively. As can be seen from Figs. 3(a) and 4(a), the interaction of two solitons over the time interval [0,60] is traveling as required, presenting the good preservation of the phase space structure. The two solitons leave the interaction region without any disturbances in their identities. This phenomenon shows that the collision is elastic. Hence, the DNLS system (1) is integrable. Figures 3(b) and 4(b) display the shape of the numerical solution at t = 0, t = 23 (collision time), and t = 60. We can easily observe that the case a = 0.03 is damping faster than the case a = 0.01. The changes of discrete conformal norm and momentum are exhibited in Figs.3(d), 3(e) and Figs. 4(d), 4(e), respectively. We also plot the residuals of CMCL as a function of the spatial grid location and the time step in Figs. 3(c) and 4(c). Figures 3(f) and 4(f) show that the dissipation rate error of momentum is O(10−12). These results show that the CMP method has excellent structure-preserving ability.

Fig. 3. The two-soliton solution with a = 0.01. (a) Propagation of solitary wave |ψ|; (b) the shape of solution at t = 0, t = 23, and t = 60; (c) the residual of ; (d) the evolution of total conformal momentum GM; (e) the evolution of total conformal norm GN; (f) error of the dissipation rate of momentum r.
Fig. 4. The two-soliton solution with a = 0.03. (a) Propagation of solitary wave |ψ|; (b) the shape of solution at t = 0, t = 23, and t = 60; (c) the residual of ; (d) the evolution of total conformal momentum GM; (e) the evolution of total conformal norm GN; (f) error of the dissipation rate of momentum r.

In Figs. 5 and 6, we make a comparison of the standard LMP method and the CMP method for DNLS equation with a = 0.01 and a = 0.03, respectively. Figures 5(a), 5(b) and 6(a), 6(b) show the comparison of the numerical solution at t = 60 between the CMP method and the standard LMP method, showing that both of the methods preserve the phase space structure well. As can be seen in Figs. 5(c), 5(d) and 6(c), 6(d), the conformal method is significantly more accurate in capturing the damped momentum than the standard method, i.e., the error in the total conformal momentum with the conformal method is O(10−12), while the standard method is O(10−8).

Fig. 5. Comparison of the standard LMP method and the CMP method for the case a = 0.01. (a) Comparison of the one-soliton solution at t = 60; (b) comparison of the two-soliton solution at t = 60; (c) the dissipation rate of momentum for the one-soliton case; (d) the dissipation rate of momentum for the two-soliton case.
Fig. 6. Comparison of the standard LMP method and the CMP method for the case a = 0.03. (a) Comparison of the one-soliton solution at t = 60; (b) comparison of the two-soliton solution at t = 60; (c) the dissipation rate of momentum for the one-soliton case; (d) the dissipation rate of momentum for the two-soliton case.
5. Conclusion

In this paper, we have derived conformal norm conservation law and conformal momentum conservation law, describing the dissipation of norm and momentum, for the DNLS equation. All these conservation laws are local. That is to say, they are independent of the boundary conditions. Based on its multi-symplectic formulation, we construct a conformal momentum-preserving method; this is proven to preserve the conformal momentum conservation law exactly. With periodic boundary conditions, this method also preserves the total conformal momentum and the dissipation rate of momentum exactly. The method proposed is applicable to other PDEs which can be rewritten as damped multi-symplectic formulation. Numerical results show that the proposed method is effective for the DNLS equation. Compared with the standard local momentum-preserving method, we find that the standard method applied to the DNLS equation destroys the rate of dissipation, while the conformal method conserves it exactly. To sum up, the conformal method is more advantageous for long time simulation than the standard method for the DNLS equation.

Reference
1Zhang HSong S HChen X DZhou W E2014Chin. Phys. B23070208
2Qian XSong S HGao ELi W B2012Chin. Phys. B21070206
3Sun J QGu X YMa Z Q2004Phys. D196311328
4Chen Y MSong S HZhu H J 2011 J. Comput. Appl. Math. 236 1354
5Cai J X 2010 Appl. Math. Comput. 216 2417
6Liao C CCui J CLiang J ZDing X H2016Chin. Phys. B25010205
7Cai J XWang J LWang Y S2015Chin. Phys. B24050205
8Cai W JWang Y SSong Y Z2014Chin. Phys. Lett.31040201
9Chen J BQin M ZTang Y F 2002 Comput. Math. Appl. 43 1095
10Fu F FKong L HWang L2009Adv. Appl. Math. Mech.1699
11Islas A LKarpeev D ASchober C M 2001 J. Comput. Phys. 173 116
12Moore B EReich S2003Future Gener. Comput. Syst.19395
13Bridges T JReich S 2001 Phys. Lett. 284 184
14Bridges T JReich S 2001 Phys. 152 491
15Chen J BQin M Z 2001 Electron. Trans. Numer. Anal. 12 193
16Hong J LLiu X YLi C 2007 J. Comput. Phys. 226 1968
17Gong Y ZCai J XWang Y S 2014 J. Comput. Phys. 279 80
18Elphick CMeron E1989Phys. Rev. A403226
19Deutsch I HAbram I1994J. Opt. Soc. Am. B112303
20Perlmutter MMcLachlan R2001J. Geom. Phys.39276
21McLachlan R IQuispel G R W2002Acta Numer.11699
22Moore B E 2009 Math. Comput. Simulat. 80 20
23Moore B ENoreña LSchober C M 2013 J. Comput. Phys. 232 80