- Email: [email protected]
A multisymplectic integrator for the periodic nonlinear Schrödinger equation
Applied Mathematics and Computation 170 (2005) 1394–1417
www.elsevier.com/locate/amc
A multisymplectic integrator for the periodic nonlinear Schro¨dinger equation Jing-Bo Chen Institute of Geology and Geophysics, Chinese Academy of Sciences, P.O. Box 9825, Beijing 100029, PR China
Abstract A multisymplectic integrator for the periodic nonlinear Schro¨dinger equation is presented in this paper. Its accuracy is proved. By introducing a norm, we investigate its nonlinear stability. We also discuss the relationship between this multisymplectic integrator and two variational integrators which are derived by using the discrete multisymplectic field theory and the finite element method. 2005 Elsevier Inc. All rights reserved. Keywords: Multisymplectic integrator; Nonlinear Schro¨dinger equation
1. Introduction We consider the nonlinear Schro¨dinger equation (NLS) written in the form 2 iwt þ wxx þ V 0 jwj w ¼ 0; ð1:1Þ
E-mail address: [email protected] 0096-3003/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.01.031
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
1395
with the periodic boundary condition w(x + L,t) = w(x, t). Here L is the period and V : R ! R is a differentiable function. We will consider the space domain [0, L] in this paper. Eq. (1.1) is one of the most important completely integrable models in the theory of solitons. Its application can be found in many areas of physics, including nonlinear optics and plasma physics [14,17]. Various numerical methods for (1.1) have been developed. See for example [1,7,8,12,15,19,21–23]. Among these numerical methods, one kind of methods are so-called symplectic methods, which take the symplectic structure of the NLS equation into account when discretizing. The basic idea is to find a finite dimensional spatial truncation using finite difference , finite element or spectral methods, so that the resulting semidiscrete system can be cast into a finite dimensional Hamiltonian system. Then we integrate the finite dimensional Hamiltonian system in time using symplectic integrators [15,19,22,23]. The symplectic methods preserve a symplectic form which is summed over the discrete space domain. Recently, Bridges and Reich introduced the concept of multisymplectic integrators for Hamiltonian partial differential equations [4]. Hamiltonian partial differential equations can be formulated as multisymplectic Hamiltonian systems which possess a local multisymplectic conservation law. The multisymplectic integrators preserve a discrete multisymplectic conservation law. We know that the symplectic integrators for partial differential equations preserve a global symplectic form which is summed over the discrete spatial domain, i.e., the local variations in symplecticity over space domain are averaged. Therefore, the symplectic integrators do not reflect the spatial distribution of symplecticity. The multisymplectic integrators overcome the limitations of symplectic integrators since they preserve a local multisymplectic conservation law. For multisymplectic integrators, see also [5–8,21]. The NLS equation can be reformulated as a multisymplectic Hamiltonian system [3]. Reich developed multisymplectic integrators for the NLS equation in [21]. Numerical experiments were performed with a multisymplectic integrator for the NLS equation in [7], which demonstrate the good performance and local conservation properties of the multisymplectic integrator. However, the accuracy and nonlinear stability of the multisymplectic integrator employed in [7] were not presented. The purpose of this paper is to investigate its accuracy and stability as well as its relationship with two variational integrators. This paper is organized as follows. In Section 2, we present a multisymplectic integrator for the periodic NLS equation. In Section 3, we investigate its accuracy and stability. In Section 4, we present the multisymplectic geometry for the NLS equation. We explore the relationship between the multisymplectic integrator and two variational integrators in Section 5. Section 6 is devoted to reporting the numerical experiments. We end this paper with some conclusions in Section 7.
1396
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
2. A multisymplectic integrator for the periodic NLS equation In this section, we will present a multisymplectic integrator for the periodic NLS equation. We begin by introducing the multisymplectic formulation of the NLS equation. Using w = p + iq, we can rewrite (1.1) as a pair of real-valued equations pt þ qxx þ V 0 p2 þ q2 q ¼ 0; ð2:1Þ qt pxx V 0 p2 þ q2 p ¼ 0: Introducing a pair of conjugate momenta v = px, w = qx, the pair of equations (2.1) can be reformulated as the following first-order system qt vx ¼ V 0 ðp2 þ q2 Þp; pt wx ¼ V 0 ðp2 þ q2 Þq;
ð2:2Þ
px ¼ v; qx ¼ w; which can be written as a multisymplectic Hamiltonian system Mzt þ Kzx ¼ rSðzÞ;
ð2:3Þ
where SðzÞ ¼ 12 v2 þ w2 þ V ðp2 þ q2 Þ ;
z ¼ ðp; q; v; wÞT ; and 0
0
B 1 B M ¼B @ 0 0
0
1
1
0 0
0 0
0 0C C C; 0 0A
0
0 0
0
0
1
B0 B K¼B @1
0 0
0 0
0
1
0
0
1
1 C C C: 0 A 0
By direct calculations, we can prove that the system (2.3) satisfies the multisymplectic conservation law o o ðdp ^ dqÞ þ ðdv ^ dp þ dw ^ dqÞ ¼ 0: ot ox
ð2:4Þ
Now we discretize (2.2) using mid-point scheme in both time and space directions and obtain qmþ12;nþ1 qmþ12;n Dt
pmþ12;nþ1 pmþ12;n Dt
vmþ1;nþ12 vm;nþ12 Dx
¼ V 1;
wmþ1;nþ12 wm;nþ12 Dx
¼ V 2;
ð2:5aÞ ð2:5bÞ
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
pmþ1;nþ12 pm;nþ12 Dx qmþ1;nþ12 qm;nþ12 Dx
1397
¼ vmþ12;nþ12 ;
ð2:5cÞ
¼ wnþ12;mþ12 ;
ð2:5dÞ
where Dx and Dt are the space step size and time step size respectively and qmþ12;nþ1 ¼ 12ðqm;nþ1 þ qmþ1;nþ1 Þ;
pm;n pðmDx; nDtÞ;
vmþ12;nþ12 ¼ 14ðvm;n þ vmþ1;n þ vm;nþ1 þ vmþ1;nþ1 Þ;
etc:;
and 2 V1 ¼V pmþ12;nþ12 þ qmþ12;nþ12 pmþ12;nþ12 ; 2 2 V 2 ¼ V 0 pmþ12;nþ12 þ qmþ12;nþ12 qmþ12;nþ12 : 0
2
The integrator (2.5a)–(2.5d) preserves the discrete multisymplectic conservation law dpmþ12;nþ1 ^ dqmþ12;nþ1 dpmþ12;n ^ dqmþ12;n Dt
þ
xmþ1;nþ12 xm;nþ12 Dx
¼ 0;
ð2:6Þ
where xl;nþ12 ¼ dvl;nþ12 ^ dpl;nþ12 þ dwl;nþ12 ^ dql;nþ12 ;
l ¼ m; m þ 1:
Summing (2.6) over the space index m and using the periodic boundary condition leads to X X dpmþ12;nþ1 ^ dqmþ12;nþ1 ¼ dpmþ12;n ^ dqmþ12n ; ð2:7Þ m
m
which shows that the integrator (2.5a)–(2.5d) is also a symplectic integrator for the periodic NLS equation. In practice, we usually need to know the values of w only. Eliminating v and w from (2.5a)–(2.5d) and using w = p + iq, we can obtain an integrator for w: i
w½m ;nþ1 w½m ;n wmþ1;nþ12 2wm;nþ12 þ wm1;nþ12 1 þ þ Gm;n ¼ 0; 2 Dt Dx2
where w½m ;r ¼ 14 wm1;r þ 2wm;r þ wmþ1;r ; 2
r ¼ n; n þ 1; 2
Gm;n ¼ V 0 ðjwm12;nþ12 j Þwm12;nþ12 þ V 0 ðjwmþ12;nþ12 j Þwmþ12;nþ12 :
ð2:8Þ
1398
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
Fig. 1. Schematic of the six-point multisymplectic integrator (2.8).
The integrator (2.8) is a six-point integrator (see Fig. 1). We also call (2.8) a multisymplectic integrator because it is equivalent to (2.5a)–(2.5d) if we apply (2.5c) and (2.5d) to obtain the values of v and w respectively when necessary. The integrator (2.8) was generalized to an integrator for the Schro¨dinger equation with variable coefficients [16].
3. Accuracy and stability of the multisymplectic integrator In this section, we will explore the accuracy and stability of the multisymplectic integrator (2.8). We first discuss its accuracy. For its accuracy, we have the following theorem: Theorem 3.1. The multisymplectic integrator (2.8) for the NLS equation has accuracy of OðDt2 þ Dx2 Þ. Proof. Using Taylor series expansion, we can obtain 1 wm1;r þ 2wm;r þ wmþ1;r ¼ wm;r þ OðDx2 Þ; 4
r ¼ n; n þ 1;
1 wmþ12;nþ12 ¼ ðwm;n þ wmþ1;n þ wm;nþ1 þ wmþ1;nþ1 Þ 4 1 1 ¼ w m þ Dx; n þ Dt þ OðDt2 þ Dx2 Þ; 2 2 1 wm12;nþ12 ¼ ðwm1;n þ wm;n þ wm1;nþ1 þ wm;nþ1 Þ 4 1 1 ¼ w m Dx; n þ Dt þ OðDt2 þ Dx2 Þ: 2 2 Using the above equations ðmDx; ðn þ 12ÞDtÞ, we have
and
Taylor
series
expanded
at
point
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
w½m ;nþ1 w½m ;n wm;nþ1 wm;n ¼ þ OðDx2 Þ Dt Dt 1 ¼ wt mDx; n þ Dt þ OðDt2 Þ þ OðDx2 Þ 2 1 ¼ wt mDx; n þ Dt þ OðDt2 þ Dx2 Þ; 2 wmþ1;nþ12 2wm;nþ12 þ wm1;nþ12 Dx2
1399
ð3:1Þ
1 ¼ wxx mDx; n þ Dt þ OðDt2 þ Dx2 Þ; 2 ð3:2Þ
and 1 F wm12;nþ12 þ F wmþ12;nþ12 2 1 1 1 ¼ F w m Dx; n þ Dt 2 2 2 1 1 1 þ F w m þ Dx; n þ Dt þ OðDt2 þ Dx2 Þ 2 2 2 1 1 1 1 ¼ F w mDx; n þ Dt Gx mDx; n þ Dt Dx þ OðDx2 Þ 2 2 4 2 1 1 1 1 þ F w mDx; n þ Dt þ Gx mDx; n þ Dt Dx 2 2 4 2 1 þ OðDx2 Þ þ OðDt2 þ Dx2 Þ ¼ F w mDx; n þ Dt 2 þ OðDt2 þ Dx2 Þ;
ð3:3Þ
where G(x, t) = F(w(x, t)). Eqs. (3.1)–(3.3) show that the multisymplectic integrator (2.8) has accuracy of OðDt2 þ Dx2 Þ. This completes the proof. h Now we discuss the stability of the multisymplectic integrator (2.8). For the linear stability, we have Theorem 3.2. The multisymplectic integrator (2.8) is unconditionally linearly stable. This theorem can be easily proved by using the von Neumann method. For details, see [7]. It is relatively difficult to explore the nonlinear stability of (2.8). In the following, we will prove its nonlinear stability in a certain sense. We first prove two lemmas.
1400
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
Lemma 3.1. For the multisymplectic integrator (2.8), the following discrete conservation law holds: 2 M 2 M 1 1 X X wmþ12;nþ1 ¼ wmþ12;n ; m¼0
ð3:4Þ
m¼0
L where M ¼ DX .
Proof. We use the idea developed in [16]. First notice that the periodic boundary condition results in the following equations: M 1 X
M 1 X wm1;n wm;n ¼ wm;n wmþ1;n ;
m¼0
m¼0
M 1 X w
m;n
m¼0
1 X 2 M 2 ¼ w mþ1;n ;
etc:
m¼0
We will use the above equations in the following derivation. Multiplying (2.8) ^ by the complex conjugate w m;nþ12 of wm;nþ12 and summing over m lead to i 1 I 1 þ 2 I 2 þ I 3 ¼ 0; Dt Dx
ð3:5Þ
where I1 ¼
M 1 X ^ w½m ;nþ1 w½m ;n w m;nþ12 m¼0
¼
M 1 2 2 1X ^ w þ 2Reðw w w Þ þ m;nþ1 m;nþ1 mþ1;nþ1 mþ1;nþ1 8 m¼0
M 1 2 1X ^ w w 2 þ 2Reðw w Þ þ m;n m;n mþ1;n mþ1;n 8 m¼0
þ
M 1 i X ^ Þ þ 2Imðw ^ ^ Imðwmþ1;nþ1 w m;n m;nþ1 wm;n Þ þ Imðwm1;nþ1 wm;n Þ ; 4 m¼0
ð3:6Þ I2 ¼
M 1 X ^ wmþ1;nþ12 2wm;nþ12 þ wm1;nþ12 w m;nþ12 m¼0
¼2
M 1 X m¼0
2 ^ Reðw w Þ w ; 1 1 1 m;nþ2 mþ1;nþ2 m;nþ2
ð3:7Þ
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
1401
and M 1 1X ^ Gm;n w m;nþ12 2 m¼0 M1 2 2 2 1X 0 0 ¼ V wm12;nþ12 þ V wmþ12;nþ12 wm;nþ12 4 m¼0 M 1 2 1X 0 ^ þ V wmþ12;nþ12 Re wm;nþ12 wmþ1;nþ12 ; 2 m¼0
I3 ¼
ð3:8Þ
where Re and Im denote the real and imaginary parts of a complex respectively. Now taking imaginary part on both sides of (3.5) and using (3.6)–(3.8), we obtain M 1 X 2 2 ^ w m;nþ1 þ 2Reðwm;nþ1 wmþ1;nþ1 Þ þ wmþ1;nþ1 m¼0
¼
M 1 X
w
m;n
2 w w Þ þ m;n mþ1;n mþ1;n
2 ^ þ 2Reðw
m¼0
which is equivalent to (3.4). h The discrete conservation law (3.4) also holds for a multisymplectic integrator for the Schro¨dinger equation with variable coefficients [16]. Lemma 3.2. Let CM be the M dimensional complex vector space over C. Define a function k kmid : CM ! R; kzkmid ¼
M 1 X zm þ zmþ1 2 2 m¼0
!1=2 ;
8z ¼ ðz0 ; z1 ; . . . ; zM1 Þ 2 CM ;
where we suppose zM = z0. Then k Æ kmid is a seminorm. In particular, k Æ kmid is a norm if M is an odd number. Proof. In order to prove that k Æ kmid is a seminorm, we need to prove: ðaÞ kzkmid P 0;
8z 2 CM ;
ðbÞ kazkmid ¼ jajkzkmid ;
8a 2 C and 8z 2 CM ;
ðcÞ kz þ ~zkmid 6 kzkmid þ k~zkmid ;
8z; ~z 2 CM :
(a) and (b) hold obviously. Now prove (c) as follows:
1402
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
kz þ ~zkmid ¼ ¼ ¼
2 !1=2 zm þ ~zm þ zmþ1 þ ~zmþ1 2 m¼0 2 !1=2 M 1 X z þ z ~ z þ ~ z m mþ1 m mþ1 þ 2 2 m¼0 !1=2 M 1 X zm þ zmþ1 ~zm þ ~zmþ1 2 ; ~rm ¼ rm ¼ jrm þ ~rm j 2 2 m¼0 M 1 X
¼ kr þ ~rkE ðk kE is the Euclid normÞ 6 krkE þ k~rkE ¼ kzkmid þ k~zkmid ; where r = (r0, r1, . . . , rM1), ~r ¼ ð~r0 ; ~r1 . . . ; ~rM1 Þ. Now we suppose that M is an odd number. From kzkmid = 0, we can obtain z0 ¼ z1 ;
z1 ¼ z2 ;
...;
zM1 ¼ z0 ;
which means z0 = (1)Mz0. Since M is an odd number, we deduce that z0 = 0. Consequently, z1 = z2 = = zM1 = z0 = 0, i.e. z = 0. Therefore, k Æ kmid is a norm if M is an odd number. Let wn = (w0,n, w1,n, . . . , wM1,n) be the solution of the multisymplectic integrator (2.8), and let M be an odd number. From Lemmas 3.1 and 3.2, we have kwnþ1 kmid ¼ kwn kmid : Thus, we have proved.
ð3:9Þ h
Theorem 3.3. The multisymplectic integrator (2.8) is nonlinearly stable with respect to the norm k Æ kmid. Remark. Since we require that M be an odd number, the nonlinear stability analysis of (2.8) is incomplete. In fact, from our numerical experiments, we see that the stability of (2.8) does not depend on the oddness of the number of the spatial points. Therefore, further work (maybe difficult) remains to be done with respect to the nonlinear stability analysis of (2.8). In principle, we can employ the norm k Æ kmid to investigate the convergence of (2.8). As pointed out in the above Remark, however, we need to suppose M is an odd number. Therefore, we will take another approach to investigate the convergence of (2.8). The basic idea is to explore the relationship between (2.8) and variational integrators for (1.1) in the framework of discrete multisymplectic field theory developed in [18]. In the next section, we will present the multisymplectic geometry for the NLS equation in both continuous and discrete versions. Then we will discuss the convergence of (2.8) in Section 5.
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
1403
4. Multisymplectic geometry for the NLS equation We will work in local coordinates and use the notion of prolongation spaces instead of that of jet bundles. As was pointed out in [20], the construction of prolongation spaces is a greatly simplified version of the theory of jet bundles. We begin with the pair of real-valued equations (2.1). The covariant configuration space for (2.1) is X · U, where X = (x, t) represents the space of independent variables and U = (p, q) represents the space of dependent variables. We define the first order prolongation of X · U as U ð1Þ ¼ X U U 1 ; where U1 = (px,qx,pt,qt) represents the space consisting of first order partial derivatives. Let u : X ! U be a section and we suppose u 2 H2[X], where H2[X] is the second order Soblev space defined onX. We denote first prolongation of u by pr1 u ¼ ðx; t; p; q; px ; qx ; pt ; qt Þ: The Lagrangian density for the nonlinear Schro¨dinger equation (2.1) is Lðpr1 ðuÞÞ ¼ Lðpr1 ðuÞÞ dx ^ dt;
ð4:1Þ
where 1 Lðpr1 ðuÞ ¼ ½p2x þ q2x þ pqt qpt V ðp2 þ q2 Þ : 2
ð4:2Þ
Corresponding the lagrangian density (4.1), the action functional is defined as follows: Z Lðpr1 ðuÞÞ;
SðuÞ ¼
M is an open set of X ;
ð4:3Þ
M
where u 2 H 2 ½M . Let V be a vector filed on X · U with the form o o o o V ¼ nðx; tÞ þ sðx; tÞ þ aðx; t; p; qÞ þ bðx; t; p; qÞ : ox ot op oq The flow exp(kV) of the vector V is a one-parameter transformation group of ~ ! U, ~ :M X · U and transforms a section u : M ! U to a family of sections u which depend on the parameter k. Now we calculate the variation of the action functional (4.3) as follows: d dS ¼ Sð~ uÞ dkk¼0 Z d 1 2 ½~ p~x þ ~ q~2x þ ~ p~ q~t ~ q~ p~t V ð~ p2 þ ~q2 Þ d~x ^ d~t ¼ dk ~ 2 M k¼0 Z Z ¼ I 1 dx ^ dt þ I 2 dx þ I 3 dt; ð4:4Þ M
oM
1404
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
where I 1 ¼ ðpt þ qxx þ V 0 ðp2 þ q2 ÞqÞðqx n þ qt s bÞ þ ðqt pxx V 0 ðp2 þ q2 ÞpÞ ða px n pt sÞ; and I 2 ¼ 12½ðpqx qpx Þn þ V ðp2 þ q2 Þ p2x q2x Þs þ qa pb ; I 3 ¼ 12½ðpqt qpt p2x q2x V ðp2 þ q2 ÞÞn 2ðpx pt þ qx qt Þs þ 2px a þ 2qx b : If n, s, and a, b have compact support on M, then the boundary integral in (4.4) vanishes. In this case, from the requirement of dS ¼ 0, we can obtain the Euler–Lagrange equation pt þ qxx þ V 0 p2 þ q2 q ¼ 0; ð4:5Þ qt pxx V 0 p2 þ q2 p ¼ 0: If the condition that n, s, and a, b have compact support on M is not imposed, then from the boundary integral in (4.4) we can define the Cartan form 1 1 HL ¼ px dp ^ dt þ q dp ^ dx þ qx dq ^ dt p dq ^ dx 2 2 1 2 ðpx þ q2x þ V ðp2 þ q2 ÞÞ dx ^ dt; 2 which satisfies Z Z I 2 dx þ I 3 dt ¼ ðpr1 uÞ ðpr1 V c HL Þ: oM
ð4:6Þ
ð4:7Þ
oM
The multisymplectic form is defined to be XL ¼ dHL . Let u be the solution of (2.1), we have Theorem 4.1 [13,18]. Let gk and fk be two one-parameter symmetry groups of Eq. (2.1), and V and W are the corresponding infinitesimal symmetries, then we have the multisymplectic form formula Z pr1 ðuÞ ðpr1 ðV Þc pr1 ðW Þc XL Þ ¼ 0: ð4:8Þ oU
Now we present a discrete version of the multisymplectic geometry for (2.1) in the framework developed in [18]. We take X ¼ ðxi ; tj Þ :¼ ði; jÞ and U ¼ ðpij ; qij Þ as the discrete versions of X and U respectively. A rectangle h of X is an ordered quadruple of the form ¼ ðði; jÞ; ði þ 1; jÞ; ði þ 1; j þ 1Þ; ði; j þ 1ÞÞ:
ð4:9Þ
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
1405
The ith component of h is the ith vertex of the rectangle, denoted by hi. A point ði; jÞ 2 X is touched by a rectangle if it is a vertex of that rectangle. If M X, then (i, j) is an interior point of M if M contains all four rectangles that touch (i, j). We denote M as the union of all rectangles touching interior points of M. A boundary point of M is a point in M which is not an interior point. If M ¼ M, we call M regular. int M is the set of the interior points of M, and oM is the set of boundary points. The discrete first order prolongation of X U is defined to be Uð1Þ ð; pi;j ; qi;j ; piþ1;j ; qiþ1;j ; piþ1;jþ1 ; qiþ1;jþ1 ; pi;jþ1 ; qi;jþ1 Þ: The first order prolongation of the discrete map ud : X ! U; ud ði; jÞ ¼ ðpi;j ; qi;j Þ is defined by pr1 ud ¼ ð; pi;j ; qi;j ; piþ1;j ; qiþ1;j ; piþ1;jþ1 ; qiþ1;jþ1 ; pi;jþ1 ; qi;jþ1 Þ:
ð4:10Þ
Let M ¼ ½a; b ½c; d be a rectangular domain and consider a uniform rectangular subdivision a ¼ x0 < x1 < < xM1 < xM ¼ b;
ð4:11Þ
c ¼ t0 < t1 < < tN 1 < tN ¼ d; where xi ¼ a þ iDx; tj ¼ c þ jDt; ba d c ; Dt ¼ ; Dx ¼ M N
i ¼ 0; 1; . . . ; M;
j ¼ 0; 1; . . . ; N :
We denote the uniform rectangular subdivision (4.11) by M ¼ fðxi ; tj Þ : i ¼ 0; 1; . . . ; M; j ¼ 0; 1; . . . ; N g:
ð4:12Þ
For a discrete Lagrangian L : Uð1Þ ! R, we define the discrete action functional as follows: SðuÞ ¼
M1 X N 1 X i¼0
Lðpi;j ; qi;j ; piþ1;j ; qiþ1;j ; piþ1;jþ1 ; qiþ1;jþ1 ; pi;jþ1 ; qi;jþ1 ÞDxDt:
j¼0
ð4:13Þ Now we calculate the variation of the discrete action functional (4.13). Let ~ qi;j ¼ qi;j þ kdqi;j ; pi;j ¼ pi;j þ kdpi;j ; ~
i ¼ 0; 1; . . . ; M; j ¼ 0; 1; . . . ; N :
1406
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
By direct calculation, we can obtain
M 1 X N 1 d X Lð~ pi;j ; ~ qi;j ; ~ piþ1;j ; ~ qiþ1;j ; ~ piþ1;jþ1 ; ~ piþ1;jþ1 ; ~ pi;jþ1 ; ~ qi;jþ1 ÞDxDt dk k¼0 i¼0 j¼0
dS ¼
N 1 X N 1 X
¼
D1 Li;j þ D3 Li1;j þ D5 Li1;j1 þ D7 Li;j1 DxDtdpi;j
i¼1 j¼1 N1 X N 1 X
þ þ
þ
D2 Li;j þ D4 Li1;j þ D6 Li1;j1 þ D8 Li;j1 DxDtdqi;j
i¼1 j¼1 M1 X
D1 Li;0 dpi;0 þ D3 Li10 dpi;0 þ D5 Li1;N 1 dpi;N þ D7 Li;N 1 dpi;N DxDt
i¼1 M1 X
D2 Li;0 dqi;0 þ D4 Li1;0 dqi;0 þ D6 Li1;N 1 dqi;N þ D8 Li;N 1 dqi;N DxDt
i¼1 N1 X þ D1 L0;j dp0;j þ D3 LM1;j dpM;j þ D5 LM1;j dpM;jþ1 þ D7 L0;j dp0;jþ1 DxDt; j¼0 N1 X þ D2 L0;j dq0;j þ D4 LM1;j dqM;j þ D6 LM1;j dqM;jþ1 þ D8 L0;j dq0;jþ1 DxDt: j¼0
ð4:14Þ
Here Di denotes the partial derivative of L with respect to the ith argument and Li;j ¼ Lðpi;j ; qi;j ; piþ1;j ; qiþ1;j ; piþ1;jþ1 ; qiþ1;jþ1 ; pi;jþ1 ; qi;jþ1 Þ; Li1;j ¼ Lðpi1;j ; qi1;j ; pi;j ; qi;j ; pi;jþ1 ; qi;jþ1 ; pi1jþ1 ; qi1jþ1 Þ; Li1;0 ¼ Lðpi1;0 ; qi1;0 ; pi;0 ; qi;0 ; pi1 ; qi1 ; pi1;1 ; qi1;1 Þ;
etc:
If dpi,0 = dpi,N = dp0,j = dpM,j = 0 and dqi,0 = dqi,N = dq0,j = dqM,j = 0, i = 0, 1, . . . , M; j = 0, 1, . . . , N, then the boundary terms in (4.14) vanish. In this case, from the requirement of dS ¼ 0, we can obtain the discrete Euler– Lagrange equation D1 Li;j þ D3 Li1;j þ D5 Li1;j1 þ D7 Li;j1 ¼ 0; ð4:15Þ D2 Li;j þ D4 Li1;j þ D6 Li1;j1 þ D8 Li;j1 ¼ 0: We call the discrete Euler–Lagrange equation (4.15) a variational integrator. If we do not impose the vanishing of the boundary variations, then from the boundary terms in (4.14) we can define the following four one forms oLi;j oLi;j oLi;j oLi;j dpi;j þ dqi;j ; H2L ¼ dpiþ1;j þ dq ; opi;j oqi;j opiþ1;j oqiþ1;j iþ1;j oLi;j oLi;j dpiþ1;jþ1 þ dq ; H3L ¼ opiþ1;jþ1 oqiþ1;jþ1 iþ1;jþ1 oLi;j oLi;j dpi;jþ1 þ dq : H4L ¼ opi;jþ1 oqi;jþ1 i;jþ1 H1L ¼
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
1407
The above P four one forms are the discrete analog of the Cartan form (4.6). 4 l i;j Note that l¼1 HL ¼ dL . The solution ui,j(=ud(i, j) of the discrete Euler– Lagrange equation (4.15) satisfies the discrete multisymplectic form formula 0 1 X X @ ðpr1 ud Þ ðpr1 V1 c pr1 V2 c XlL ÞA ¼ 0; ð4:16Þ ;\oM6¼;
l;l 2oM
where XlL ¼ dHlL ; l ¼ 1; . . . ; 4 and V1 and V2 are solutions of the linearized equation of (4.15). For the proof of the discrete multisymplectic form formula (4.16), see [18].
5. Relationship between the multisymplectic and variational integrators Now we consider a concrete discretization of the Lagrangian (4.2). Let Lðpi;j ; qi;j ; piþ1;j ; qiþ1;j ; piþ1;jþ1 ; qiþ1;jþ1 ; pi;jþ1 ; qi;jþ1 Þ h i 2 2 ¼ 12 ðDx pi;j Þ þ ðDx qi;j Þ þ piþ12;jþ12 Dt qi;j qiþ12;jþ12 Dt pi;j 2 2 12V ðpiþ12;jþ12 Þ þ ðqiþ12;jþ12 Þ ;
ð5:1Þ
where Dx pi;j ¼ Dt pi;j ¼
piþ1;jþ12 pi;jþ12 Dx piþ12;jþ1 piþ12;j Dt
;
Dx qi;j ¼
;
Dt qi;j ¼
qiþ1;jþ12 qi;jþ12 Dx qiþ12;jþ1 qiþ12;j Dt
; :
For the discrete Lagrangian (5.1), the discrete Euler–Lagrange equation (4.15) takes the form p½i ;jþ1 p½i ;j1 qiþ1;½j 2qi;½j þ qi1½j 1 þ þ N 1 ¼ 0; 4 2Dt Dx2 q½i ;jþ1 q½i ;j1 piþ1;½j 2pi;½j þ pi1;½j 1 N 2 ¼ 0; 4 2Dt Dx2 where pi1;l þ 2pi;l þ piþ1;l ; 4 l ¼ j 1; j þ 1;
p½i ;l ¼
ps;j1 þ 2ps;j þ ps;jþ1 ; 4 s ¼ i 1; i; i þ 1;
ps;½j ¼
q½i ;l ¼
qi1;l þ 2qi;l þ qiþ1;l ; 4
qs;½j ¼
qs;j1 þ 2qs;j þ qs;jþ1 ; 4
ð5:2Þ
1408
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417 2
2
2
2
N 1 ¼ V 0 ððpiþ12;jþ12 Þ þ ðqiþ12;jþ12 Þ Þqiþ12;jþ12 þ V 0 ððpi12;jþ12 Þ þ ðqi12;jþ12 Þ Þqi12;jþ12 2
2
2
þ V 0 ððpiþ12;j12 Þ þ ðqiþ12;j12 Þ Þqiþ12;j12 þ V 0 ððpi12;j12 Þ þ ðqi12;j12 Þqi12;j12 ; N 2 ¼ V 0 ððpiþ12;jþ12 Þ2 þ ðqiþ12;jþ12 Þ2 Þpiþ12;jþ12 þ V 0 ððpi12;jþ12 Þ2 þ ðqi12;jþ12 Þ2 Þpi12;jþ12 2
2
2
þ V 0 ððpiþ12;j12 Þ þ ðqiþ12;j12 Þ Þpiþ12;j12 þ V 0 ððpi12;j12 Þ þ ðqi12;j12 Þpi12;j12 : Using w = p + iq, (5.2) becomes an integrator for w: i
w½m ;nþ1 w½m ;n1 wmþ1;½n 2wm;½n þ wm1;½n 1 þ þ N ¼ 0; 4 2Dt Dx2
ð5:3Þ
where w½m ;r ¼
wm1;r þ 2wm;r þ wmþ1;r ; 4
r ¼ n 1; n þ 1;
wl;n1 þ 2wl;n þ wl;nþ1 ; l ¼ m 1; m; m þ 1; 4 2 2 N ¼ V 0 jwmþ12;nþ12 j wmþ12;nþ12 þ V 0 jwm12;nþ12 j wm12;nþ12 2 2 þ V 0 jwmþ12n12 j wmþ12;n12 þ V 0 jwm12;n12 j wm12;n12 : wl;½n ¼
Now we consider the relationship between the integrators (5.3) and (2.8). In fact, the integrator (2.8) is a reduction of (5.3). To see this, we replace n by n 1 in (2.8) and obtain i
w½m ;n w½m ;n1 wmþ1;n12 2wm;n12 þ wm1;n12 1 þ þ Gm;n1 ¼ 0; 2 Dt Dx2
ð5:4Þ
Adding (5.4) to (2.8) and multiplying the resulting equation by 12 lead to (5.3). Therefore, the solution of (2.8) satisfies (5.3). Now we consider another discretization of the Lagrangian (4.2) based on the finite element method [9]. To begin with, let us go back to the variational problem of the action functional (4.3). The finite element method is an approximate method for solving the problem. Instead of solving the variational problem in the infinite dimensional function space H 2 ½M , the finite element method solves the problem in a finite dimensional subspace W ½M . W ½M consists of piecewise polynomials interpolating u 2 H 2 ½M . For details of the finite element method, we refer the reader to [2,10,11]. We consider bilinear Lagrange elements for the uniform rectangle subdivision (4.12): M ¼ fðxi ; tj Þ : i ¼ 0; 1; . . . ; M; j ¼ 0; 1; . . . ; N g. We begin by constructing the subspace W 1 ½M which consists of piecewise bilinear
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
1409
polynomials interpolating u 2 H 2 ½M . To do so, we first construct bilinear interpolation basis functions /mn(x, t) defined on M which satisfy /m;n ðxi ; tj Þ ¼ dmi dnj ;
i ¼ 0; 1; . . . ; M;
j ¼ 0; 1; . . . ; N :
ð5:5Þ
Note that /m,n(x, t) can be obtained from x-direction and t-direction linear interpolation basis functions. Therefore, we now construct x-direction and t-direction linear interpolation basis functions respectively. First consider x-direction linear interpolation basis functions /mð1Þ ðxÞ on [a, b] which satisfy /mð1Þ ðxi Þ ¼ dmi ;
i ¼ 0; 1; . . . ; M:
We list the basis functions as follows: 0 1 xx ; x0 6 x 6 x1 ; ð1Þ Dx /0 ðxÞ ¼ 0; otherwise; ð1Þ /M ðxÞ
¼
M 1 þ xx ; xM1 6 x 6 xM ; Dx 0; otherwise;
and for m = 1, 2, . . . , M 1 8 xxm > < 1 þ Dx ; xm1 6 x 6 xm ; m ; xm 6 x 6 xmþ1 ; /mð1Þ ðtÞ ¼ 1 xx Dx > : 0; otherwise:
ð5:6Þ
ð5:7Þ
ð5:8Þ
ð5:9Þ
Second consider t-direction linear interpolation basis functions /nð2Þ ðtÞ on [c, d] which satisfy /nð2Þ ðtj Þ ¼ dnj ;
j ¼ 0; 1; . . . ; N :
We list the basis functions as follows: 0 ; t0 6 t 6 t1 ; 1 tt ð2Þ Dt /0 ðtÞ ¼ 0; otherwise; ð2Þ
/N ðtÞ ¼
N 1 þ tt ; tN 1 6 t 6 tN ; Dt
0;
otherwise;
and for n = 1, 2, . . . , N 1 8 ttn > < 1 þ Dt ; tn1 6 t 6 tn ; ð2Þ n /n ðtÞ ¼ 1 tt ; tn 6 t 6 tnþ1 ; Dt > : 0; otherwise:
ð5:10Þ
ð5:11Þ
ð5:12Þ
ð5:13Þ
Using (5.7)–(5.13), we obtain ð2Þ /m;n ðx; tÞ ¼ /ð1Þ m ðxÞ/n ðtÞ;
m ¼ 0; 1; . . . ; M; n ¼ 0; 1; . . . ; N :
ð5:14Þ
1410
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
Using (5.14), we can obtain the expression of P ðx; tÞ; Qðx; tÞ 2 W 1 ½M : P ðx; tÞ ¼
M X N X
pm;n /m;n ðx; tÞ;
m¼0 n¼0
Qðx; tÞ ¼
M X N X
ð5:15Þ qm;n /m;n ðx; tÞ;
m¼0 n¼0
where pm,n p(xm,tn), qm,n q(xm,tn). In the subspace W 1 ½M , the action functional (4.3) becomes S ððP ðx; tÞ; Qðx; tÞÞÞ Z b Z d 1 2 P x þ Q2x þ PQt QP t V ðP 2 þ Q2 Þ dx ^ dt ¼ 2 a c M 1 X N 1 Z xiþ1 Z tjþ1 X 1 2 2 2 2 P þ Qx þ PQt QP t V ðP þ Q Þ dx ^ dt ¼ 2 x xi tj i¼0 j¼0 ¼
M 1 X N 1 X i¼0
Lðpi;j ; qi;j ; piþ1;j ; qiþ1;j ; piþ1;jþ1 ; qiþ1;jþ1 ; pi;jþ1 ; qi;jþ1 ÞDxDt;
j¼0
ð5:16Þ where Lðpi;j ; qi;j ; piþ1;j ; qiþ1;j ; piþ1;jþ1 ; qiþ1;jþ1 ; pi;jþ1 ; qi;jþ1 Þ Z xiþ1 Z tjþ1 1 1 2 ½P x þ Q2x þ PQt QP t V ðP 2 þ Q2 Þ dx ^ dt ¼ DxDt xi 2 tj 1 ¼ ðI 1 þ I 2 I 3 Þ; 2
ð5:17Þ
where I1 ¼
1 DxDt
Z
xiþ1
Z
xi
tjþ1
P 2x þ Q2x dx ^ dt
tj
2 t t t tj j 1 ðpiþ1;j pi;j Þ þ ðpiþ1;jþ1 pi;jþ1 Þ Dt Dt tj 2 t tj t tj þ 1 dt ðqiþ1;j qi;j Þ þ ðqiþ1;jþ1 qi;jþ1 Þ Dt Dt
1 ¼ 2 Dx Dt
Z
tjþ1
2
¼
2
ðpiþ1;j pi;j Þ þ ðpiþ1;j pi;j Þðpiþ1;jþ1 pi;jþ1 Þ þ ðpiþ1;jþ1 pi;jþ1 Þ 3Dx2 2 ðqiþ1;j qi;j Þ þ ðqiþ1;j qi;j Þðqiþ1;jþ1 qi;jþ1 Þ þ ðqiþ1;jþ1 qi;jþ1 Þ2 þ ; 3Dx2 ð5:18Þ
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
1 I2 ¼ DxDt
Z
xiþ1
Z
xi
1411
tjþ1
ðPQt QP t Þ dx ^ dt
tj
Z xiþ1 Z tjþ1 1 x xi P ½i;j 1 ðqijþ1 qij Þ DxDt xi Dx tj x x i þ ðqiþ1jþ1 qiþ1j Þ dx dt Dx Z xiþ1 Z tjþ1 1 x xi Q½i;j 1 ðpijþ1 pij Þ DxDt xi Dx tj x x i þ ðpiþ1jþ1 piþ1j Þ dx dt Dx 1 ðq qi;j Þð2pi;j þ 2pi;jþ1 þ piþ1j þ piþ1;jþ1 Þ ¼ 12Dt i;jþ1 1 ðq qiþ1;j Þðpi;j þ pi;jþ1 þ 2piþ1j þ 2piþ1;jþ1 Þ þ 12Dt iþ1;jþ1 1 ðp pi;j Þð2qi;j þ 2qi;jþ1 þ qiþ1j þ qiþ1;jþ1 Þ 12Dt i;jþ1 1 ðp piþ1;j Þðqi;j þ qi;jþ1 þ 2qiþ1j þ 2qiþ1;jþ1 Þ; 12Dt iþ1;jþ1
¼
ð5:19Þ
and I3 ¼
Z
xiþ1
Z
xiþ1
Z
xi
¼
Z
xi
tjþ1
V ðP 2 þ Q2 Þ dx ^ dt
tj tjþ1
V ððP ½i;j Þ2 þ ðQ½i;j Þ2 Þ dx ^ dt:
ð5:20Þ
tj
Here P[i,j] and Q[i,j] denote the restriction of P(x, t) and Q(x, t) on [xi, xi+1; t j, tj+1] respectively. They are given by x x x xi t tj t tj i P ½i;j ¼ 1 1 pij þ 1 piþ1j Dx Dt Dx Dt x x x xi t tj t tjþ1 i þ 1 pijþ1 þ 1þ piþ1jþ1 ; Dx Dt Dx Dt x x x xi t tj t tj i Q½i;j ¼ 1 1 qij þ 1 qiþ1j Dx Dt Dx Dt x xi t t j x xi t tj þ 1 qijþ1 þ qiþ1jþ1 : Dx Dt Dx Dt For the discrete Lagrangian (5.17), the discrete Euler–Lagrange equation (4.15) takes the form
1412
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
pðiÞ;jþ1 pðiÞ;j1 qiþ1;ðjÞ 2qi;ðjÞ þ qi1;ðjÞ 1 þ þ N 1 ¼ 0; 4 2Dt Dx2 qðiÞ;jþ1 qðiÞ;j1 piþ1;ðjÞ 2pi;ðjÞ þ pi1;ðjÞ 1 N 2 ¼ 0; 4 2Dt Dx2
ð5:21Þ
where pi1;l þ 4pi;l þ piþ1;l qi1;l þ 4qi;l þ qiþ1;l ; qðiÞ;l ¼ ; l ¼ j 1; j þ 1; 6 6 pr;j1 þ 4pr;j þ pr;jþ1 qr;j1 þ 4qr;j þ qr;jþ1 ; qr;ðjÞ ¼ ; r ¼ i 1; i; i þ 1; pr;ðjÞ ¼ 6 6 ! oI 3 ðwi;j Þ oI 3 ðwi1;j Þ oI 3 ðwi1;j1 Þ oI 3 ðwi;j1 Þ N1 ¼ þ þ þ ; opi;j opi;j opi;j opi;j ! oI 3 ðwi;j Þ oI 3 ðwi1;j Þ oI 3 ðwi1;j1 Þ oI 3 ðwi;j1 Þ N2 ¼ þ þ þ ; oqi;j oqi;j oqi;j oqi;j wi;j ¼ pi;j ; qi;j ; piþ1;j ; qiþ1;j ; piþ1;jþ1 ; qiþ1;jþ1 ; pi;jþ1 ; qi;jþ1 : pðiÞ;l ¼
Using w = p + iq, (5.21) becomes an integrator for w: i
wðmÞ;nþ1 wðmÞ;n1 wmþ1;ðnÞ 2wm;ðnÞ þ wm1;ðnÞ 1 þ þ N ¼ 0; 4 2Dt Dx2
ð5:22Þ
where wm1;r þ 4wm;r þ wmþ1;r ; r ¼ n 1; n þ 1; 6 wl;n1 þ 4wl;n þ wl;nþ1 ; l ¼ m 1; m; m þ 1; wl;ðnÞ ¼ 6 N ¼ N 2 þ iN 1 : wðmÞ;r ¼
Now we explore the relationship between the integrator (5.3) and the integrator (5.22). Using mid-point numerical integration formula in (5.18)– (5.20), we obtain piþ1;jþ12 pi;jþ12 2 qiþ1;jþ12 qi;jþ12 2 I1 þ ; Dx Dx qiþ1;jþ1 qiþ12;j piþ1;jþ1 piþ12;j qiþ12;jþ12 2 ; I 2 piþ12;jþ12 2 Dt Dt I 3 V 0 ðpiþ12;jþ12 Þ2 þ ðqiþ12;jþ12 Þ2 : Under these approximations, the discrete Lagrangian (5.17) becomes the discrete Lagrangian (5.1). Correspondingly, the variational integrator (5.22) becomes the variational integrator (5.3). This fact indicates that the multisymplectic integrator is closely related to the finite element method, which
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
1413
encourages to further explore the properties of the multisymplectic integrator by using the results from the finite element methods. In fact, we can prove the convergence of the integrator (2.8) by the convergence of the integrator (5.22) which is guaranteed by the error estimation from finite element analysis. 6. Numerical experiments For the numerical experiments, we consider the cubic NLS equation 2
iwt þ wxx þ 2jwj w ¼ 0;
ð6:1Þ
with periodic boundary condition w(0,t) = w(L,t) and the initial condition wðx; 0Þ ¼ Að1 þ 0:05 cosðlxÞÞ: ð6:2Þ pffiffiffi Here L ¼ 4 2p, l = 2p/L and A is a parameter. The initial condition (6.2) is in the vicinity of the homoclinic orbit associated with the NLS equation. The homoclinic structure causes severe sensitivities to numerical discretizations. The parameter A serves to change the complexity of the homoclinic structure. Larger A will result in higher complexity of the homoclinic structure. For details, see [15] and the references therein. Now we perform numerical experiments using the multisymplectic integrator (2.8). In our first experiment, we take A = 0.75. The space step size L Dx ¼ 127 and the time step size D t = 0.01 are used. We employ the iteration method introduced in [24] to solve (2.8). The iteration tolerance 1010 is used. We first check the discrete conservation law (3.4) numerically. Let 2 PM1 En ¼ m¼0 wmþ12;n . Then from (3.4), we have Enþ1 En ¼ 0:
ð6:3Þ
From our experiment, we see that the equation (6.3) holds exactly (in the sense of machine accuracy) in spite of the iteration error 1010 (see Fig. 2). The temporal development of waveforms is shown in Fig. 3. From the plots, we see that the spatial symmetry is well preserved. The waveforms are smooth and no instability occurs. Now we increase A to A = 1, which increases the
1
x 10
En+1 En
0.5
0
0
1000
2000
3000
4000
n
Fig. 2. The discrete conservation law (6.3).
5000
1414
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417 t=5 3
2
2
|ψ|
|ψ|
t=0 3
1
0
1
0
0
8.9
17.8
0
3
3
2
2
1
0
1
8.9
17.8
0
x 3
2
2
|ψ|
|ψ|
17.8
t=25
3
1
0
1
0
0
8.9
17.8
0
8.9
x
x
t=30
t=35
3
3
2
2
|ψ|
|ψ|
8.9
x
t=20
1
0
17.8
1
0
0
8.9
17.8
0
x
8.9
17.8
x
t=40
t=45
3
3
2
2
|ψ|
|ψ|
17.8
0
0
1
0
1
0
0
8.9
17.8
0
8.9
x
x
t=47
t=50
3
3
2
2
|ψ|
|ψ|
8.9
x t=15
|ψ|
|ψ|
x t=10
1
0
17.8
1
0
0
8.9
x
17.8
0
8.9
x
Fig. 3. The temporal development of waveforms. A = 0.75.
17.8
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417 t=5
3
3
2
2
|ψ|
|ψ|
t=0
1
0
8.9
17.8
0
x t=15
3
3
2
2
1
0
1
8.9
17.8
0
x
17.8
t=25
3
3
2
2
|ψ|
|ψ|
8.9
x
t=20
1
0
1
0
0
8.9
17.8
0
x
8.9
17.8
x
t=30
t=35
3
3
2
2
|ψ|
|ψ|
17.8
0
0
1
0
1
0
0
8.9
17.8
0
8.9
x
x
t=40
t=45
3
3
2
2
|ψ|
|ψ|
8.9
x t=10
|ψ|
|ψ|
1
0
0
1
0
17.8
1
0
0
8.9
17.8
0
x
8.9
17.8
x
t=47
t=50
3
3
2
2
|ψ|
|ψ|
1415
1
0
1
0
0
8.9
x
17.8
0
8.9
x
Fig. 4. The temporal development of waveforms. A = 1.0.
17.8
1416
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
complexity and sensitivity of the homoclinic structure. In this case, we use the L space step size Dx ¼ 256 and the time step size Dt = 0.01. Waveforms at different times are shown in Fig. 4. Again, the spatial symmetry is well preserved and no instability occurs. Notice that, in the first experiment, we use an odd number of spatial points while in the second experiment, an even number of spatial points are used. The numerical experiments indicate that the stability of the multisymplectic integrator (2.8) does not depend on the oddness of the spatial points.
7. Conclusions In this paper, we present a multisymplectic integrator which is derived from the multisymplectic formulation of the periodic nonlinear Schro¨dinger equation. The integrator satisfies a discrete multisymplectic conservation law. Its accuracy is proved. By introducing a norm, we investigate its nonlinear stability in some sense. We prove that this multisymplectic integrator is a reduction of a variational integrator which is derived from the MPS framework of discrete multisymplectic field theory. Furthermore, we prove that this integrator is also an approximation of another variational integrator which is derived from the finite element method. Numerical experiments are reported. Our results in this paper can be naturally generalized to those for higher dimensional NLS equation: 2 iwt þ wxx þ wyy þ V 0 jwj w ¼ 0; and 2 iwt þ wxx þ wyy þ wzz þ V 0 jwj w ¼ 0; with corresponding periodic boundary conditions.
Acknowledgements This work is supported by National Natural Science Foundation of China under Grant No. 40474047 and also partially by Chinese Academy of Sciences with Key Project of Knowledge Innovation KZCX1-SW-18.
References [1] M.J. Ablowitz, J.F. Ladik, A nonlinear difference scheme and inverse scattering, Stud. Appl. Math 55 (1976) 213–229.
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
1417
[2] S.C. Brenner, L.R. Scott, The Mathematical Theory of Finite Element Methods, SpringerVerlag, New York, 1994. [3] T.J. Bridges, Multi-symplectic structures and wave propagation, Math. Proc. Cam. Phil. Soc. 121 (1997) 147–190. [4] T.J. Bridges, S. Reich, Multi-symplectic integrators: numerical schemes for Hamiltonian PDEs that conserve symplecticity, Phys. Lett. A 284 (2001) 184–193. [5] J.B. Chen, Total variation in discrete multisymplectic field theory and multisymplectic energy momentum integrators, Lett. Math. Phys. 61 (2002) 63–73. [6] J.B. Chen, H.Y. Guo, K. Wu, Total variation in Hamiltonian formalism and symplecticenergy integrators, J. Math. Phys. 44 (2003) 1688–1702. [7] J.B. Chen, M.Z. Qin, Y.F. Tang, Symplectic and multisymplectic methods for the nonlinear Schro¨dinger equation, Comput. Math. Appl. 43 (2002) 1095–1106. [8] J.B. Chen, A multisymplectic variational integrator for the nonlinear Schro¨dinger equation, Numer. Meth. Part. Diff. Eq. 18 (2002) 523–536. [9] J.B. Chen, Variational integrators for discrete multisymplectic field theory based on the finite element method, preprint. [10] P.G. Ciarlet, Numerical Analysis of the Finite Element Method, Les Presses, de LUniversite´ de Montreal, Quebec, 1976. [11] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, New York, Oxford, 1978. [12] M. Delfour, M. Fortin, G. Payre, Finite-difference solutions of a nonlinear Schro¨dinger equation, J. Comput. Phys. 44 (1981) 277–288. [13] P.L. Garcia, The Poincare´–Cartan invariant in the calculus of variations, Symp. Math. 14 (1974) 219–246. [14] A. Hasegawa, Optical Solitons in Fibers, Springer-Verlag, Berlin, 1989. [15] B.M. Herbst, F. Varadi, M.J. Ablowitz, Symplectic methods for the nonlinear Schro¨dinger equation, Math. Comput. Simul. 37 (1994) 353–369. [16] J.L. Hong, Y. Liu, H. Munthe-Kass, A. Zanna, On a multisymplectic scheme for Schro¨dinger equations with variable coefficients, preprint. [17] G.L. Lamb, Elements of Soliton Theory, John Wiley & Sons, New York, 1980. [18] J.E. Marsden, G.W. Patrick, S. Shkoller, Multisymplectic geometry, variational integrators, and nonlinear PDEs, Comm. Math. Phys. 199 (1998) 351–395. [19] R. McLachlan, Symplectic integration of Hamiltonian wave equations, Numer. Math. 66 (1994) 465–492. [20] P.J. Olver, Applications of Lie Groups to Differential Equations, second ed., Springer-Verlag, New York, 1993. [21] S. Reich, Multi-symplectic Runge–Kutta collocation methods for Hamiltonian wave equations, J. Comput. Phys. 157 (2000) 473–499. [22] C.M. Schober, Symplectic integrators for Ablowitz–Ladik discrete nonlinear Schro¨dinger equation, Phys. Lett. A 259 (1999) 140–151. [23] Y.F. Tang, L. Va´zquez, F. Zhang, V.M. Pe´rez-Garcı´a, Symplectic methods for the nonlinear Schro¨dinger equation, Comput. Math. Appl. 32 (1996) 73–83. [24] S.B. Wineberg, J.F. McGrath, E.F. Gabl, L.R. Scott, C.E. Southwell, Implicit spectral methods for wave propagation problems, J. Comput. Phys. 97 (1991) 311–336.
www.elsevier.com/locate/amc
A multisymplectic integrator for the periodic nonlinear Schro¨dinger equation Jing-Bo Chen Institute of Geology and Geophysics, Chinese Academy of Sciences, P.O. Box 9825, Beijing 100029, PR China
Abstract A multisymplectic integrator for the periodic nonlinear Schro¨dinger equation is presented in this paper. Its accuracy is proved. By introducing a norm, we investigate its nonlinear stability. We also discuss the relationship between this multisymplectic integrator and two variational integrators which are derived by using the discrete multisymplectic field theory and the finite element method. 2005 Elsevier Inc. All rights reserved. Keywords: Multisymplectic integrator; Nonlinear Schro¨dinger equation
1. Introduction We consider the nonlinear Schro¨dinger equation (NLS) written in the form 2 iwt þ wxx þ V 0 jwj w ¼ 0; ð1:1Þ
E-mail address: [email protected] 0096-3003/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.01.031
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
1395
with the periodic boundary condition w(x + L,t) = w(x, t). Here L is the period and V : R ! R is a differentiable function. We will consider the space domain [0, L] in this paper. Eq. (1.1) is one of the most important completely integrable models in the theory of solitons. Its application can be found in many areas of physics, including nonlinear optics and plasma physics [14,17]. Various numerical methods for (1.1) have been developed. See for example [1,7,8,12,15,19,21–23]. Among these numerical methods, one kind of methods are so-called symplectic methods, which take the symplectic structure of the NLS equation into account when discretizing. The basic idea is to find a finite dimensional spatial truncation using finite difference , finite element or spectral methods, so that the resulting semidiscrete system can be cast into a finite dimensional Hamiltonian system. Then we integrate the finite dimensional Hamiltonian system in time using symplectic integrators [15,19,22,23]. The symplectic methods preserve a symplectic form which is summed over the discrete space domain. Recently, Bridges and Reich introduced the concept of multisymplectic integrators for Hamiltonian partial differential equations [4]. Hamiltonian partial differential equations can be formulated as multisymplectic Hamiltonian systems which possess a local multisymplectic conservation law. The multisymplectic integrators preserve a discrete multisymplectic conservation law. We know that the symplectic integrators for partial differential equations preserve a global symplectic form which is summed over the discrete spatial domain, i.e., the local variations in symplecticity over space domain are averaged. Therefore, the symplectic integrators do not reflect the spatial distribution of symplecticity. The multisymplectic integrators overcome the limitations of symplectic integrators since they preserve a local multisymplectic conservation law. For multisymplectic integrators, see also [5–8,21]. The NLS equation can be reformulated as a multisymplectic Hamiltonian system [3]. Reich developed multisymplectic integrators for the NLS equation in [21]. Numerical experiments were performed with a multisymplectic integrator for the NLS equation in [7], which demonstrate the good performance and local conservation properties of the multisymplectic integrator. However, the accuracy and nonlinear stability of the multisymplectic integrator employed in [7] were not presented. The purpose of this paper is to investigate its accuracy and stability as well as its relationship with two variational integrators. This paper is organized as follows. In Section 2, we present a multisymplectic integrator for the periodic NLS equation. In Section 3, we investigate its accuracy and stability. In Section 4, we present the multisymplectic geometry for the NLS equation. We explore the relationship between the multisymplectic integrator and two variational integrators in Section 5. Section 6 is devoted to reporting the numerical experiments. We end this paper with some conclusions in Section 7.
1396
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
2. A multisymplectic integrator for the periodic NLS equation In this section, we will present a multisymplectic integrator for the periodic NLS equation. We begin by introducing the multisymplectic formulation of the NLS equation. Using w = p + iq, we can rewrite (1.1) as a pair of real-valued equations pt þ qxx þ V 0 p2 þ q2 q ¼ 0; ð2:1Þ qt pxx V 0 p2 þ q2 p ¼ 0: Introducing a pair of conjugate momenta v = px, w = qx, the pair of equations (2.1) can be reformulated as the following first-order system qt vx ¼ V 0 ðp2 þ q2 Þp; pt wx ¼ V 0 ðp2 þ q2 Þq;
ð2:2Þ
px ¼ v; qx ¼ w; which can be written as a multisymplectic Hamiltonian system Mzt þ Kzx ¼ rSðzÞ;
ð2:3Þ
where SðzÞ ¼ 12 v2 þ w2 þ V ðp2 þ q2 Þ ;
z ¼ ðp; q; v; wÞT ; and 0
0
B 1 B M ¼B @ 0 0
0
1
1
0 0
0 0
0 0C C C; 0 0A
0
0 0
0
0
1
B0 B K¼B @1
0 0
0 0
0
1
0
0
1
1 C C C: 0 A 0
By direct calculations, we can prove that the system (2.3) satisfies the multisymplectic conservation law o o ðdp ^ dqÞ þ ðdv ^ dp þ dw ^ dqÞ ¼ 0: ot ox
ð2:4Þ
Now we discretize (2.2) using mid-point scheme in both time and space directions and obtain qmþ12;nþ1 qmþ12;n Dt
pmþ12;nþ1 pmþ12;n Dt
vmþ1;nþ12 vm;nþ12 Dx
¼ V 1;
wmþ1;nþ12 wm;nþ12 Dx
¼ V 2;
ð2:5aÞ ð2:5bÞ
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
pmþ1;nþ12 pm;nþ12 Dx qmþ1;nþ12 qm;nþ12 Dx
1397
¼ vmþ12;nþ12 ;
ð2:5cÞ
¼ wnþ12;mþ12 ;
ð2:5dÞ
where Dx and Dt are the space step size and time step size respectively and qmþ12;nþ1 ¼ 12ðqm;nþ1 þ qmþ1;nþ1 Þ;
pm;n pðmDx; nDtÞ;
vmþ12;nþ12 ¼ 14ðvm;n þ vmþ1;n þ vm;nþ1 þ vmþ1;nþ1 Þ;
etc:;
and 2 V1 ¼V pmþ12;nþ12 þ qmþ12;nþ12 pmþ12;nþ12 ; 2 2 V 2 ¼ V 0 pmþ12;nþ12 þ qmþ12;nþ12 qmþ12;nþ12 : 0
2
The integrator (2.5a)–(2.5d) preserves the discrete multisymplectic conservation law dpmþ12;nþ1 ^ dqmþ12;nþ1 dpmþ12;n ^ dqmþ12;n Dt
þ
xmþ1;nþ12 xm;nþ12 Dx
¼ 0;
ð2:6Þ
where xl;nþ12 ¼ dvl;nþ12 ^ dpl;nþ12 þ dwl;nþ12 ^ dql;nþ12 ;
l ¼ m; m þ 1:
Summing (2.6) over the space index m and using the periodic boundary condition leads to X X dpmþ12;nþ1 ^ dqmþ12;nþ1 ¼ dpmþ12;n ^ dqmþ12n ; ð2:7Þ m
m
which shows that the integrator (2.5a)–(2.5d) is also a symplectic integrator for the periodic NLS equation. In practice, we usually need to know the values of w only. Eliminating v and w from (2.5a)–(2.5d) and using w = p + iq, we can obtain an integrator for w: i
w½m ;nþ1 w½m ;n wmþ1;nþ12 2wm;nþ12 þ wm1;nþ12 1 þ þ Gm;n ¼ 0; 2 Dt Dx2
where w½m ;r ¼ 14 wm1;r þ 2wm;r þ wmþ1;r ; 2
r ¼ n; n þ 1; 2
Gm;n ¼ V 0 ðjwm12;nþ12 j Þwm12;nþ12 þ V 0 ðjwmþ12;nþ12 j Þwmþ12;nþ12 :
ð2:8Þ
1398
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
Fig. 1. Schematic of the six-point multisymplectic integrator (2.8).
The integrator (2.8) is a six-point integrator (see Fig. 1). We also call (2.8) a multisymplectic integrator because it is equivalent to (2.5a)–(2.5d) if we apply (2.5c) and (2.5d) to obtain the values of v and w respectively when necessary. The integrator (2.8) was generalized to an integrator for the Schro¨dinger equation with variable coefficients [16].
3. Accuracy and stability of the multisymplectic integrator In this section, we will explore the accuracy and stability of the multisymplectic integrator (2.8). We first discuss its accuracy. For its accuracy, we have the following theorem: Theorem 3.1. The multisymplectic integrator (2.8) for the NLS equation has accuracy of OðDt2 þ Dx2 Þ. Proof. Using Taylor series expansion, we can obtain 1 wm1;r þ 2wm;r þ wmþ1;r ¼ wm;r þ OðDx2 Þ; 4
r ¼ n; n þ 1;
1 wmþ12;nþ12 ¼ ðwm;n þ wmþ1;n þ wm;nþ1 þ wmþ1;nþ1 Þ 4 1 1 ¼ w m þ Dx; n þ Dt þ OðDt2 þ Dx2 Þ; 2 2 1 wm12;nþ12 ¼ ðwm1;n þ wm;n þ wm1;nþ1 þ wm;nþ1 Þ 4 1 1 ¼ w m Dx; n þ Dt þ OðDt2 þ Dx2 Þ: 2 2 Using the above equations ðmDx; ðn þ 12ÞDtÞ, we have
and
Taylor
series
expanded
at
point
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
w½m ;nþ1 w½m ;n wm;nþ1 wm;n ¼ þ OðDx2 Þ Dt Dt 1 ¼ wt mDx; n þ Dt þ OðDt2 Þ þ OðDx2 Þ 2 1 ¼ wt mDx; n þ Dt þ OðDt2 þ Dx2 Þ; 2 wmþ1;nþ12 2wm;nþ12 þ wm1;nþ12 Dx2
1399
ð3:1Þ
1 ¼ wxx mDx; n þ Dt þ OðDt2 þ Dx2 Þ; 2 ð3:2Þ
and 1 F wm12;nþ12 þ F wmþ12;nþ12 2 1 1 1 ¼ F w m Dx; n þ Dt 2 2 2 1 1 1 þ F w m þ Dx; n þ Dt þ OðDt2 þ Dx2 Þ 2 2 2 1 1 1 1 ¼ F w mDx; n þ Dt Gx mDx; n þ Dt Dx þ OðDx2 Þ 2 2 4 2 1 1 1 1 þ F w mDx; n þ Dt þ Gx mDx; n þ Dt Dx 2 2 4 2 1 þ OðDx2 Þ þ OðDt2 þ Dx2 Þ ¼ F w mDx; n þ Dt 2 þ OðDt2 þ Dx2 Þ;
ð3:3Þ
where G(x, t) = F(w(x, t)). Eqs. (3.1)–(3.3) show that the multisymplectic integrator (2.8) has accuracy of OðDt2 þ Dx2 Þ. This completes the proof. h Now we discuss the stability of the multisymplectic integrator (2.8). For the linear stability, we have Theorem 3.2. The multisymplectic integrator (2.8) is unconditionally linearly stable. This theorem can be easily proved by using the von Neumann method. For details, see [7]. It is relatively difficult to explore the nonlinear stability of (2.8). In the following, we will prove its nonlinear stability in a certain sense. We first prove two lemmas.
1400
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
Lemma 3.1. For the multisymplectic integrator (2.8), the following discrete conservation law holds: 2 M 2 M 1 1 X X wmþ12;nþ1 ¼ wmþ12;n ; m¼0
ð3:4Þ
m¼0
L where M ¼ DX .
Proof. We use the idea developed in [16]. First notice that the periodic boundary condition results in the following equations: M 1 X
M 1 X wm1;n wm;n ¼ wm;n wmþ1;n ;
m¼0
m¼0
M 1 X w
m;n
m¼0
1 X 2 M 2 ¼ w mþ1;n ;
etc:
m¼0
We will use the above equations in the following derivation. Multiplying (2.8) ^ by the complex conjugate w m;nþ12 of wm;nþ12 and summing over m lead to i 1 I 1 þ 2 I 2 þ I 3 ¼ 0; Dt Dx
ð3:5Þ
where I1 ¼
M 1 X ^ w½m ;nþ1 w½m ;n w m;nþ12 m¼0
¼
M 1 2 2 1X ^ w þ 2Reðw w w Þ þ m;nþ1 m;nþ1 mþ1;nþ1 mþ1;nþ1 8 m¼0
M 1 2 1X ^ w w 2 þ 2Reðw w Þ þ m;n m;n mþ1;n mþ1;n 8 m¼0
þ
M 1 i X ^ Þ þ 2Imðw ^ ^ Imðwmþ1;nþ1 w m;n m;nþ1 wm;n Þ þ Imðwm1;nþ1 wm;n Þ ; 4 m¼0
ð3:6Þ I2 ¼
M 1 X ^ wmþ1;nþ12 2wm;nþ12 þ wm1;nþ12 w m;nþ12 m¼0
¼2
M 1 X m¼0
2 ^ Reðw w Þ w ; 1 1 1 m;nþ2 mþ1;nþ2 m;nþ2
ð3:7Þ
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
1401
and M 1 1X ^ Gm;n w m;nþ12 2 m¼0 M1 2 2 2 1X 0 0 ¼ V wm12;nþ12 þ V wmþ12;nþ12 wm;nþ12 4 m¼0 M 1 2 1X 0 ^ þ V wmþ12;nþ12 Re wm;nþ12 wmþ1;nþ12 ; 2 m¼0
I3 ¼
ð3:8Þ
where Re and Im denote the real and imaginary parts of a complex respectively. Now taking imaginary part on both sides of (3.5) and using (3.6)–(3.8), we obtain M 1 X 2 2 ^ w m;nþ1 þ 2Reðwm;nþ1 wmþ1;nþ1 Þ þ wmþ1;nþ1 m¼0
¼
M 1 X
w
m;n
2 w w Þ þ m;n mþ1;n mþ1;n
2 ^ þ 2Reðw
m¼0
which is equivalent to (3.4). h The discrete conservation law (3.4) also holds for a multisymplectic integrator for the Schro¨dinger equation with variable coefficients [16]. Lemma 3.2. Let CM be the M dimensional complex vector space over C. Define a function k kmid : CM ! R; kzkmid ¼
M 1 X zm þ zmþ1 2 2 m¼0
!1=2 ;
8z ¼ ðz0 ; z1 ; . . . ; zM1 Þ 2 CM ;
where we suppose zM = z0. Then k Æ kmid is a seminorm. In particular, k Æ kmid is a norm if M is an odd number. Proof. In order to prove that k Æ kmid is a seminorm, we need to prove: ðaÞ kzkmid P 0;
8z 2 CM ;
ðbÞ kazkmid ¼ jajkzkmid ;
8a 2 C and 8z 2 CM ;
ðcÞ kz þ ~zkmid 6 kzkmid þ k~zkmid ;
8z; ~z 2 CM :
(a) and (b) hold obviously. Now prove (c) as follows:
1402
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
kz þ ~zkmid ¼ ¼ ¼
2 !1=2 zm þ ~zm þ zmþ1 þ ~zmþ1 2 m¼0 2 !1=2 M 1 X z þ z ~ z þ ~ z m mþ1 m mþ1 þ 2 2 m¼0 !1=2 M 1 X zm þ zmþ1 ~zm þ ~zmþ1 2 ; ~rm ¼ rm ¼ jrm þ ~rm j 2 2 m¼0 M 1 X
¼ kr þ ~rkE ðk kE is the Euclid normÞ 6 krkE þ k~rkE ¼ kzkmid þ k~zkmid ; where r = (r0, r1, . . . , rM1), ~r ¼ ð~r0 ; ~r1 . . . ; ~rM1 Þ. Now we suppose that M is an odd number. From kzkmid = 0, we can obtain z0 ¼ z1 ;
z1 ¼ z2 ;
...;
zM1 ¼ z0 ;
which means z0 = (1)Mz0. Since M is an odd number, we deduce that z0 = 0. Consequently, z1 = z2 = = zM1 = z0 = 0, i.e. z = 0. Therefore, k Æ kmid is a norm if M is an odd number. Let wn = (w0,n, w1,n, . . . , wM1,n) be the solution of the multisymplectic integrator (2.8), and let M be an odd number. From Lemmas 3.1 and 3.2, we have kwnþ1 kmid ¼ kwn kmid : Thus, we have proved.
ð3:9Þ h
Theorem 3.3. The multisymplectic integrator (2.8) is nonlinearly stable with respect to the norm k Æ kmid. Remark. Since we require that M be an odd number, the nonlinear stability analysis of (2.8) is incomplete. In fact, from our numerical experiments, we see that the stability of (2.8) does not depend on the oddness of the number of the spatial points. Therefore, further work (maybe difficult) remains to be done with respect to the nonlinear stability analysis of (2.8). In principle, we can employ the norm k Æ kmid to investigate the convergence of (2.8). As pointed out in the above Remark, however, we need to suppose M is an odd number. Therefore, we will take another approach to investigate the convergence of (2.8). The basic idea is to explore the relationship between (2.8) and variational integrators for (1.1) in the framework of discrete multisymplectic field theory developed in [18]. In the next section, we will present the multisymplectic geometry for the NLS equation in both continuous and discrete versions. Then we will discuss the convergence of (2.8) in Section 5.
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
1403
4. Multisymplectic geometry for the NLS equation We will work in local coordinates and use the notion of prolongation spaces instead of that of jet bundles. As was pointed out in [20], the construction of prolongation spaces is a greatly simplified version of the theory of jet bundles. We begin with the pair of real-valued equations (2.1). The covariant configuration space for (2.1) is X · U, where X = (x, t) represents the space of independent variables and U = (p, q) represents the space of dependent variables. We define the first order prolongation of X · U as U ð1Þ ¼ X U U 1 ; where U1 = (px,qx,pt,qt) represents the space consisting of first order partial derivatives. Let u : X ! U be a section and we suppose u 2 H2[X], where H2[X] is the second order Soblev space defined onX. We denote first prolongation of u by pr1 u ¼ ðx; t; p; q; px ; qx ; pt ; qt Þ: The Lagrangian density for the nonlinear Schro¨dinger equation (2.1) is Lðpr1 ðuÞÞ ¼ Lðpr1 ðuÞÞ dx ^ dt;
ð4:1Þ
where 1 Lðpr1 ðuÞ ¼ ½p2x þ q2x þ pqt qpt V ðp2 þ q2 Þ : 2
ð4:2Þ
Corresponding the lagrangian density (4.1), the action functional is defined as follows: Z Lðpr1 ðuÞÞ;
SðuÞ ¼
M is an open set of X ;
ð4:3Þ
M
where u 2 H 2 ½M . Let V be a vector filed on X · U with the form o o o o V ¼ nðx; tÞ þ sðx; tÞ þ aðx; t; p; qÞ þ bðx; t; p; qÞ : ox ot op oq The flow exp(kV) of the vector V is a one-parameter transformation group of ~ ! U, ~ :M X · U and transforms a section u : M ! U to a family of sections u which depend on the parameter k. Now we calculate the variation of the action functional (4.3) as follows: d dS ¼ Sð~ uÞ dkk¼0 Z d 1 2 ½~ p~x þ ~ q~2x þ ~ p~ q~t ~ q~ p~t V ð~ p2 þ ~q2 Þ d~x ^ d~t ¼ dk ~ 2 M k¼0 Z Z ¼ I 1 dx ^ dt þ I 2 dx þ I 3 dt; ð4:4Þ M
oM
1404
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
where I 1 ¼ ðpt þ qxx þ V 0 ðp2 þ q2 ÞqÞðqx n þ qt s bÞ þ ðqt pxx V 0 ðp2 þ q2 ÞpÞ ða px n pt sÞ; and I 2 ¼ 12½ðpqx qpx Þn þ V ðp2 þ q2 Þ p2x q2x Þs þ qa pb ; I 3 ¼ 12½ðpqt qpt p2x q2x V ðp2 þ q2 ÞÞn 2ðpx pt þ qx qt Þs þ 2px a þ 2qx b : If n, s, and a, b have compact support on M, then the boundary integral in (4.4) vanishes. In this case, from the requirement of dS ¼ 0, we can obtain the Euler–Lagrange equation pt þ qxx þ V 0 p2 þ q2 q ¼ 0; ð4:5Þ qt pxx V 0 p2 þ q2 p ¼ 0: If the condition that n, s, and a, b have compact support on M is not imposed, then from the boundary integral in (4.4) we can define the Cartan form 1 1 HL ¼ px dp ^ dt þ q dp ^ dx þ qx dq ^ dt p dq ^ dx 2 2 1 2 ðpx þ q2x þ V ðp2 þ q2 ÞÞ dx ^ dt; 2 which satisfies Z Z I 2 dx þ I 3 dt ¼ ðpr1 uÞ ðpr1 V c HL Þ: oM
ð4:6Þ
ð4:7Þ
oM
The multisymplectic form is defined to be XL ¼ dHL . Let u be the solution of (2.1), we have Theorem 4.1 [13,18]. Let gk and fk be two one-parameter symmetry groups of Eq. (2.1), and V and W are the corresponding infinitesimal symmetries, then we have the multisymplectic form formula Z pr1 ðuÞ ðpr1 ðV Þc pr1 ðW Þc XL Þ ¼ 0: ð4:8Þ oU
Now we present a discrete version of the multisymplectic geometry for (2.1) in the framework developed in [18]. We take X ¼ ðxi ; tj Þ :¼ ði; jÞ and U ¼ ðpij ; qij Þ as the discrete versions of X and U respectively. A rectangle h of X is an ordered quadruple of the form ¼ ðði; jÞ; ði þ 1; jÞ; ði þ 1; j þ 1Þ; ði; j þ 1ÞÞ:
ð4:9Þ
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
1405
The ith component of h is the ith vertex of the rectangle, denoted by hi. A point ði; jÞ 2 X is touched by a rectangle if it is a vertex of that rectangle. If M X, then (i, j) is an interior point of M if M contains all four rectangles that touch (i, j). We denote M as the union of all rectangles touching interior points of M. A boundary point of M is a point in M which is not an interior point. If M ¼ M, we call M regular. int M is the set of the interior points of M, and oM is the set of boundary points. The discrete first order prolongation of X U is defined to be Uð1Þ ð; pi;j ; qi;j ; piþ1;j ; qiþ1;j ; piþ1;jþ1 ; qiþ1;jþ1 ; pi;jþ1 ; qi;jþ1 Þ: The first order prolongation of the discrete map ud : X ! U; ud ði; jÞ ¼ ðpi;j ; qi;j Þ is defined by pr1 ud ¼ ð; pi;j ; qi;j ; piþ1;j ; qiþ1;j ; piþ1;jþ1 ; qiþ1;jþ1 ; pi;jþ1 ; qi;jþ1 Þ:
ð4:10Þ
Let M ¼ ½a; b ½c; d be a rectangular domain and consider a uniform rectangular subdivision a ¼ x0 < x1 < < xM1 < xM ¼ b;
ð4:11Þ
c ¼ t0 < t1 < < tN 1 < tN ¼ d; where xi ¼ a þ iDx; tj ¼ c þ jDt; ba d c ; Dt ¼ ; Dx ¼ M N
i ¼ 0; 1; . . . ; M;
j ¼ 0; 1; . . . ; N :
We denote the uniform rectangular subdivision (4.11) by M ¼ fðxi ; tj Þ : i ¼ 0; 1; . . . ; M; j ¼ 0; 1; . . . ; N g:
ð4:12Þ
For a discrete Lagrangian L : Uð1Þ ! R, we define the discrete action functional as follows: SðuÞ ¼
M1 X N 1 X i¼0
Lðpi;j ; qi;j ; piþ1;j ; qiþ1;j ; piþ1;jþ1 ; qiþ1;jþ1 ; pi;jþ1 ; qi;jþ1 ÞDxDt:
j¼0
ð4:13Þ Now we calculate the variation of the discrete action functional (4.13). Let ~ qi;j ¼ qi;j þ kdqi;j ; pi;j ¼ pi;j þ kdpi;j ; ~
i ¼ 0; 1; . . . ; M; j ¼ 0; 1; . . . ; N :
1406
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
By direct calculation, we can obtain
M 1 X N 1 d X Lð~ pi;j ; ~ qi;j ; ~ piþ1;j ; ~ qiþ1;j ; ~ piþ1;jþ1 ; ~ piþ1;jþ1 ; ~ pi;jþ1 ; ~ qi;jþ1 ÞDxDt dk k¼0 i¼0 j¼0
dS ¼
N 1 X N 1 X
¼
D1 Li;j þ D3 Li1;j þ D5 Li1;j1 þ D7 Li;j1 DxDtdpi;j
i¼1 j¼1 N1 X N 1 X
þ þ
þ
D2 Li;j þ D4 Li1;j þ D6 Li1;j1 þ D8 Li;j1 DxDtdqi;j
i¼1 j¼1 M1 X
D1 Li;0 dpi;0 þ D3 Li10 dpi;0 þ D5 Li1;N 1 dpi;N þ D7 Li;N 1 dpi;N DxDt
i¼1 M1 X
D2 Li;0 dqi;0 þ D4 Li1;0 dqi;0 þ D6 Li1;N 1 dqi;N þ D8 Li;N 1 dqi;N DxDt
i¼1 N1 X þ D1 L0;j dp0;j þ D3 LM1;j dpM;j þ D5 LM1;j dpM;jþ1 þ D7 L0;j dp0;jþ1 DxDt; j¼0 N1 X þ D2 L0;j dq0;j þ D4 LM1;j dqM;j þ D6 LM1;j dqM;jþ1 þ D8 L0;j dq0;jþ1 DxDt: j¼0
ð4:14Þ
Here Di denotes the partial derivative of L with respect to the ith argument and Li;j ¼ Lðpi;j ; qi;j ; piþ1;j ; qiþ1;j ; piþ1;jþ1 ; qiþ1;jþ1 ; pi;jþ1 ; qi;jþ1 Þ; Li1;j ¼ Lðpi1;j ; qi1;j ; pi;j ; qi;j ; pi;jþ1 ; qi;jþ1 ; pi1jþ1 ; qi1jþ1 Þ; Li1;0 ¼ Lðpi1;0 ; qi1;0 ; pi;0 ; qi;0 ; pi1 ; qi1 ; pi1;1 ; qi1;1 Þ;
etc:
If dpi,0 = dpi,N = dp0,j = dpM,j = 0 and dqi,0 = dqi,N = dq0,j = dqM,j = 0, i = 0, 1, . . . , M; j = 0, 1, . . . , N, then the boundary terms in (4.14) vanish. In this case, from the requirement of dS ¼ 0, we can obtain the discrete Euler– Lagrange equation D1 Li;j þ D3 Li1;j þ D5 Li1;j1 þ D7 Li;j1 ¼ 0; ð4:15Þ D2 Li;j þ D4 Li1;j þ D6 Li1;j1 þ D8 Li;j1 ¼ 0: We call the discrete Euler–Lagrange equation (4.15) a variational integrator. If we do not impose the vanishing of the boundary variations, then from the boundary terms in (4.14) we can define the following four one forms oLi;j oLi;j oLi;j oLi;j dpi;j þ dqi;j ; H2L ¼ dpiþ1;j þ dq ; opi;j oqi;j opiþ1;j oqiþ1;j iþ1;j oLi;j oLi;j dpiþ1;jþ1 þ dq ; H3L ¼ opiþ1;jþ1 oqiþ1;jþ1 iþ1;jþ1 oLi;j oLi;j dpi;jþ1 þ dq : H4L ¼ opi;jþ1 oqi;jþ1 i;jþ1 H1L ¼
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
1407
The above P four one forms are the discrete analog of the Cartan form (4.6). 4 l i;j Note that l¼1 HL ¼ dL . The solution ui,j(=ud(i, j) of the discrete Euler– Lagrange equation (4.15) satisfies the discrete multisymplectic form formula 0 1 X X @ ðpr1 ud Þ ðpr1 V1 c pr1 V2 c XlL ÞA ¼ 0; ð4:16Þ ;\oM6¼;
l;l 2oM
where XlL ¼ dHlL ; l ¼ 1; . . . ; 4 and V1 and V2 are solutions of the linearized equation of (4.15). For the proof of the discrete multisymplectic form formula (4.16), see [18].
5. Relationship between the multisymplectic and variational integrators Now we consider a concrete discretization of the Lagrangian (4.2). Let Lðpi;j ; qi;j ; piþ1;j ; qiþ1;j ; piþ1;jþ1 ; qiþ1;jþ1 ; pi;jþ1 ; qi;jþ1 Þ h i 2 2 ¼ 12 ðDx pi;j Þ þ ðDx qi;j Þ þ piþ12;jþ12 Dt qi;j qiþ12;jþ12 Dt pi;j 2 2 12V ðpiþ12;jþ12 Þ þ ðqiþ12;jþ12 Þ ;
ð5:1Þ
where Dx pi;j ¼ Dt pi;j ¼
piþ1;jþ12 pi;jþ12 Dx piþ12;jþ1 piþ12;j Dt
;
Dx qi;j ¼
;
Dt qi;j ¼
qiþ1;jþ12 qi;jþ12 Dx qiþ12;jþ1 qiþ12;j Dt
; :
For the discrete Lagrangian (5.1), the discrete Euler–Lagrange equation (4.15) takes the form p½i ;jþ1 p½i ;j1 qiþ1;½j 2qi;½j þ qi1½j 1 þ þ N 1 ¼ 0; 4 2Dt Dx2 q½i ;jþ1 q½i ;j1 piþ1;½j 2pi;½j þ pi1;½j 1 N 2 ¼ 0; 4 2Dt Dx2 where pi1;l þ 2pi;l þ piþ1;l ; 4 l ¼ j 1; j þ 1;
p½i ;l ¼
ps;j1 þ 2ps;j þ ps;jþ1 ; 4 s ¼ i 1; i; i þ 1;
ps;½j ¼
q½i ;l ¼
qi1;l þ 2qi;l þ qiþ1;l ; 4
qs;½j ¼
qs;j1 þ 2qs;j þ qs;jþ1 ; 4
ð5:2Þ
1408
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417 2
2
2
2
N 1 ¼ V 0 ððpiþ12;jþ12 Þ þ ðqiþ12;jþ12 Þ Þqiþ12;jþ12 þ V 0 ððpi12;jþ12 Þ þ ðqi12;jþ12 Þ Þqi12;jþ12 2
2
2
þ V 0 ððpiþ12;j12 Þ þ ðqiþ12;j12 Þ Þqiþ12;j12 þ V 0 ððpi12;j12 Þ þ ðqi12;j12 Þqi12;j12 ; N 2 ¼ V 0 ððpiþ12;jþ12 Þ2 þ ðqiþ12;jþ12 Þ2 Þpiþ12;jþ12 þ V 0 ððpi12;jþ12 Þ2 þ ðqi12;jþ12 Þ2 Þpi12;jþ12 2
2
2
þ V 0 ððpiþ12;j12 Þ þ ðqiþ12;j12 Þ Þpiþ12;j12 þ V 0 ððpi12;j12 Þ þ ðqi12;j12 Þpi12;j12 : Using w = p + iq, (5.2) becomes an integrator for w: i
w½m ;nþ1 w½m ;n1 wmþ1;½n 2wm;½n þ wm1;½n 1 þ þ N ¼ 0; 4 2Dt Dx2
ð5:3Þ
where w½m ;r ¼
wm1;r þ 2wm;r þ wmþ1;r ; 4
r ¼ n 1; n þ 1;
wl;n1 þ 2wl;n þ wl;nþ1 ; l ¼ m 1; m; m þ 1; 4 2 2 N ¼ V 0 jwmþ12;nþ12 j wmþ12;nþ12 þ V 0 jwm12;nþ12 j wm12;nþ12 2 2 þ V 0 jwmþ12n12 j wmþ12;n12 þ V 0 jwm12;n12 j wm12;n12 : wl;½n ¼
Now we consider the relationship between the integrators (5.3) and (2.8). In fact, the integrator (2.8) is a reduction of (5.3). To see this, we replace n by n 1 in (2.8) and obtain i
w½m ;n w½m ;n1 wmþ1;n12 2wm;n12 þ wm1;n12 1 þ þ Gm;n1 ¼ 0; 2 Dt Dx2
ð5:4Þ
Adding (5.4) to (2.8) and multiplying the resulting equation by 12 lead to (5.3). Therefore, the solution of (2.8) satisfies (5.3). Now we consider another discretization of the Lagrangian (4.2) based on the finite element method [9]. To begin with, let us go back to the variational problem of the action functional (4.3). The finite element method is an approximate method for solving the problem. Instead of solving the variational problem in the infinite dimensional function space H 2 ½M , the finite element method solves the problem in a finite dimensional subspace W ½M . W ½M consists of piecewise polynomials interpolating u 2 H 2 ½M . For details of the finite element method, we refer the reader to [2,10,11]. We consider bilinear Lagrange elements for the uniform rectangle subdivision (4.12): M ¼ fðxi ; tj Þ : i ¼ 0; 1; . . . ; M; j ¼ 0; 1; . . . ; N g. We begin by constructing the subspace W 1 ½M which consists of piecewise bilinear
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
1409
polynomials interpolating u 2 H 2 ½M . To do so, we first construct bilinear interpolation basis functions /mn(x, t) defined on M which satisfy /m;n ðxi ; tj Þ ¼ dmi dnj ;
i ¼ 0; 1; . . . ; M;
j ¼ 0; 1; . . . ; N :
ð5:5Þ
Note that /m,n(x, t) can be obtained from x-direction and t-direction linear interpolation basis functions. Therefore, we now construct x-direction and t-direction linear interpolation basis functions respectively. First consider x-direction linear interpolation basis functions /mð1Þ ðxÞ on [a, b] which satisfy /mð1Þ ðxi Þ ¼ dmi ;
i ¼ 0; 1; . . . ; M:
We list the basis functions as follows: 0 1 xx ; x0 6 x 6 x1 ; ð1Þ Dx /0 ðxÞ ¼ 0; otherwise; ð1Þ /M ðxÞ
¼
M 1 þ xx ; xM1 6 x 6 xM ; Dx 0; otherwise;
and for m = 1, 2, . . . , M 1 8 xxm > < 1 þ Dx ; xm1 6 x 6 xm ; m ; xm 6 x 6 xmþ1 ; /mð1Þ ðtÞ ¼ 1 xx Dx > : 0; otherwise:
ð5:6Þ
ð5:7Þ
ð5:8Þ
ð5:9Þ
Second consider t-direction linear interpolation basis functions /nð2Þ ðtÞ on [c, d] which satisfy /nð2Þ ðtj Þ ¼ dnj ;
j ¼ 0; 1; . . . ; N :
We list the basis functions as follows: 0 ; t0 6 t 6 t1 ; 1 tt ð2Þ Dt /0 ðtÞ ¼ 0; otherwise; ð2Þ
/N ðtÞ ¼
N 1 þ tt ; tN 1 6 t 6 tN ; Dt
0;
otherwise;
and for n = 1, 2, . . . , N 1 8 ttn > < 1 þ Dt ; tn1 6 t 6 tn ; ð2Þ n /n ðtÞ ¼ 1 tt ; tn 6 t 6 tnþ1 ; Dt > : 0; otherwise:
ð5:10Þ
ð5:11Þ
ð5:12Þ
ð5:13Þ
Using (5.7)–(5.13), we obtain ð2Þ /m;n ðx; tÞ ¼ /ð1Þ m ðxÞ/n ðtÞ;
m ¼ 0; 1; . . . ; M; n ¼ 0; 1; . . . ; N :
ð5:14Þ
1410
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
Using (5.14), we can obtain the expression of P ðx; tÞ; Qðx; tÞ 2 W 1 ½M : P ðx; tÞ ¼
M X N X
pm;n /m;n ðx; tÞ;
m¼0 n¼0
Qðx; tÞ ¼
M X N X
ð5:15Þ qm;n /m;n ðx; tÞ;
m¼0 n¼0
where pm,n p(xm,tn), qm,n q(xm,tn). In the subspace W 1 ½M , the action functional (4.3) becomes S ððP ðx; tÞ; Qðx; tÞÞÞ Z b Z d 1 2 P x þ Q2x þ PQt QP t V ðP 2 þ Q2 Þ dx ^ dt ¼ 2 a c M 1 X N 1 Z xiþ1 Z tjþ1 X 1 2 2 2 2 P þ Qx þ PQt QP t V ðP þ Q Þ dx ^ dt ¼ 2 x xi tj i¼0 j¼0 ¼
M 1 X N 1 X i¼0
Lðpi;j ; qi;j ; piþ1;j ; qiþ1;j ; piþ1;jþ1 ; qiþ1;jþ1 ; pi;jþ1 ; qi;jþ1 ÞDxDt;
j¼0
ð5:16Þ where Lðpi;j ; qi;j ; piþ1;j ; qiþ1;j ; piþ1;jþ1 ; qiþ1;jþ1 ; pi;jþ1 ; qi;jþ1 Þ Z xiþ1 Z tjþ1 1 1 2 ½P x þ Q2x þ PQt QP t V ðP 2 þ Q2 Þ dx ^ dt ¼ DxDt xi 2 tj 1 ¼ ðI 1 þ I 2 I 3 Þ; 2
ð5:17Þ
where I1 ¼
1 DxDt
Z
xiþ1
Z
xi
tjþ1
P 2x þ Q2x dx ^ dt
tj
2 t t t tj j 1 ðpiþ1;j pi;j Þ þ ðpiþ1;jþ1 pi;jþ1 Þ Dt Dt tj 2 t tj t tj þ 1 dt ðqiþ1;j qi;j Þ þ ðqiþ1;jþ1 qi;jþ1 Þ Dt Dt
1 ¼ 2 Dx Dt
Z
tjþ1
2
¼
2
ðpiþ1;j pi;j Þ þ ðpiþ1;j pi;j Þðpiþ1;jþ1 pi;jþ1 Þ þ ðpiþ1;jþ1 pi;jþ1 Þ 3Dx2 2 ðqiþ1;j qi;j Þ þ ðqiþ1;j qi;j Þðqiþ1;jþ1 qi;jþ1 Þ þ ðqiþ1;jþ1 qi;jþ1 Þ2 þ ; 3Dx2 ð5:18Þ
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
1 I2 ¼ DxDt
Z
xiþ1
Z
xi
1411
tjþ1
ðPQt QP t Þ dx ^ dt
tj
Z xiþ1 Z tjþ1 1 x xi P ½i;j 1 ðqijþ1 qij Þ DxDt xi Dx tj x x i þ ðqiþ1jþ1 qiþ1j Þ dx dt Dx Z xiþ1 Z tjþ1 1 x xi Q½i;j 1 ðpijþ1 pij Þ DxDt xi Dx tj x x i þ ðpiþ1jþ1 piþ1j Þ dx dt Dx 1 ðq qi;j Þð2pi;j þ 2pi;jþ1 þ piþ1j þ piþ1;jþ1 Þ ¼ 12Dt i;jþ1 1 ðq qiþ1;j Þðpi;j þ pi;jþ1 þ 2piþ1j þ 2piþ1;jþ1 Þ þ 12Dt iþ1;jþ1 1 ðp pi;j Þð2qi;j þ 2qi;jþ1 þ qiþ1j þ qiþ1;jþ1 Þ 12Dt i;jþ1 1 ðp piþ1;j Þðqi;j þ qi;jþ1 þ 2qiþ1j þ 2qiþ1;jþ1 Þ; 12Dt iþ1;jþ1
¼
ð5:19Þ
and I3 ¼
Z
xiþ1
Z
xiþ1
Z
xi
¼
Z
xi
tjþ1
V ðP 2 þ Q2 Þ dx ^ dt
tj tjþ1
V ððP ½i;j Þ2 þ ðQ½i;j Þ2 Þ dx ^ dt:
ð5:20Þ
tj
Here P[i,j] and Q[i,j] denote the restriction of P(x, t) and Q(x, t) on [xi, xi+1; t j, tj+1] respectively. They are given by x x x xi t tj t tj i P ½i;j ¼ 1 1 pij þ 1 piþ1j Dx Dt Dx Dt x x x xi t tj t tjþ1 i þ 1 pijþ1 þ 1þ piþ1jþ1 ; Dx Dt Dx Dt x x x xi t tj t tj i Q½i;j ¼ 1 1 qij þ 1 qiþ1j Dx Dt Dx Dt x xi t t j x xi t tj þ 1 qijþ1 þ qiþ1jþ1 : Dx Dt Dx Dt For the discrete Lagrangian (5.17), the discrete Euler–Lagrange equation (4.15) takes the form
1412
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
pðiÞ;jþ1 pðiÞ;j1 qiþ1;ðjÞ 2qi;ðjÞ þ qi1;ðjÞ 1 þ þ N 1 ¼ 0; 4 2Dt Dx2 qðiÞ;jþ1 qðiÞ;j1 piþ1;ðjÞ 2pi;ðjÞ þ pi1;ðjÞ 1 N 2 ¼ 0; 4 2Dt Dx2
ð5:21Þ
where pi1;l þ 4pi;l þ piþ1;l qi1;l þ 4qi;l þ qiþ1;l ; qðiÞ;l ¼ ; l ¼ j 1; j þ 1; 6 6 pr;j1 þ 4pr;j þ pr;jþ1 qr;j1 þ 4qr;j þ qr;jþ1 ; qr;ðjÞ ¼ ; r ¼ i 1; i; i þ 1; pr;ðjÞ ¼ 6 6 ! oI 3 ðwi;j Þ oI 3 ðwi1;j Þ oI 3 ðwi1;j1 Þ oI 3 ðwi;j1 Þ N1 ¼ þ þ þ ; opi;j opi;j opi;j opi;j ! oI 3 ðwi;j Þ oI 3 ðwi1;j Þ oI 3 ðwi1;j1 Þ oI 3 ðwi;j1 Þ N2 ¼ þ þ þ ; oqi;j oqi;j oqi;j oqi;j wi;j ¼ pi;j ; qi;j ; piþ1;j ; qiþ1;j ; piþ1;jþ1 ; qiþ1;jþ1 ; pi;jþ1 ; qi;jþ1 : pðiÞ;l ¼
Using w = p + iq, (5.21) becomes an integrator for w: i
wðmÞ;nþ1 wðmÞ;n1 wmþ1;ðnÞ 2wm;ðnÞ þ wm1;ðnÞ 1 þ þ N ¼ 0; 4 2Dt Dx2
ð5:22Þ
where wm1;r þ 4wm;r þ wmþ1;r ; r ¼ n 1; n þ 1; 6 wl;n1 þ 4wl;n þ wl;nþ1 ; l ¼ m 1; m; m þ 1; wl;ðnÞ ¼ 6 N ¼ N 2 þ iN 1 : wðmÞ;r ¼
Now we explore the relationship between the integrator (5.3) and the integrator (5.22). Using mid-point numerical integration formula in (5.18)– (5.20), we obtain piþ1;jþ12 pi;jþ12 2 qiþ1;jþ12 qi;jþ12 2 I1 þ ; Dx Dx qiþ1;jþ1 qiþ12;j piþ1;jþ1 piþ12;j qiþ12;jþ12 2 ; I 2 piþ12;jþ12 2 Dt Dt I 3 V 0 ðpiþ12;jþ12 Þ2 þ ðqiþ12;jþ12 Þ2 : Under these approximations, the discrete Lagrangian (5.17) becomes the discrete Lagrangian (5.1). Correspondingly, the variational integrator (5.22) becomes the variational integrator (5.3). This fact indicates that the multisymplectic integrator is closely related to the finite element method, which
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
1413
encourages to further explore the properties of the multisymplectic integrator by using the results from the finite element methods. In fact, we can prove the convergence of the integrator (2.8) by the convergence of the integrator (5.22) which is guaranteed by the error estimation from finite element analysis. 6. Numerical experiments For the numerical experiments, we consider the cubic NLS equation 2
iwt þ wxx þ 2jwj w ¼ 0;
ð6:1Þ
with periodic boundary condition w(0,t) = w(L,t) and the initial condition wðx; 0Þ ¼ Að1 þ 0:05 cosðlxÞÞ: ð6:2Þ pffiffiffi Here L ¼ 4 2p, l = 2p/L and A is a parameter. The initial condition (6.2) is in the vicinity of the homoclinic orbit associated with the NLS equation. The homoclinic structure causes severe sensitivities to numerical discretizations. The parameter A serves to change the complexity of the homoclinic structure. Larger A will result in higher complexity of the homoclinic structure. For details, see [15] and the references therein. Now we perform numerical experiments using the multisymplectic integrator (2.8). In our first experiment, we take A = 0.75. The space step size L Dx ¼ 127 and the time step size D t = 0.01 are used. We employ the iteration method introduced in [24] to solve (2.8). The iteration tolerance 1010 is used. We first check the discrete conservation law (3.4) numerically. Let 2 PM1 En ¼ m¼0 wmþ12;n . Then from (3.4), we have Enþ1 En ¼ 0:
ð6:3Þ
From our experiment, we see that the equation (6.3) holds exactly (in the sense of machine accuracy) in spite of the iteration error 1010 (see Fig. 2). The temporal development of waveforms is shown in Fig. 3. From the plots, we see that the spatial symmetry is well preserved. The waveforms are smooth and no instability occurs. Now we increase A to A = 1, which increases the
1
x 10
En+1 En
0.5
0
0
1000
2000
3000
4000
n
Fig. 2. The discrete conservation law (6.3).
5000
1414
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417 t=5 3
2
2
|ψ|
|ψ|
t=0 3
1
0
1
0
0
8.9
17.8
0
3
3
2
2
1
0
1
8.9
17.8
0
x 3
2
2
|ψ|
|ψ|
17.8
t=25
3
1
0
1
0
0
8.9
17.8
0
8.9
x
x
t=30
t=35
3
3
2
2
|ψ|
|ψ|
8.9
x
t=20
1
0
17.8
1
0
0
8.9
17.8
0
x
8.9
17.8
x
t=40
t=45
3
3
2
2
|ψ|
|ψ|
17.8
0
0
1
0
1
0
0
8.9
17.8
0
8.9
x
x
t=47
t=50
3
3
2
2
|ψ|
|ψ|
8.9
x t=15
|ψ|
|ψ|
x t=10
1
0
17.8
1
0
0
8.9
x
17.8
0
8.9
x
Fig. 3. The temporal development of waveforms. A = 0.75.
17.8
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417 t=5
3
3
2
2
|ψ|
|ψ|
t=0
1
0
8.9
17.8
0
x t=15
3
3
2
2
1
0
1
8.9
17.8
0
x
17.8
t=25
3
3
2
2
|ψ|
|ψ|
8.9
x
t=20
1
0
1
0
0
8.9
17.8
0
x
8.9
17.8
x
t=30
t=35
3
3
2
2
|ψ|
|ψ|
17.8
0
0
1
0
1
0
0
8.9
17.8
0
8.9
x
x
t=40
t=45
3
3
2
2
|ψ|
|ψ|
8.9
x t=10
|ψ|
|ψ|
1
0
0
1
0
17.8
1
0
0
8.9
17.8
0
x
8.9
17.8
x
t=47
t=50
3
3
2
2
|ψ|
|ψ|
1415
1
0
1
0
0
8.9
x
17.8
0
8.9
x
Fig. 4. The temporal development of waveforms. A = 1.0.
17.8
1416
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
complexity and sensitivity of the homoclinic structure. In this case, we use the L space step size Dx ¼ 256 and the time step size Dt = 0.01. Waveforms at different times are shown in Fig. 4. Again, the spatial symmetry is well preserved and no instability occurs. Notice that, in the first experiment, we use an odd number of spatial points while in the second experiment, an even number of spatial points are used. The numerical experiments indicate that the stability of the multisymplectic integrator (2.8) does not depend on the oddness of the spatial points.
7. Conclusions In this paper, we present a multisymplectic integrator which is derived from the multisymplectic formulation of the periodic nonlinear Schro¨dinger equation. The integrator satisfies a discrete multisymplectic conservation law. Its accuracy is proved. By introducing a norm, we investigate its nonlinear stability in some sense. We prove that this multisymplectic integrator is a reduction of a variational integrator which is derived from the MPS framework of discrete multisymplectic field theory. Furthermore, we prove that this integrator is also an approximation of another variational integrator which is derived from the finite element method. Numerical experiments are reported. Our results in this paper can be naturally generalized to those for higher dimensional NLS equation: 2 iwt þ wxx þ wyy þ V 0 jwj w ¼ 0; and 2 iwt þ wxx þ wyy þ wzz þ V 0 jwj w ¼ 0; with corresponding periodic boundary conditions.
Acknowledgements This work is supported by National Natural Science Foundation of China under Grant No. 40474047 and also partially by Chinese Academy of Sciences with Key Project of Knowledge Innovation KZCX1-SW-18.
References [1] M.J. Ablowitz, J.F. Ladik, A nonlinear difference scheme and inverse scattering, Stud. Appl. Math 55 (1976) 213–229.
J.-B. Chen / Appl. Math. Comput. 170 (2005) 1394–1417
1417
[2] S.C. Brenner, L.R. Scott, The Mathematical Theory of Finite Element Methods, SpringerVerlag, New York, 1994. [3] T.J. Bridges, Multi-symplectic structures and wave propagation, Math. Proc. Cam. Phil. Soc. 121 (1997) 147–190. [4] T.J. Bridges, S. Reich, Multi-symplectic integrators: numerical schemes for Hamiltonian PDEs that conserve symplecticity, Phys. Lett. A 284 (2001) 184–193. [5] J.B. Chen, Total variation in discrete multisymplectic field theory and multisymplectic energy momentum integrators, Lett. Math. Phys. 61 (2002) 63–73. [6] J.B. Chen, H.Y. Guo, K. Wu, Total variation in Hamiltonian formalism and symplecticenergy integrators, J. Math. Phys. 44 (2003) 1688–1702. [7] J.B. Chen, M.Z. Qin, Y.F. Tang, Symplectic and multisymplectic methods for the nonlinear Schro¨dinger equation, Comput. Math. Appl. 43 (2002) 1095–1106. [8] J.B. Chen, A multisymplectic variational integrator for the nonlinear Schro¨dinger equation, Numer. Meth. Part. Diff. Eq. 18 (2002) 523–536. [9] J.B. Chen, Variational integrators for discrete multisymplectic field theory based on the finite element method, preprint. [10] P.G. Ciarlet, Numerical Analysis of the Finite Element Method, Les Presses, de LUniversite´ de Montreal, Quebec, 1976. [11] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, New York, Oxford, 1978. [12] M. Delfour, M. Fortin, G. Payre, Finite-difference solutions of a nonlinear Schro¨dinger equation, J. Comput. Phys. 44 (1981) 277–288. [13] P.L. Garcia, The Poincare´–Cartan invariant in the calculus of variations, Symp. Math. 14 (1974) 219–246. [14] A. Hasegawa, Optical Solitons in Fibers, Springer-Verlag, Berlin, 1989. [15] B.M. Herbst, F. Varadi, M.J. Ablowitz, Symplectic methods for the nonlinear Schro¨dinger equation, Math. Comput. Simul. 37 (1994) 353–369. [16] J.L. Hong, Y. Liu, H. Munthe-Kass, A. Zanna, On a multisymplectic scheme for Schro¨dinger equations with variable coefficients, preprint. [17] G.L. Lamb, Elements of Soliton Theory, John Wiley & Sons, New York, 1980. [18] J.E. Marsden, G.W. Patrick, S. Shkoller, Multisymplectic geometry, variational integrators, and nonlinear PDEs, Comm. Math. Phys. 199 (1998) 351–395. [19] R. McLachlan, Symplectic integration of Hamiltonian wave equations, Numer. Math. 66 (1994) 465–492. [20] P.J. Olver, Applications of Lie Groups to Differential Equations, second ed., Springer-Verlag, New York, 1993. [21] S. Reich, Multi-symplectic Runge–Kutta collocation methods for Hamiltonian wave equations, J. Comput. Phys. 157 (2000) 473–499. [22] C.M. Schober, Symplectic integrators for Ablowitz–Ladik discrete nonlinear Schro¨dinger equation, Phys. Lett. A 259 (1999) 140–151. [23] Y.F. Tang, L. Va´zquez, F. Zhang, V.M. Pe´rez-Garcı´a, Symplectic methods for the nonlinear Schro¨dinger equation, Comput. Math. Appl. 32 (1996) 73–83. [24] S.B. Wineberg, J.F. McGrath, E.F. Gabl, L.R. Scott, C.E. Southwell, Implicit spectral methods for wave propagation problems, J. Comput. Phys. 97 (1991) 311–336.