Multi-symplectic variational integrators for nonlinear Schrödinger equations with variable coefficients
Liao Cui-Cui 1, †, , Cui Jin-Chao 1 , Liang Jiu-Zhen 2 , Ding Xiao-Hua 3
Department of Information and Computing Science, College of Science, Jiangnan University, Wuxi 214122, China
College of Internet of Things, Jiangnan University, Wuxi 214122, China
Department of Mathematics, Harbin Institute of Technology. Harbin 150001, China

 

Project supported by the National Natural Science Foundation of China (Grant No. 11401259) and the Fundamental Research Funds for the Central Universities, China (Grant No. JUSRR11407).

Abstract
Abstract

In this paper, we propose a variational integrator for nonlinear Schrödinger equations with variable coefficients. It is shown that our variational integrator is naturally multi-symplectic. The discrete multi-symplectic structure of the integrator is presented by a multi-symplectic form formula that can be derived from the discrete Lagrangian boundary function. As two examples of nonlinear Schrödinger equations with variable coefficients, cubic nonlinear Schrödinger equations and Gross–Pitaevskii equations are extensively studied by the proposed integrator. Our numerical simulations demonstrate that the integrator is capable of preserving the mass, momentum, and energy conservation during time evolutions. Convergence tests are presented to verify that our integrator has second-order accuracy both in time and space.

1. Introduction

Nonlinear Schrödinger equations have wide applications in many areas such as quantum mechanics, nonlinear optics, [ 1 ] and plasma physics, [ 2 ] etc. Extensive efforts have been devoted to studying the equation theoretically [ 3 ] and numerically, [ 4 14 ] due to its broad and important applications. Many cases of practical interest, such as the dispersion-managed optical fibers and soliton lasers, certain inhomogeneous optical fibers, arterial mechanics, and laser-atom interactions, [ 15 17 ] lead to nonlinear Schrödinger equations with variable coefficients.

In this paper, we consider nonlinear Schrödinger equations with variable coefficients of the form

where , variable coefficients α ( t ), ν ( x ), β ( t ) are bounded real functions, α ( t ) is related to the second-order dispersion coefficient, and ψ ( x , t ) is a complex-valued wave function. [ 18 , 19 ]

Most of the real physical processes with negligible dissipation can be cast in suitable Hamiltonian formulation in phase space with multi-symplectic geometric structure. [ 20 24 ] Such an intrinsic geometry structure is preserved by the solution flow of Eq. ( 1 ). By the Legendre transforms, [ 25 ] p t = ψ t and p x = ψ x , equation ( 1 ) can be reformulated as multi-symplectic Hamiltonian equations [ 20 , 21 ]

where Z = ( ψ , p t , p x ) T and S ( Z , x , t ) = H ( ψ , p t , p x , x , t ) with the momenta of time p t and the momenta of space p x . Here, M and K ( t ) are antisymmetric matrices on ℝ 3 and H ( ψ , p t , p x , x , t ) is the Hamiltonian function. [ 26 ] An intrinsic conservative property of Hamiltonian systems, multi-symplecticity, is presented by

for the differential forms ω ( U , V ) = 〈 MU , V 〉 and κ ( U , V ) = 〈 K ( t ) U , V 〉, where 〈·, ·〉 is the inner product. Here, U ( x , t ) and V ( x , t ) are the solutions of the variational equation of Eq. ( 2 ),

where D zz S is the second-order derivative of S with respect to Z , D zz S is a symmetric matrix function.

Alternatively, equation ( 1 ) can be reformulated as Euler–Lagrangian equations

where L ( ψ , ψ t , ψ x ) is a Lagrangian function. Exact solution flows of nonlinear Schrödinger equations ( 1 ) preserve the multi-symplectic geometric structures, i.e., multi-symplectic form formulas ( 5 ) in the following proposition.

Proposition 1 [ 23 , 24 ] If ψ is a solution of Euler–Lagrangian equation, and V and W are first variations of ψ , then for any subset U of the space of independent variables, the following multi-symplectic form formula holds:

where * is a pullback operator, j 1 ψ is the first jet of ψ , and Ω L is the Lagrangian multi-symplectic form.

In coordinates, j 1 ψ is given by x μ ↦ ( x μ , ψ A ( x μ ), μ ψ A ( x μ )). More details on the notations in Eq. ( 5 ) are referred to Refs. [ 23 ] and [ 24 ] and references therein.

It follows from the Legendre transforms [ 25 ] that the Hamiltonian equations ( 2 ) and the Lagrangian equations ( 4 ) are equivalent, and the multi-symplectic conservation law ( 3 ) is actually the multi-symplectic form formula ( 5 ). The readers are referred to Refs. [ 23 ] and [ 24 ] for detailed differential geometric connection proofs.

In addition to the geometric multi-symplectic structure, the nonlinear Schrödinger equation ( 1 ) admits the following global conservation laws.

Proposition 2 [ 18 , 19 ] If the wave function ψ is the solution of Eq. ( 1 ), then this wave function satisfies the following relations:

Mass conservation

Momentum conservation

where ℜ and 𝔍 stand for the real part and the imaginary part, respectively.

Energy conservation

If α ( t ) and β ( t ) are independent of t (i.e., α ( t ) ≡ α and β ( t ) ≡ β ) then

where α and β are constants.

The nonlinear Schrödinger equation has been extensively studied by various numerical methods, such as finite element methods, [ 27 ] finite difference methods, [ 28 ] spectral methods, [ 29 , 30 ] splitting methods, [ 6 , 7 , 18 , 19 , 31 33 ] etc. Among these numerical methods of different categories, the multi-symplectic method has attracted special attention for its better numerical stability for long-time computations and perfect performance in preserving the intrinsic properties and conservation laws of nonlinear Schrödinger equations.

As an intrinsic geometric structure of the original system, the multi-symplectic structure, i.e., multi-symplectic conservation law ( 3 ) or multi-symplectic form formula ( 5 ) in Proposition 1, is perfectly preserved by multi-symplectic methods. Thanks to such an advantage, the methods are known for their better numerical stabilities on long-time computations. Furthermore, multi-symplectic methods preserve the quantities of conservation laws, e.g., the mass conservation, momentum conservation, and energy conservation in Proposition 2.

Many works have been devoted to the derivation of multi-symplectic numerical methods for nonlinear Schrödinger equations. There are mainly two ways to construct multi-symplectic schemes. One is based on the Hamiltonian point of view. Hong et al. [ 26 , 34 ] proposed a numerical scheme for nonlinear Schrödinger equations with variable coefficients by means of Preissman integrators. For this scheme, they derived a discrete multi-symplectic structure, i.e., discrete multi-symplectic conservation law. [ 21 , 22 ] In addition, they developed the discrete normal conservation law and a global energy transit formula in temporal direction. Chen et al. studied multi-symplectic methods for Hamiltonian equations of nonlinear Schrödinger equations. [ 35 , 36 ] After applying a numerical discretization to Hamiltonian equations, one needs to re-derive multi-symplectic conservation law since it is unclear what is geometrically conserved by this discretization.

Alternatively, the multi-symplectic numerical schemes, based on the Lagrangian viewpoint and variational principle, [ 37 , 38 ] lead in a natural way to multi-symplectic integrators and the discrete multi-symplectic structures at the same time. From the Lagrangian viewpoint, Chen et al. [ 39 , 40 ] have elaborately studied the variational multi-symplectic integrators for the nonlinear Schrödinger equations. By the discrete variational principle with the discrete Lagrangian function, the discrete variational integrator is derived, and the corresponding multi-symplectic form formula, due to Marsden, [ 23 , 24 ] is also obtained from the variational principle.

In our previous work, [ 41 ] we took a new approach to provide a criterion of multi-symplecticity which can be applied to any set of PDEs where a Dirichlet-to-Neumann map can be defined. We defined an analogue of Jacobi’s solution to the Hamilton–Jacobi equation for field theories, and showed that, by taking variational derivatives of this functional, an isotropic submanifold of the space of Cauchy data, described by the so-called multisymplectic form formula, can be obtained as well. A similar framework of generating functions for discrete field theories was presented in the end of the paper.

In this work, we follow this Lagrangian viewpoint to study the multi-symplectic methods for the nonlinear Schrödinger equations with variable coefficients ( 1 ). One of main aims of this work is to construct variational integrators for Eq. ( 1 ), and derive discrete multi-symplectic form formulas by the definition of discrete Lagrangian boundary function. [ 41 ] Equations ( 1 ) are a general form of the following equations.

Cubic nonlinear Schrödinger equation (cubic NLSE)

Theoretical investigations of this equation are referred to Ref. [ 42 ] and references therein. For this case, we assume that the solution ψ exists globally and satisfies

Gross–Pitaevskii equation (GPE)

The existence of stationary solutions to this equation was studied in Ref. [ 43 ] and the references therein. For this case, we assume that the solution ψ satisfies

The numerical conservation behaviors of the proposed multi-symplectic methods for cubic NLSE ( 9 ) and GPE ( 10 ) are presented in Section 3.

2. Variational integrator for nonlinear Schrödinger equations with variables coefficients

We consider a nonlinear Schrödinger equation with variable coefficients (see Eq. ( 11 ))

This equation can be reformulated as a Euler–Lagrange equation

with a Lagrangian function

where is the conjugation of ψ .

Assume that we are given the regular quadrangular mesh in the base space, with mesh lengths Δ x and Δ t . The nodes in this mesh are denoted by ( j , k ) ∈ ℤ × ℤ, corresponding to the points ( x j , t k ) ≔ ( j Δ x , k Δ t ) in ℝ 2 . We denote by the value of the field ψ at the node ( j , k ) by . If we consider the square discretization, we denote a square at ( j , k ) with four ordered quaternion (( j , k ), ( j + 1, k ), ( j + 1, k + 1), ( j , k + 1)) by ☐ jk , and define X to be the set of all such squares. Then, the discrete jet bundle [ 23 , 41 ] is defined as follows:

which is equal to X × ℝ 4 .

Now, the discrete Lagrangian L d on is a discrete version of the Lagrangian functional, which is discretized by a finite difference method as follows:

where , , and . Similarly, we give the definitions of the discrete Lagrangian L d on the other three squares adjacent to ,

and

Given a finite subset U X of the space of squares, we form the discrete action sum as

Here, ψ is a discrete field, assigning to every node ( j , k ) a field value , and j 1 ψ is its first jet extension, defined by

We now focus on a particular configuration U , consisting of four adjacent squares, cf. Fig. 1 . The action sum for this U is explicitly given by

Fig. 1. Four adjacent squares with a common vertex .

By keeping the values of the field on the boundary fixed, and taking variations with respect to , we obtain the following discrete Euler–Lagrange equations: [ 23 , 24 , 41 ]

where D i means the derivative with respect to the i -th variable.

For nonlinear Schrödinger equations with variable coefficients ( 1 ) and the discrete Lagrangian L d (Eqs. ( 13 )–( 16 )), the discrete Euler–Lagrange equations ( 18 ) result in a second-order scheme

We discuss the convergence order of the proposed variational integrator ( 19 ). By the Taylor expansion, we have

From these equations, we can readily observe that the variational integrator (Eq. ( 19 )) has truncation error . To validate this conclusion, we investigate the numerical convergence order in numerical experiments.

We have mentioned in the Introduction section that, the advantages of deriving the multi-symplectic numerical schemes from discrete variational principle are that they are naturally multi-symplectic and the discrete multi-symplectic structures are also generated in the variational principle. It is meaningful to show the multi-symplectic structure of this discrete variational integrator using the definition of boundary Lagrangian.

Because of the square discretization, we here focus on four adjacent squares around , and denote this area by U . For the sake of simplicity, we denote the values of the field on the boundary by ψ ∂U , i.e.,

Along the line of Ref. [ 41 ], the discrete boundary Lagrangian L ∂U is given by

where the extremum is taken over all . Alternatively, L ∂U can be defined by solving the discrete Euler–Lagrange equations ( 18 ) for given the values of ψ ∂U as boundary data and substituting this value into the action sum ( 17 ).

We now derive the discrete multi-symplectic form formula by taking twice the exterior derivative of the boundary Lagrangian L ∂U . By taking the exterior derivative of both sides of Eq. ( 20 ) and using the definition of discrete Poincaré–Cartan [ 23 ] forms, we obtain

where discrete Cartan forms , , , and are defined by

and similarly for , , and . We can rewrite the discrete Euler–Lagrange equations ( 18 ) in terms of the discrete Cartan forms as

Using the discrete Euler–Lagrange equations ( 22 ), we can rewrite Eq. ( 21 ) as

and by taking another exterior derivative of both sides. Using that d 2 L ∂U ≡ 0, we finally get the discrete multisymplectic form formula with the following form: [ 41 ]

where discrete Poincaré–Cartan forms for n = 1, 2, 3, 4. This is the discrete version of Eq. ( 5 ) in Proposition 1. Numerical solutions of variational integrators ( 19 ) preserve the discrete multi-symplectic structure ( 23 ).

3. Numerical tests and simulations
3.1. Example 1: Cubic nonlinear Schrödinger equation

We consider a cubic nonlinear Schrödinger equation (cubic NLSE) ( 9 ) with parameters α ( t ) = cos( t )/2, β ( t ) = cos( t )/(sin( t ) + 3), and ν ( x ) = 0. The exact solution of this equation is given by [ 18 ]

We use the multi-symplectic variational integrator ( 19 ) to solve this problem over [−40, 40].

3.1.1. Convergence order tests

To investigate the numerical convergence of the proposed scheme ( 19 ), we conduct a series of numerical tests with varying mesh sizes. We take Δ t = 0.1Δ x to solve the cubic nonlinear Schrödinger equation ( 9 ), and compute errors at time t = 1 with different mesh sizes Δ x . In this test, we consider the l -error and l 2 -error that are defined by

The convergence order is calculated using the formula [ 44 ]

The l -error and l 2 -error at time t = 1 and the computed convergence orders are listed in Table 1 . It is clear that the error decreases as the mesh size goes to zero, indicating the convergence of our numerical scheme. Moreover, the numerical orders clearly exhibit second-order convergence when the mesh size decreases with fixing Δ t = 0.1Δ x .

Table 1.

The errors and convergence orders of the variational integrator ( 19 ) for the cubic NLSE ( 9 ) at time t = 1 with Δ t = 0.1 Δ x .

.

Figure 2 is the lg–lg plot of N ( N = 80/d x ) versus the l -error of variational integrator ( 19 ) for cubic NLSE ( 9 ). The number N varies from 200 to 900. The discrete points in the figure are close to the line with slope −2, showing that the order of accuracy of scheme ( 19 ) is second order for the cubic NLSE.

Fig. 2. Lg–lg plot of N versus l -error of integrator ( 19 ) for cubic NLSE ( 9 ).
3.1.2. Numerical simulations

The multi-symplectic variational integrator ( 19 ) is employed to compute the solution from t = 0 to t = 40 in [−40, 40] using Δ t = 0.1 and Δ x = 0.1. Figure 3 shows the numerical waveform of the solution. From the figure, one can observe that the multi-symplectic variational integrator ( 19 ) displays the numerical properties of the periodic solitary-wave clearly and precisely. In Fig. 4 , we show the l -error of numerical solution. The figure displays that the numerical solution has a good approximation to the exact solution.

Fig. 3. Numerical waveform of cubic NLSE by scheme ( 19 ) with Δ t = 0.1 and Δ x = 0.1.
Fig. 4. The l -error of integrator ( 19 ) for cubic NLSE with Δ t = 0.1 and Δ x = 0.1.

Now, we test the numerical conservation property of multi-symplectic variational integrator ( 19 ) for the mass conservation law ( 6 ) and momentum conservation law ( 7 ) of cubic NLSE. The discrete version of mass conservation law ( 6 ) can be written as

and the discrete version of momentum conservation law ( 7 )

where ℜ means the real part and 𝔍 is the imaginary part.

The numerical conservation performances of Eqs. ( 27 ) and ( 28 ) by the proposed variational integrator ( 19 ) are shown in Figs. 5 and 6 respectively. l -error of mass conservation and momentum conservation with step Δ t = 0.1 and Δ x = 0.1 are displayed in the two figures. We find that multi-symplectic variational integrator ( 19 ) preserves the mass conservation law and momentum conservation law pretty well with very small periodic oscillations.

Fig. 5. Error of mass conservation law in integrator ( 19 ) for cubic NLSE.
Fig. 6. Error of momentum conservation law in integrator ( 19 ) for cubic NLSE.
3.2. Example 2: Gross–Pitaevskii equation

For Gross–Pitaevskii equation (GPE) in one dimension, we can consider Eq. ( 1 ) with parameters α ( t ) = 1/2, β ( t ) = −1, and ν ( x ) = −cos 2 ( x ). For this equation, the exact solution [ 18 , 19 ] is given by

We solved this equation by the proposed multi-symplectic variational integrator ( 19 ) with homogeneous boundary conditions over [− π , π ].

3.2.1. Convergence order tests

To verify that the proposed scheme ( 19 ) has the expected convergence order in time and space, we first test the discretization error in space. We choose a very small time step Δ t = 0.001 such that the error from the time discretization is negligible compared to the spatial discretization error, and solve the GPE using the proposed method with various spatial mesh sizes Δ x .

Table 2 lists the numerical l -error and l 2 -error, respectively, at time t = 1 for various spatial mesh sizes Δ x . The convergence orders in space computed by formula ( 26 ) are also listed in Table 2 .

Table 2.

The convergence rate in time in the integrator ( 19 ) for GPE at time t = 1 with Δ t = 0.001.

.

We also test the discretization error in time. Table 3 shows the numerical l -error and l 2 -error at t = 1 with a very small space step size Δ x = 2 π /4096 for different time steps Δ t . The convergence order in time is computed and shown in this table.

Table 3.

The convergence rate in time in the integrator ( 19 ) for GPE at time t = 1 with Δ x = 2 π /4096.

.

From Tables 2 and 3 , one can observe that the numerical scheme ( 19 ) indeed has the expected order, i.e., second order in space and second order in time for GPE. The error decreases quickly when the mesh refines.

The lg–lg plot of N ( N = 2 π /d x ) versus the l -error for GPE is displayed in Fig. 7 . The number N varies from 20 to 170. The discrete points in Fig. 7 are close to the line with slope −2, showing that the order of accuracy of the variational integrator ( 19 ) is second-order for the GPE.

Fig. 7. Lg–lg plot of N versus l -error of integrator ( 19 ) for GPE.
3.2.2. Numerical simulations

The proposed multi-symplectic variational integrator ( 19 ) is used to solve GPE problem from t = 0 to t = 20 in [− π , π ] using Δ t = 0.01 and Δ x = 2 π /64. Figure 8 shows the numerical waveform of this problem. One can observe the multi-symplectic variational integrator ( 19 ) displays the waveform of GPE very well. In Fig. 9 , we show the l -error of numerical solutions.

Fig. 8. Numerical waveform of GPE by integrator ( 19 ) with Δ t = 0.01 and Δ x = 2 π /64.
Fig. 9. The l -error of scheme ( 19 ) for GPE with Δ t = 0.01 and Δ x = 2 π /64.

The exact solution of GPE preserves mass, momentum, and energy conservation laws. The discrete version of mass and momentum conservation law are defined as Eqs. ( 27 ) and ( 28 ). The discrete version of energy conservation law ( 8 ) is written as

The numerical performance of proposed scheme ( 19 ) for mass conservation ( 27 ), momentum conservation ( 28 ), and energy conservation ( 29 ) of GPE are shown in Figs. 10 12 . l -error of mass conservation, momentum conservation, and every conservation with step Δ t = 0.01 and Δ x = 2 π /64 are displayed in those three figures. Form the figures, one can observe that the multi-symplectic variational integrator ( 19 ) displays a good numerical behavior in preserving those three conservation laws. It preserves mass and energy conservation law pretty well with very small periodic oscillations.

Fig. 10. Error of mass conservation law in integrator ( 19 ) for GPE.
Fig. 11. Error of momentum conservation law in integrator ( 19 ) for GPE.
Fig. 12. Error of energy conservation law in integrator ( 19 ) for GPE.
4. Conclusions

We have proposed a second-order multi-symplectic variational integrator for nonlinear Schrödinger equations with variable coefficients. The natural multi-symplecticity is presented by a discrete multi-symplectic form formula, which is derived from the discrete Lagrangian boundary function. Cubic nonlinear Schrödinger equations and Gross–Pitaevskii equations have been extensively studied by the proposed integrator, as two examples of nonlinear Schrödinger equations with variable coefficients. Our numerical simulations demonstrate that the integrator is capable of preserving the mass, momentum, and energy conservation during time evolutions. Convergence tests are presented to verify the second-order accuracy both in time and space.

The multi-symplecticity variational integrator is efficient for solving nonlinear Schrödinger equations with variable coefficients. It leads to a good numerical performance in preserving intrinsic properties and conservation laws. We are going to study the second nonlinear Schrödinger equations and the coupled nonlinear Schrödinger equations for future work, try to extend the theory of variational error analysis to the setting of discrete field theories, and prove the convergence of the integrator.

Reference
1 Saleh B Teich M 1991 Fundamentals of Fotonics New York Wiley-Interscience
2 Stenflo L Yu M 1997 IEEE Trans. Plasma. Sci. 25 1155
3 Lu J 2004 Chin. Phys. B. 13 0811
4 Bao W Jaksch D Markowich P 2003 J. Comput. Phys. 187 318
5 Bao W Cai Y 2013 Kinet. Relat. Mod. 6 1
6 Liang X Khaliq A Sheng Q 2014 Appl. Math. Comput. 235 235
7 Sheng Q Khaliq A Al-Said E 2001 J. Comput. Phys. 166 400
8 Cai W Wang Y 2014 Chin. Phys. B 23 030204
9 Cai J Yang B Liang H 2013 Chin. Phys. B 22 030209
10 Cai J Wang J Wang Y 2015 Chin. Phys. B 24 050205
11 Cai J Wang Y 2013 Chin. Phys. B 22 060207
12 Cai J Qin Z Bai C 2013 Chin. Phys. Lett. 30 070202
13 Cai W Wang Y Song Y 2014 Chin. Phys. Lett. 31 040201
14 Wang Y Lv Z Zhang L 2014 Chin. Phys. B 23 120203
15 Lappas D L’Huillier A 1998 Phys. Rev. A 58 4140
16 Serkin V Hasegawa A 2000 Phys. Rev. Lett. 85 4502
17 Zhu H Chen Y Song S Hu H 2011 Appl. Numer. Math. 61 308
18 Dehghan M Taleei A 2010 CP Communications 181 43
19 Qian X Chen Y Song S 2014 Commun. Comput. Phys. 15 692
20 Bridges T Reich S 2001 Phys. Lett. A 284 184
21 Bridges T Reich S 2006 J. Phys. A: Math. Gen. 39 5287
22 Reich S 2000 J. Comput. Phys. 157 473
23 Marsden J Patrick G Shkoller S 1998 Commun. Math. Phys. 199 351
24 Marsden J West M 2001 Acta. Numerica 10 357
25 Hairer E Lubich C Wanner G 2006 Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations Berlin Springer
26 Hong J Liu Y Munthe-Kaas H Zanna A 2006 Appl. Numer. Math. 56 814
27 Bai D Zhang L 2011 Int. J. Comput. Math. 88 1714
28 Ran M Zhang C 2015 Int. J. Comput. Math.
29 Chen S Li Z Wang Y Chien C 2013 Int. J. Comput. Math. 90 651
30 Yuan Y Sun C Zhao F 2015 Acta Phys. Sin. 64 034202 (in Chinese)
31 Zhou S Cheng X 2010 Math. Comput. Simulat. 80 2362
32 Sheng Q Sun H 2014 Comput. Math. Appl. 68 1341
33 Deng D 2014 Int. J. Comput. Math. 91 1755
34 Hong J Liu Y 2003 Appl. Math. Lett. 16 759
35 Chen J Qin M Tang Y 2002 Comput. Math. Appl. 43 1095
36 Chen J Qin M 2001 ETNA 12 193
37 Leok M Zhang J 2011 IMA J. Numer. Anal. 31 1497
38 Liu S Liu C Guo Y 2011 Chin. Phys. B 20 034501
39 Chen J Qin M 2002 Numer. Meth. Part. D E 18 523
40 Chen J 2005 Appl. Math. Comput. 170 1394
41 Vankerschaver J Liao C Leok M 2013 J. Math. Phys. 54 082901
42 Sulem C Sulem P 1999 The Nonlinear Schrödinger Equation: Self-fousing and Wave Collapse New York Springer
43 Sacchetti A 2004 J. Evol. Equ. 4 345
44 Zhou S Cheng X 2011 Int. J. Comput. Math. 88 795