- Email: [email protected]
Symplectic integrators for discrete nonlinear Schrödinger systems
Mathematics and Computers in Simulation 56 (2001) 145–156
Symplectic integrators for discrete nonlinear Schrödinger systems D.A. Karpeev a , C.M. Schober b,∗ a
b
Department of Computer Science, Old Dominion University, Norfolk, VA 23529, USA Department of Mathematics and Statistics, Old Dominion University, Norfolk, VA 23529, USA
Abstract Symplectic methods for integrating canonical and non-canonical Hamiltonian systems are examined. A general form for higher order symplectic schemes is developed for non-canonical Hamiltonian systems using generating functions and is directly applied to the Ablowitz–Ladik discrete nonlinear Schrödinger system. The implicit midpoint scheme, which is symplectic for canonical systems, is applied to a standard Hamiltonian discretization. The symplectic integrators are compared with an explicit Runge–Kutta scheme of the same order. The relative performance of the integrators as the dimension of the system is varied is discussed. © 2001 IMACS. Published by Elsevier Science B.V. All rights reserved. Keywords: Symplectic integrators; Schrödinger systems; Hamiltonian systems
1. Introduction In recent years there has been wide interest in the use of symplectic methods for integrating Hamiltonian systems. Standard numerical schemes often neglect important features of the Hamiltonian dynamics in contrast with symplectic methods which are specifically designed to preserve the symplectic structure. While symplectic integrators have been shown to capture the dynamics of low-dimensional Hamiltonian systems very accurately for very long times, less has been established on their practicality for high-dimensional systems. To begin to address this issue, we examine the performance of symplectic integrators in solving two different discretizations of the nonlinear Schrödinger (NLS) system: iqt + qxx + 2q 2 p = 0,
−ipt + pxx + 2p 2 q = 0,
(1)
under periodic boundary conditions, (p(x + L, t), q(x + L, t)) = (p(x, t), q(x, t)), as the dimension of the system is varied. In modeling a variety of physically significant nonlinear phenomena, the condition p = ±q ∗ is frequently imposed and in this case the system reduces to the standard cubic NLS equation. Zakharov and ∗ Corresponding author. E-mail addresses: [email protected] (D.A. Karpeev), [email protected] (C.M. Schober).
0378-4754/01/$20.00 © 2001 IMACS. Published by Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 4 7 5 4 ( 0 1 ) 0 0 2 8 6 - 5
146
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
Shabat showed in 1972 that the cubic NLS equation is a completely integrable Hamiltonian system in the sense of the inverse scattering transform (IST) and thus possesses infinitely many integrals of motion [1]. In this paper, the conjugacy relation between p and q is not imposed and we consider two quite distinct lattice systems for integrating the NLS system. One discretization, the Ablowitz–Ladik (AL) system (Eq. (2)), inherits all the properties of the original PDE system, but it has a non-canonical Poisson bracket; the other, the DNLS (Eq. (5)) is simply a canonical Hamiltonian discretization. In addition to being used as numerical schemes, both discrete systems have important physical applications (see e.g. [8]). For a general Hamiltonian system, symplectic schemes can be constructed by noting that any transformation derived from a generating function is a symplectic map. In this paper we generalize the generating function method to include systems which are non-canonical. When the Hamiltonian is of potential form, an effective approach is to use the Hamilton–Jacobi equations to solve for the generating function [6]. Since the AL system does not have any such special structure to exploit, it is easier, as in [4,9], to work directly with the equations of motion to develop the appropriate generating function. On the other hand, the DNLS is canonical and standard symplectic integrators can be applied; we use the implicit midpoint rule, a member of the Gauss family of implicit Runge–Kutta methods which are symplectic [5]. It has been shown that the implicit midpoint method exhibits stability problems when dealing with stiff systems with symmetry [7]. One question then is whether the generating function symplectic schemes that we develop (which are implicit) will also develop stability problems. The symplectic schemes are compared with an explicit (non-symplectic) Runge–Kutta scheme of the same order. In doing so, it is important to bear in mind this additional issue of the explicit versus implicit nature of the schemes. Given this fact, it is significant that the symplectic schemes perform as well as they do. The paper is organized as follows. In Section 2 we present the equations under consideration along with their fundamental properties. In Section 3 we provide an overview of the generating function method and develop algorithms for constructing symplectic schemes for systems in a non-canonical form. The generating function symplectic scheme as well as standard time-stepping schemes, are provided in Section 4. The results of the numerical experiments comparing the symplectic and non-symplectic schemes on the bases of accuracy and preservation of the invariants is presented in Section 5. The relative performance of the integrators as the dimension of the system is varied is also discussed. 2. Lattice nonlinear Schrödinger systems The following semi-discrete version of the nonlinear Schrödinger system i
d qn−1 + qn+1 − 2qn qn + + pn qn (qn−1 + qn+1 ) = 0, dt h2 d pn−1 + pn+1 − 2pn −i pn + + pn qn (pn−1 + pn+1 ) = 0, dt h2
(2)
which we will denote as the AL system, was designed by Ablowitz and Ladik to possess an integrable structure equivalent to that of the continuous NLS system [2]. The Hamiltonian structure of the AL system is given by the Hamiltonian, N i 2 p, q ) = − 3 H (p [h pn (qn−1 + qn+1 ) − 2 ln(1 + h2 qn pn )], h n=1
(3)
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
147
p , q ) by together with the non-canonical Poisson brackets defined in coordinates (p {pm , qn } =
1 rn δm,n , h
{pm , pn } = {qm , qn } = 0,
(4)
where rn = 1 + h2 qn pn . Then the Hamiltonian equations of motion p , H }, p˙ = {p
q˙ = {qq , H }
are equivalent to Eq. (2). p, q ) = N A second invariant for system (2) is given by I (p n=1 pn (qn−1 + qn+1 ). There are N − 2 additional independent invariants (on the periodic interval or infinitely many on the infinite interval). As for the PDE system, when pn = ±qn∗ , it is possible to assign physical meaning to the invariants and derive the N-soliton solution for rapidly decreasing boundary conditions, pn , qn → 0 as n → ±∞, via the discrete IST [2] as well as quasi-periodic Riemnann theta function solutions for periodic boundary conditions (pj +N , qj +N ) = (pj , qj ) [3]. The second system we consider is the diagonal discrete NLS (DNLS) i
d qn−1 + qn+1 − 2qn qn + + 2qn2 pn = 0, dt h2
−i
d pn−1 + pn+1 − 2pn pn + + 2pn2 qn = 0, dt h2
(5)
which is a Hamiltonian system with respect to the canonical Poisson bracket {pm , qn } = δm,n ,
{pm , pn } = {qm , qn } = 0.
(6)
The DNLS possesses two constants of motion, the Hamiltonian p , q ) = −i H (p
N (pn+1 − pn )(qn+1 − qn ) n=1
h2
−
pn2 qn2
,
(7)
p, q ) = N and the L2 norm, I (p n=1 pn qn . Note that it can be proven that the solutions of AL and DNLS converge to those of the NLS PDE with O(h2 ) convergence as h → 0. Also, the scaled Hamiltonians H and the Poisson brackets {,} of AL and DNLS approach the Hamiltonian and the Poisson bracket, respectively, for the NLS PDE.
3. Symplectic discretizations of NLS lattice systems In this section we construct symplectic discretizations of AL and DNLS. The technique is the method for calculating the generating function of the exact flow to an arbitrary order. 3.1. Symplectic structure of Hamiltonian systems The phase space of any Hamiltonian system q˙n = {qn , H },
p˙ n = {pn , H }
(8)
148
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
with a non-degenerate Poisson bracket carries a natural symplectic structure given by dual to the bracket closed non-degenerate two-form ω. If {pm , pn } = {qm , qn } = 0, then the bracket tensor is 0 D , J = −D 0 and the corresponding symplectic form p, q ) = ω(p ωjk dpj ∧ dqk , has Ω as its tensor: 0 Ω ≡ (ωjk ) = D −1
−D −1 0
.
P (t), Q(t)) be the flow generated by the solution of Eq. (8) with the initial data (p p , q ), then the Poisson Let (P p , q ) → (P P , Q ): bracket and the corresponding form ω are invariant under the coordinate transformation (p p , q ) = ω(P P , Q ). ω(p
(9)
Since ω is closed, and thus locally exact, there is a primitive one-form θ, such that ω = dθ. Here θ is defined up to an exact differential so that θ = θ + dF is also a primitive of ω for any smooth F. Now Eq. (9) can be written as P , Q) = 0. p , q ) − dθ (P dθ(p P , Q ) is also closed and locally exact, therefore p , q ) − θ (P It follows that θ (p P , Q ) = dG p , q ) − θ (P θ(p
(10)
for some function G. In general, Eq. (10) characterizes any symplectic map and the function G is called the generating function of the transformation. P , Q ) for small values of t to obtain an In the case of a Hamiltonian flow, Eq. (10) can be solved for (P P (t), Q (t)). For the procedure to work the generating function explicit local coordinate representation of (P G must be known. Below we show how G can be obtained for AL, but the results readily generalize to other systems with similar symplectic structures. 3.2. Generating function of AL flow The Poisson bracket of AL is non-degenerate with the tensor 1 + h2 pn qn 0 D J = , D = diag , −D 0 h and the dual symplectic form ω=
N n=1
−
h dpn ∧ dqn . 1 + h 2 pn q n
(11)
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
149
Upon choosing two primitives, θ=
N ln(1 + h2 pn qn ) − dqn = f dqq , hq n n=1
θ =
fn
N ln(1 + h2 Pn Qn ) n=1
hPn −gn
P, dPn = −gg dP
the transformation equation (10) becomes P = dG, f dqq + g dP or in the component form ∂G p , q ), = fn (p ∂qn
∂G P , Q ). = gn (P ∂Pn
(12)
P (t), Q (t)) be the solution of AL with the initial data (p p , q ). The first equation of Eq. (12) is solved Let (P P ∂qq is non-singular. At for P locally using the implicit function theorem since the Jacobi matrix ∂ 2 G/∂P f dqq ) = −d(gg dp p ) is non-degenerate, and for t near zero the t = 0 this follows from the fact that ω = d(f P , q ) can be taken to be local coordinates in the neighborhood of same holds by continuity. Therefore, (P p , q ) and all expressions can be evaluated at (P P , q ). (p Observe that the right-hand side of the equations of motion for AL Rn ∂H P , Q ), P˙n = (P h ∂Qn
˙ n = − Rn ∂H (P P , Q ), Q h ∂Pn
Rn = 1 + h2 Pn Qn
(13)
P , Q ). This justifies local Taylor expansions in t is smooth (in fact analytic) in a C2N -neighborhood of (P P , q ). Likewise, smoothness of fn and gn and Eq. (12) imply for Qn (t), in a neighborhood of any point (P the same for G(t) so we have Qn (t) = qn +
∞ m t m=1
m!
P , q ), Qm,n (P
!qn
G(t) =
∞ m t m=0
m!
P , q ), Gm (P
(14)
whenever t is sufficiently small (and therefore P is close to p ). Substituting Qn = qn + !qn allows Rn , gn and the derivatives of H to be expanded in power series in P , q) t with all coefficients evaluated at (P Rn = 1 + h2 Pn qn +
R0,n
∞ s t
h2 Pn Qs,n ,
s! s=1
(15)
Rs,n
N s ∞ Qs,jk ls Q1,j1 l1 ∂ k+1 H 1 ∂H s = t ... , 1! s! ∂Pn ∂Pn (∂qj1 , . . . , ∂qjk ) l ≥0 l1 !, . . . , ls ! s=0 k=0 j ,... ,j =1 1
k
i ili =s li =k
Hs ,Pn
(16)
150
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
N s ∞ Qs,jk ls Q1 , j1 l1 ∂H ∂ k+1 H 1 s = t ... , ∂Qn ∂qn (∂qj1 , . . . , ∂qjk ) l ≥0 l1 !, . . . , ls ! 1! s! s=0 k=0 j ,... ,j =1 1
ili =s li =k
and
(17)
i
k
Hs,Qn
k ∞ s 1 s h2 Pn k gn = ln(1 + h Pn qn ) + t (−1) (k − 1)! hPn s=1 k=1 1 + h2 Pn qn 2
×
li ≥0 ili =s li =k
1 l1 !, . . . , ls !
Q1,n 1!
l1
Qm,n ... s!
ls
.
(18)
P , q ) the full time derivative of Qn is In the local coordinates (P N
∂Qn ˙ n = ∂Qn + Q P˙j . ∂t ∂P j j =1 Then solving for ∂Qn /∂t and using the equations of motion (13) yields N
∂Qn Rj ∂H ∂Qn Rn ∂H − . =− ∂t h ∂Pn j =1 ∂Pj h ∂Qj
(19)
Substituting Eqs. (15)–(17) and the power series for Qn into Eqs. (18) and (19) and the expansion for G into the second of the transformation equations, then equating powers of t yields a coupled recursive system ∂G0 1 =− ln(1 + h2 Pn qn ), ∂Pn hPn k m ∂Gm m! 1 h2 Pn Qm,n lm Q1,n l1 k = (−1) (k − 1)! ... , ∂Pn hPn k=1 1 + h2 Pn qn l !, . . . , lm ! 1! m! l ≥0 1
(20) (21)
i
ili =m li =k
and
Qm+1,n
N 1 1 1 ∂Qs,n m! Rs,n Hk,Pn + =− Rk,j Hl,Qj . h s+k=m s! s! k! j =1 ∂Pj s+k+l=m s,k≥0
(22)
s≥1 k,l≥0
Expressions for Gm are easily identifiable as full derivatives and solved. Then G is completely specified m P ,q) = ∞ P , q ). in the form G(P m=0 (t /m!)Gm (P
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
151
3.3. Symplectic schemes ˜ = Using a truncated expansion G p, q ) = fn (p
˜ ∂G (P˜ , q ), ∂qn
r
m=0 (t
m
/m!)Gm , equations
˜ ˜ ) = ∂ G (P˜ , q), gn (P˜ , Q ∂ P˜n
(23)
˜ ) = (P P , Q ) + O(t r+1 ), which follows from the fact that generate a coordinate transformation with (P˜ , Q ˜ ) agree to the same order. Therefore, (p ˜ ) is a P , Q ) and (P˜ , Q p , q ) → (P˜ , Q the implicit equations for (P symplectic map preserving dynamics up to rth order. 4. Numerical schemes In this section we describe symplectic schemes as well as other standard time integrators that are used in the numerical experiments. For AL a second-order symplectic scheme is derived as above using an approximate generating function ˜ = G0 (P P , q ) + tG1 (P P , q ) + 21 t 2 G2 (P P , q ), G and the transformation equation (23). Solving Eqs. (21) and (22) for Gm , m = 0, 1, 2 we obtain N qn 1 + h2 Pj qj ∂H ∂H 1 ln(1 + h2 Pn ξ ) dξ, G1 = H, G2 = − . G0 = − hP h ∂P ∂q n j j 0 j n=1 ˜ yields ˜ into Eq. (23) and solving for P˜ and Q Substituting the expression for G ∂En 1 ˜n = Q (1 + h2 P˜n qn )exp −hP˜n −1 , h2 P˜n ∂ P˜n ∂E t2 1 n 2 G2 . (1 + h − 1 , E p q )exp hq = tG + P˜n = n n n n 1 h2 qn ∂qn 2
(24)
Table 1 √ Hamiltonian deviation, AL: h = L/N , L = 2 2π , T = 500
S2 R2
N = 4; t = 10−2
N = 4; t = 10−3
N = 16; t = 10−2
N = 16; t = 10−3
N = 32; t = 10−3
N = 32; t = 10−4
N = 64; t = 5 × 10−4
N = 64; t = 10−4
7.1 × 10−6 1.3 × 10−5
7.1 × 10−8 3.5 × 10−8
2.7 × 10−4 2.2 × 10−4
2.7 × 10−6 5.2 × 10−7
3.5 × 10−6 6.0 × 10−7
5.1 × 10−8 4.1 × 10−9
9.6 × 10−7 1.2 × 10−7
9.9 × 10−8 4.1 × 10−9
Table 2 √ Hamiltonian deviation, DNLS: h = L/N , L = 2 2π, T = 500
S2 R2
N = 4; t = 10−2
N = 4; t = 10−3
N = 16; t = 10−2
N = 16; t = 10−3
N = 32; t = 10−3
N = 32; t = 10−4
N = 64; t = 5 × 10−4
N = 64; t = 10−4
1.7 × 10−5 1.6 × 10−4
1.7 × 10−7 2.7 × 10−7
2.4 × 10−4 9.4 × 10−2
3.4 × 10−6 1.0 × 10−4
1.5 × 10−6 2.5 × 10−6
1.6 × 10−8 1.6 × 10−8
7.2 × 10−7 7.2 × 10−7
3.0 × 10−8 2.3 × 10−8
152
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
For DNLS a second-order symplectic integrator is furnished by an implicit Runge–Kutta method or the implicit midpoint method (see e.g. [5]). It is defined as follows for a dynamical system: z˙ = F (zz ), Given the initial data z we compute an approximation Z˜ at time t as Z˜ = z + tF z + 21 tF(Z˜ ) .
(25)
Lastly, the standard explicit second-order Runge–Kutta scheme is defined, using the notation above, by Z˜ = z + tF(zz + 21 tF(zz )).
(26)
Fig. 1. AL Hamiltonian deviation, using a = 0.5 for N = 4 with (a) t = 10−2 , (b) t = 10−3 , and for N = 16 with (c) t = 10−2 , (d) t = 10−3 (left to right, top to bottom).
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
153
5. Numerical experiments In this section we examine the performance of symplectic methods in solving the AL and DNLS systems under periodic boundary conditions (pj +N , qj +N ) (pj , qj ). We numerically integrate the AL system (2) over the time interval T = 0–500 using the following schemes: (a) the implicit second-order symplectic scheme (24), denoted by S2 and (b) the explicit second-order Runge–Kutta scheme (26), denoted as R2. The DNLS system (5) is integrated over the time interval T = 0–500 using the following schemes: (a) the implicit midpoint method (25), denoted by SR2 and (b) the explicit Runge–Kutta scheme R2. Since the symplectic integrators S2 and SR2 are implicit, to advance the solution one time-step, a fixed point iteration procedure is used until the relative error is in the order of 10−12 . A fixed time-step is used in all schemes throughout the integration; we do not optimize the Runge–Kutta schemes with a variable step size. We consider initial data of the following form: qn (0) = pn∗ (0) = a(1 − - cos(bxn ))
(27) √ for xn = −L/2 +(n−1)h, h = L/N, n = 1, 2, . . . , N +1, where - = 10−2 , b = 2π/L and L = 2 2π , corresponding to a perturbation of the spatially uniform plane wave which is invariant under the phase flows. Note that although qn (0) = pn∗ (0) and the semi-discrete AL and DNLS flows preserve conjugacy, this condition is not imposed throughout the time evolution as it would violate the symplecticity of the numerical schemes. Even so, it is interesting to consider initial data of this form to determine which of the schemes minimizes the deviation from conjugacy, i.e. |pn − qn∗ |. We consider lattices with a varying number of mesh sites and evaluate the numerical schemes by monitoring the Hamiltonian H, the norm I and the amplitude of the waveform of the solution. In Tables 1 and 2 we present the maximum deviations in H, for a = 0.5 and T = 0–500, for mesh sizes N = 4, 16, 32 and 64, each for two step sizes, for AL and DNLS, respectively. The preservation of the second invariant I is not presented as it is qualitatively similar to H. In all of the numerical experiments with AL and DNLS, the symplectic integrators consistently preserved the invariants, as well as the waveform,
Fig. 2. AL amplitude of q1 using a = 0.5 for N = 16 with (a) t = 10−2 and (b) t = 10−3 .
154
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
in a more controlled fashion. However, as the dimension of the system is increased, the benefits of symplecticity are manifested on longer time-scales. Fig. 1a–d show the deviation in the Hamiltonian of the AL system with initial data (27), a = 0.5, obtained using S2 and R2 for N = 4 with (a) t = 10−2 , (b) t = 10−3 , and for N = 16 with (c) t = 10−2 , (d) t = 10−3 . The magnitude of the deviations in the Hamiltonian obtained using S2 can actually be slightly larger than that obtained with R2 (see Table 1). However, in all the cases examined, one notices that a linear drift in the error occurs with R2 which is not present for S2. Using the symplectic scheme, the error in the Hamiltonian oscillates in a bounded fashion and does not grow as it does with R2. For long time simulations the error obtained using R2 eventually exceeds that of S2 and R2 is thus sufficiently accurate on a shorter time-scale than S2. Even so, it is important to note that the linear drift in R2 becomes less significant as the dimension of the system is increased. When N is small, e.g. N = 4 and 16, the error in H obtained with R2 rapidly overtakes the error obtained with S2 (by T = 200 and 500, respectively). For larger N, the linear drift in H using R2 and the bounded oscillations with negligible growth in H using
Fig. 3. DNLS Hamiltonian deviation, using a = 0.5 for N = 4 with (a) t = 10−2 , (b) t = 10−3 , and for N = 32 with (c) t = 10−3 , (d) t = 10−4 .
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
155
Fig. 4. Conjugacy deviation |p1 − q1∗ |, using a = 0.5, N = 16 and t = 10−3 for (a) AL system and (b) DNLS system.
S2 are still observed, but it takes substantially longer for the error obtained with R2 to exceed S2 (see Table 1, N = 32 and 64). Fig. 2a and b show the amplitude of q1 of the solution obtained with the two integrators R2 and S2 using a = 0.5 for N = 16 with (a) t = 10−2 and (b) t = 10−3 . The amplitudes of the other lattice sites show similar qualitative behavior. Solutions of the AL system exhibit regular quasi-periodic motion due to the fact that the AL flow occurs in general on an N-torus. For t = 0.01, a phase lag develops using R2 which becomes more pronounced as the system evolves. Using t = 0.001 the solutions from the two integrators are virtually indistinguishable on the time-scale examined. A comparable relative performance between SR2 and R2 is observed for the DNLS system (see Fig. 3a–d and Table 2 for deviations in the Hamiltonian). Comparing Figs. 1 and 3, the difference in the performance of the symplectic and non-symplectic scheme is more significant than for the AL system. One feature which must be stressed is that for both systems R2 preserves the conjugacy relation, pn = qn∗ , much better than the symplectic schemes (see Fig. 4). For N = 16 the deviation in |p1 − q1∗ | is O(10−6 ) using S2, whereas with R2 it is in the order of round-off (the deviation in |pn − qn∗ | is comparable for general n). For both systems, stability issues are exhibited as can be seen from the N = 4 and 16 cases. Keeping the time-step fixed and varying N (equivalently h), as h decreases the performance of both schemes degrades. This suggests that t/ h < M, for some M, is required for stability. The instability is more pronounced for the explicit scheme R2 than for either of the symplectic schemes. It is surprising then that R2 preserves conjugacy better, indicating that instabilities of R2 lie in the p = q ∗ subspace, whereas for S2 and SR2 they are transverse to it.
Acknowledgements This work was partially supported by the NSF, Grant number DMS-9803567.
156
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
References [1] V.E. Zakharov, A.B. Shabat, Exact theory of two-dimensional self-focusing and one-dimensional self-modulation of waves in nonlinear media, Sov. Phys. JETP 34 (1972) 62–69. [2] M.J. Ablowitz, J.F. Ladik, A nonlinear difference scheme and inverse scattering, Stud. Appl. Math. 55 (1976) 213–229. [3] N.N. Bogolyubov, A.K. Prikarpatskii, The inverse periodic problem for a discrete approximation of a nonlinear Schrödinger equation, Sov. Phys. Dokl. 27 (1982) 412–416. [4] P.J. Channell, C. Scovel, Symplectic integration of Hamiltonian systems, Nonlinearity 3 (1990) 1–13. [5] J.M. Sanz-Serna, M.P. Calvo, Numerical Hamiltonian Problems, Chapman & Hall, London, 1994. [6] X. Lu, R. Schmid, Symplectic integration of Sine–Gordon type systems, Math. Comput. Simul. 50 (1999) 255–263. [7] O. Gonzalez, On the stability of symplectic and energy–momentum algorithms for non-linear Hamiltonian systems with symmetry, Comput. Methods Appl. Mech. Eng. 134 (1996) 197–222. [8] Ch. Claude, Yu.S. Kivshar, O. Kluth, K.H. Spatchek, Moving modes in localized nonlinear lattices, Phys. Rev. B 47 (1993) 14228–14232. [9] C.M. Schober, Symplectic integrators for the Ablowitz–Ladik discrete nonlinear Schrödinger equation, Phys. Lett. A 259 (1999) 140–151.
Symplectic integrators for discrete nonlinear Schrödinger systems D.A. Karpeev a , C.M. Schober b,∗ a
b
Department of Computer Science, Old Dominion University, Norfolk, VA 23529, USA Department of Mathematics and Statistics, Old Dominion University, Norfolk, VA 23529, USA
Abstract Symplectic methods for integrating canonical and non-canonical Hamiltonian systems are examined. A general form for higher order symplectic schemes is developed for non-canonical Hamiltonian systems using generating functions and is directly applied to the Ablowitz–Ladik discrete nonlinear Schrödinger system. The implicit midpoint scheme, which is symplectic for canonical systems, is applied to a standard Hamiltonian discretization. The symplectic integrators are compared with an explicit Runge–Kutta scheme of the same order. The relative performance of the integrators as the dimension of the system is varied is discussed. © 2001 IMACS. Published by Elsevier Science B.V. All rights reserved. Keywords: Symplectic integrators; Schrödinger systems; Hamiltonian systems
1. Introduction In recent years there has been wide interest in the use of symplectic methods for integrating Hamiltonian systems. Standard numerical schemes often neglect important features of the Hamiltonian dynamics in contrast with symplectic methods which are specifically designed to preserve the symplectic structure. While symplectic integrators have been shown to capture the dynamics of low-dimensional Hamiltonian systems very accurately for very long times, less has been established on their practicality for high-dimensional systems. To begin to address this issue, we examine the performance of symplectic integrators in solving two different discretizations of the nonlinear Schrödinger (NLS) system: iqt + qxx + 2q 2 p = 0,
−ipt + pxx + 2p 2 q = 0,
(1)
under periodic boundary conditions, (p(x + L, t), q(x + L, t)) = (p(x, t), q(x, t)), as the dimension of the system is varied. In modeling a variety of physically significant nonlinear phenomena, the condition p = ±q ∗ is frequently imposed and in this case the system reduces to the standard cubic NLS equation. Zakharov and ∗ Corresponding author. E-mail addresses: [email protected] (D.A. Karpeev), [email protected] (C.M. Schober).
0378-4754/01/$20.00 © 2001 IMACS. Published by Elsevier Science B.V. All rights reserved. PII: S 0 3 7 8 - 4 7 5 4 ( 0 1 ) 0 0 2 8 6 - 5
146
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
Shabat showed in 1972 that the cubic NLS equation is a completely integrable Hamiltonian system in the sense of the inverse scattering transform (IST) and thus possesses infinitely many integrals of motion [1]. In this paper, the conjugacy relation between p and q is not imposed and we consider two quite distinct lattice systems for integrating the NLS system. One discretization, the Ablowitz–Ladik (AL) system (Eq. (2)), inherits all the properties of the original PDE system, but it has a non-canonical Poisson bracket; the other, the DNLS (Eq. (5)) is simply a canonical Hamiltonian discretization. In addition to being used as numerical schemes, both discrete systems have important physical applications (see e.g. [8]). For a general Hamiltonian system, symplectic schemes can be constructed by noting that any transformation derived from a generating function is a symplectic map. In this paper we generalize the generating function method to include systems which are non-canonical. When the Hamiltonian is of potential form, an effective approach is to use the Hamilton–Jacobi equations to solve for the generating function [6]. Since the AL system does not have any such special structure to exploit, it is easier, as in [4,9], to work directly with the equations of motion to develop the appropriate generating function. On the other hand, the DNLS is canonical and standard symplectic integrators can be applied; we use the implicit midpoint rule, a member of the Gauss family of implicit Runge–Kutta methods which are symplectic [5]. It has been shown that the implicit midpoint method exhibits stability problems when dealing with stiff systems with symmetry [7]. One question then is whether the generating function symplectic schemes that we develop (which are implicit) will also develop stability problems. The symplectic schemes are compared with an explicit (non-symplectic) Runge–Kutta scheme of the same order. In doing so, it is important to bear in mind this additional issue of the explicit versus implicit nature of the schemes. Given this fact, it is significant that the symplectic schemes perform as well as they do. The paper is organized as follows. In Section 2 we present the equations under consideration along with their fundamental properties. In Section 3 we provide an overview of the generating function method and develop algorithms for constructing symplectic schemes for systems in a non-canonical form. The generating function symplectic scheme as well as standard time-stepping schemes, are provided in Section 4. The results of the numerical experiments comparing the symplectic and non-symplectic schemes on the bases of accuracy and preservation of the invariants is presented in Section 5. The relative performance of the integrators as the dimension of the system is varied is also discussed. 2. Lattice nonlinear Schrödinger systems The following semi-discrete version of the nonlinear Schrödinger system i
d qn−1 + qn+1 − 2qn qn + + pn qn (qn−1 + qn+1 ) = 0, dt h2 d pn−1 + pn+1 − 2pn −i pn + + pn qn (pn−1 + pn+1 ) = 0, dt h2
(2)
which we will denote as the AL system, was designed by Ablowitz and Ladik to possess an integrable structure equivalent to that of the continuous NLS system [2]. The Hamiltonian structure of the AL system is given by the Hamiltonian, N i 2 p, q ) = − 3 H (p [h pn (qn−1 + qn+1 ) − 2 ln(1 + h2 qn pn )], h n=1
(3)
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
147
p , q ) by together with the non-canonical Poisson brackets defined in coordinates (p {pm , qn } =
1 rn δm,n , h
{pm , pn } = {qm , qn } = 0,
(4)
where rn = 1 + h2 qn pn . Then the Hamiltonian equations of motion p , H }, p˙ = {p
q˙ = {qq , H }
are equivalent to Eq. (2). p, q ) = N A second invariant for system (2) is given by I (p n=1 pn (qn−1 + qn+1 ). There are N − 2 additional independent invariants (on the periodic interval or infinitely many on the infinite interval). As for the PDE system, when pn = ±qn∗ , it is possible to assign physical meaning to the invariants and derive the N-soliton solution for rapidly decreasing boundary conditions, pn , qn → 0 as n → ±∞, via the discrete IST [2] as well as quasi-periodic Riemnann theta function solutions for periodic boundary conditions (pj +N , qj +N ) = (pj , qj ) [3]. The second system we consider is the diagonal discrete NLS (DNLS) i
d qn−1 + qn+1 − 2qn qn + + 2qn2 pn = 0, dt h2
−i
d pn−1 + pn+1 − 2pn pn + + 2pn2 qn = 0, dt h2
(5)
which is a Hamiltonian system with respect to the canonical Poisson bracket {pm , qn } = δm,n ,
{pm , pn } = {qm , qn } = 0.
(6)
The DNLS possesses two constants of motion, the Hamiltonian p , q ) = −i H (p
N (pn+1 − pn )(qn+1 − qn ) n=1
h2
−
pn2 qn2
,
(7)
p, q ) = N and the L2 norm, I (p n=1 pn qn . Note that it can be proven that the solutions of AL and DNLS converge to those of the NLS PDE with O(h2 ) convergence as h → 0. Also, the scaled Hamiltonians H and the Poisson brackets {,} of AL and DNLS approach the Hamiltonian and the Poisson bracket, respectively, for the NLS PDE.
3. Symplectic discretizations of NLS lattice systems In this section we construct symplectic discretizations of AL and DNLS. The technique is the method for calculating the generating function of the exact flow to an arbitrary order. 3.1. Symplectic structure of Hamiltonian systems The phase space of any Hamiltonian system q˙n = {qn , H },
p˙ n = {pn , H }
(8)
148
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
with a non-degenerate Poisson bracket carries a natural symplectic structure given by dual to the bracket closed non-degenerate two-form ω. If {pm , pn } = {qm , qn } = 0, then the bracket tensor is 0 D , J = −D 0 and the corresponding symplectic form p, q ) = ω(p ωjk dpj ∧ dqk , has Ω as its tensor: 0 Ω ≡ (ωjk ) = D −1
−D −1 0
.
P (t), Q(t)) be the flow generated by the solution of Eq. (8) with the initial data (p p , q ), then the Poisson Let (P p , q ) → (P P , Q ): bracket and the corresponding form ω are invariant under the coordinate transformation (p p , q ) = ω(P P , Q ). ω(p
(9)
Since ω is closed, and thus locally exact, there is a primitive one-form θ, such that ω = dθ. Here θ is defined up to an exact differential so that θ = θ + dF is also a primitive of ω for any smooth F. Now Eq. (9) can be written as P , Q) = 0. p , q ) − dθ (P dθ(p P , Q ) is also closed and locally exact, therefore p , q ) − θ (P It follows that θ (p P , Q ) = dG p , q ) − θ (P θ(p
(10)
for some function G. In general, Eq. (10) characterizes any symplectic map and the function G is called the generating function of the transformation. P , Q ) for small values of t to obtain an In the case of a Hamiltonian flow, Eq. (10) can be solved for (P P (t), Q (t)). For the procedure to work the generating function explicit local coordinate representation of (P G must be known. Below we show how G can be obtained for AL, but the results readily generalize to other systems with similar symplectic structures. 3.2. Generating function of AL flow The Poisson bracket of AL is non-degenerate with the tensor 1 + h2 pn qn 0 D J = , D = diag , −D 0 h and the dual symplectic form ω=
N n=1
−
h dpn ∧ dqn . 1 + h 2 pn q n
(11)
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
149
Upon choosing two primitives, θ=
N ln(1 + h2 pn qn ) − dqn = f dqq , hq n n=1
θ =
fn
N ln(1 + h2 Pn Qn ) n=1
hPn −gn
P, dPn = −gg dP
the transformation equation (10) becomes P = dG, f dqq + g dP or in the component form ∂G p , q ), = fn (p ∂qn
∂G P , Q ). = gn (P ∂Pn
(12)
P (t), Q (t)) be the solution of AL with the initial data (p p , q ). The first equation of Eq. (12) is solved Let (P P ∂qq is non-singular. At for P locally using the implicit function theorem since the Jacobi matrix ∂ 2 G/∂P f dqq ) = −d(gg dp p ) is non-degenerate, and for t near zero the t = 0 this follows from the fact that ω = d(f P , q ) can be taken to be local coordinates in the neighborhood of same holds by continuity. Therefore, (P p , q ) and all expressions can be evaluated at (P P , q ). (p Observe that the right-hand side of the equations of motion for AL Rn ∂H P , Q ), P˙n = (P h ∂Qn
˙ n = − Rn ∂H (P P , Q ), Q h ∂Pn
Rn = 1 + h2 Pn Qn
(13)
P , Q ). This justifies local Taylor expansions in t is smooth (in fact analytic) in a C2N -neighborhood of (P P , q ). Likewise, smoothness of fn and gn and Eq. (12) imply for Qn (t), in a neighborhood of any point (P the same for G(t) so we have Qn (t) = qn +
∞ m t m=1
m!
P , q ), Qm,n (P
!qn
G(t) =
∞ m t m=0
m!
P , q ), Gm (P
(14)
whenever t is sufficiently small (and therefore P is close to p ). Substituting Qn = qn + !qn allows Rn , gn and the derivatives of H to be expanded in power series in P , q) t with all coefficients evaluated at (P Rn = 1 + h2 Pn qn +
R0,n
∞ s t
h2 Pn Qs,n ,
s! s=1
(15)
Rs,n
N s ∞ Qs,jk ls Q1,j1 l1 ∂ k+1 H 1 ∂H s = t ... , 1! s! ∂Pn ∂Pn (∂qj1 , . . . , ∂qjk ) l ≥0 l1 !, . . . , ls ! s=0 k=0 j ,... ,j =1 1
k
i ili =s li =k
Hs ,Pn
(16)
150
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
N s ∞ Qs,jk ls Q1 , j1 l1 ∂H ∂ k+1 H 1 s = t ... , ∂Qn ∂qn (∂qj1 , . . . , ∂qjk ) l ≥0 l1 !, . . . , ls ! 1! s! s=0 k=0 j ,... ,j =1 1
ili =s li =k
and
(17)
i
k
Hs,Qn
k ∞ s 1 s h2 Pn k gn = ln(1 + h Pn qn ) + t (−1) (k − 1)! hPn s=1 k=1 1 + h2 Pn qn 2
×
li ≥0 ili =s li =k
1 l1 !, . . . , ls !
Q1,n 1!
l1
Qm,n ... s!
ls
.
(18)
P , q ) the full time derivative of Qn is In the local coordinates (P N
∂Qn ˙ n = ∂Qn + Q P˙j . ∂t ∂P j j =1 Then solving for ∂Qn /∂t and using the equations of motion (13) yields N
∂Qn Rj ∂H ∂Qn Rn ∂H − . =− ∂t h ∂Pn j =1 ∂Pj h ∂Qj
(19)
Substituting Eqs. (15)–(17) and the power series for Qn into Eqs. (18) and (19) and the expansion for G into the second of the transformation equations, then equating powers of t yields a coupled recursive system ∂G0 1 =− ln(1 + h2 Pn qn ), ∂Pn hPn k m ∂Gm m! 1 h2 Pn Qm,n lm Q1,n l1 k = (−1) (k − 1)! ... , ∂Pn hPn k=1 1 + h2 Pn qn l !, . . . , lm ! 1! m! l ≥0 1
(20) (21)
i
ili =m li =k
and
Qm+1,n
N 1 1 1 ∂Qs,n m! Rs,n Hk,Pn + =− Rk,j Hl,Qj . h s+k=m s! s! k! j =1 ∂Pj s+k+l=m s,k≥0
(22)
s≥1 k,l≥0
Expressions for Gm are easily identifiable as full derivatives and solved. Then G is completely specified m P ,q) = ∞ P , q ). in the form G(P m=0 (t /m!)Gm (P
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
151
3.3. Symplectic schemes ˜ = Using a truncated expansion G p, q ) = fn (p
˜ ∂G (P˜ , q ), ∂qn
r
m=0 (t
m
/m!)Gm , equations
˜ ˜ ) = ∂ G (P˜ , q), gn (P˜ , Q ∂ P˜n
(23)
˜ ) = (P P , Q ) + O(t r+1 ), which follows from the fact that generate a coordinate transformation with (P˜ , Q ˜ ) agree to the same order. Therefore, (p ˜ ) is a P , Q ) and (P˜ , Q p , q ) → (P˜ , Q the implicit equations for (P symplectic map preserving dynamics up to rth order. 4. Numerical schemes In this section we describe symplectic schemes as well as other standard time integrators that are used in the numerical experiments. For AL a second-order symplectic scheme is derived as above using an approximate generating function ˜ = G0 (P P , q ) + tG1 (P P , q ) + 21 t 2 G2 (P P , q ), G and the transformation equation (23). Solving Eqs. (21) and (22) for Gm , m = 0, 1, 2 we obtain N qn 1 + h2 Pj qj ∂H ∂H 1 ln(1 + h2 Pn ξ ) dξ, G1 = H, G2 = − . G0 = − hP h ∂P ∂q n j j 0 j n=1 ˜ yields ˜ into Eq. (23) and solving for P˜ and Q Substituting the expression for G ∂En 1 ˜n = Q (1 + h2 P˜n qn )exp −hP˜n −1 , h2 P˜n ∂ P˜n ∂E t2 1 n 2 G2 . (1 + h − 1 , E p q )exp hq = tG + P˜n = n n n n 1 h2 qn ∂qn 2
(24)
Table 1 √ Hamiltonian deviation, AL: h = L/N , L = 2 2π , T = 500
S2 R2
N = 4; t = 10−2
N = 4; t = 10−3
N = 16; t = 10−2
N = 16; t = 10−3
N = 32; t = 10−3
N = 32; t = 10−4
N = 64; t = 5 × 10−4
N = 64; t = 10−4
7.1 × 10−6 1.3 × 10−5
7.1 × 10−8 3.5 × 10−8
2.7 × 10−4 2.2 × 10−4
2.7 × 10−6 5.2 × 10−7
3.5 × 10−6 6.0 × 10−7
5.1 × 10−8 4.1 × 10−9
9.6 × 10−7 1.2 × 10−7
9.9 × 10−8 4.1 × 10−9
Table 2 √ Hamiltonian deviation, DNLS: h = L/N , L = 2 2π, T = 500
S2 R2
N = 4; t = 10−2
N = 4; t = 10−3
N = 16; t = 10−2
N = 16; t = 10−3
N = 32; t = 10−3
N = 32; t = 10−4
N = 64; t = 5 × 10−4
N = 64; t = 10−4
1.7 × 10−5 1.6 × 10−4
1.7 × 10−7 2.7 × 10−7
2.4 × 10−4 9.4 × 10−2
3.4 × 10−6 1.0 × 10−4
1.5 × 10−6 2.5 × 10−6
1.6 × 10−8 1.6 × 10−8
7.2 × 10−7 7.2 × 10−7
3.0 × 10−8 2.3 × 10−8
152
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
For DNLS a second-order symplectic integrator is furnished by an implicit Runge–Kutta method or the implicit midpoint method (see e.g. [5]). It is defined as follows for a dynamical system: z˙ = F (zz ), Given the initial data z we compute an approximation Z˜ at time t as Z˜ = z + tF z + 21 tF(Z˜ ) .
(25)
Lastly, the standard explicit second-order Runge–Kutta scheme is defined, using the notation above, by Z˜ = z + tF(zz + 21 tF(zz )).
(26)
Fig. 1. AL Hamiltonian deviation, using a = 0.5 for N = 4 with (a) t = 10−2 , (b) t = 10−3 , and for N = 16 with (c) t = 10−2 , (d) t = 10−3 (left to right, top to bottom).
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
153
5. Numerical experiments In this section we examine the performance of symplectic methods in solving the AL and DNLS systems under periodic boundary conditions (pj +N , qj +N ) (pj , qj ). We numerically integrate the AL system (2) over the time interval T = 0–500 using the following schemes: (a) the implicit second-order symplectic scheme (24), denoted by S2 and (b) the explicit second-order Runge–Kutta scheme (26), denoted as R2. The DNLS system (5) is integrated over the time interval T = 0–500 using the following schemes: (a) the implicit midpoint method (25), denoted by SR2 and (b) the explicit Runge–Kutta scheme R2. Since the symplectic integrators S2 and SR2 are implicit, to advance the solution one time-step, a fixed point iteration procedure is used until the relative error is in the order of 10−12 . A fixed time-step is used in all schemes throughout the integration; we do not optimize the Runge–Kutta schemes with a variable step size. We consider initial data of the following form: qn (0) = pn∗ (0) = a(1 − - cos(bxn ))
(27) √ for xn = −L/2 +(n−1)h, h = L/N, n = 1, 2, . . . , N +1, where - = 10−2 , b = 2π/L and L = 2 2π , corresponding to a perturbation of the spatially uniform plane wave which is invariant under the phase flows. Note that although qn (0) = pn∗ (0) and the semi-discrete AL and DNLS flows preserve conjugacy, this condition is not imposed throughout the time evolution as it would violate the symplecticity of the numerical schemes. Even so, it is interesting to consider initial data of this form to determine which of the schemes minimizes the deviation from conjugacy, i.e. |pn − qn∗ |. We consider lattices with a varying number of mesh sites and evaluate the numerical schemes by monitoring the Hamiltonian H, the norm I and the amplitude of the waveform of the solution. In Tables 1 and 2 we present the maximum deviations in H, for a = 0.5 and T = 0–500, for mesh sizes N = 4, 16, 32 and 64, each for two step sizes, for AL and DNLS, respectively. The preservation of the second invariant I is not presented as it is qualitatively similar to H. In all of the numerical experiments with AL and DNLS, the symplectic integrators consistently preserved the invariants, as well as the waveform,
Fig. 2. AL amplitude of q1 using a = 0.5 for N = 16 with (a) t = 10−2 and (b) t = 10−3 .
154
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
in a more controlled fashion. However, as the dimension of the system is increased, the benefits of symplecticity are manifested on longer time-scales. Fig. 1a–d show the deviation in the Hamiltonian of the AL system with initial data (27), a = 0.5, obtained using S2 and R2 for N = 4 with (a) t = 10−2 , (b) t = 10−3 , and for N = 16 with (c) t = 10−2 , (d) t = 10−3 . The magnitude of the deviations in the Hamiltonian obtained using S2 can actually be slightly larger than that obtained with R2 (see Table 1). However, in all the cases examined, one notices that a linear drift in the error occurs with R2 which is not present for S2. Using the symplectic scheme, the error in the Hamiltonian oscillates in a bounded fashion and does not grow as it does with R2. For long time simulations the error obtained using R2 eventually exceeds that of S2 and R2 is thus sufficiently accurate on a shorter time-scale than S2. Even so, it is important to note that the linear drift in R2 becomes less significant as the dimension of the system is increased. When N is small, e.g. N = 4 and 16, the error in H obtained with R2 rapidly overtakes the error obtained with S2 (by T = 200 and 500, respectively). For larger N, the linear drift in H using R2 and the bounded oscillations with negligible growth in H using
Fig. 3. DNLS Hamiltonian deviation, using a = 0.5 for N = 4 with (a) t = 10−2 , (b) t = 10−3 , and for N = 32 with (c) t = 10−3 , (d) t = 10−4 .
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
155
Fig. 4. Conjugacy deviation |p1 − q1∗ |, using a = 0.5, N = 16 and t = 10−3 for (a) AL system and (b) DNLS system.
S2 are still observed, but it takes substantially longer for the error obtained with R2 to exceed S2 (see Table 1, N = 32 and 64). Fig. 2a and b show the amplitude of q1 of the solution obtained with the two integrators R2 and S2 using a = 0.5 for N = 16 with (a) t = 10−2 and (b) t = 10−3 . The amplitudes of the other lattice sites show similar qualitative behavior. Solutions of the AL system exhibit regular quasi-periodic motion due to the fact that the AL flow occurs in general on an N-torus. For t = 0.01, a phase lag develops using R2 which becomes more pronounced as the system evolves. Using t = 0.001 the solutions from the two integrators are virtually indistinguishable on the time-scale examined. A comparable relative performance between SR2 and R2 is observed for the DNLS system (see Fig. 3a–d and Table 2 for deviations in the Hamiltonian). Comparing Figs. 1 and 3, the difference in the performance of the symplectic and non-symplectic scheme is more significant than for the AL system. One feature which must be stressed is that for both systems R2 preserves the conjugacy relation, pn = qn∗ , much better than the symplectic schemes (see Fig. 4). For N = 16 the deviation in |p1 − q1∗ | is O(10−6 ) using S2, whereas with R2 it is in the order of round-off (the deviation in |pn − qn∗ | is comparable for general n). For both systems, stability issues are exhibited as can be seen from the N = 4 and 16 cases. Keeping the time-step fixed and varying N (equivalently h), as h decreases the performance of both schemes degrades. This suggests that t/ h < M, for some M, is required for stability. The instability is more pronounced for the explicit scheme R2 than for either of the symplectic schemes. It is surprising then that R2 preserves conjugacy better, indicating that instabilities of R2 lie in the p = q ∗ subspace, whereas for S2 and SR2 they are transverse to it.
Acknowledgements This work was partially supported by the NSF, Grant number DMS-9803567.
156
D.A. Karpeev, C.M. Schober / Mathematics and Computers in Simulation 56 (2001) 145–156
References [1] V.E. Zakharov, A.B. Shabat, Exact theory of two-dimensional self-focusing and one-dimensional self-modulation of waves in nonlinear media, Sov. Phys. JETP 34 (1972) 62–69. [2] M.J. Ablowitz, J.F. Ladik, A nonlinear difference scheme and inverse scattering, Stud. Appl. Math. 55 (1976) 213–229. [3] N.N. Bogolyubov, A.K. Prikarpatskii, The inverse periodic problem for a discrete approximation of a nonlinear Schrödinger equation, Sov. Phys. Dokl. 27 (1982) 412–416. [4] P.J. Channell, C. Scovel, Symplectic integration of Hamiltonian systems, Nonlinearity 3 (1990) 1–13. [5] J.M. Sanz-Serna, M.P. Calvo, Numerical Hamiltonian Problems, Chapman & Hall, London, 1994. [6] X. Lu, R. Schmid, Symplectic integration of Sine–Gordon type systems, Math. Comput. Simul. 50 (1999) 255–263. [7] O. Gonzalez, On the stability of symplectic and energy–momentum algorithms for non-linear Hamiltonian systems with symmetry, Comput. Methods Appl. Mech. Eng. 134 (1996) 197–222. [8] Ch. Claude, Yu.S. Kivshar, O. Kluth, K.H. Spatchek, Moving modes in localized nonlinear lattices, Phys. Rev. B 47 (1993) 14228–14232. [9] C.M. Schober, Symplectic integrators for the Ablowitz–Ladik discrete nonlinear Schrödinger equation, Phys. Lett. A 259 (1999) 140–151.