Explicit Multi-Symplectic Extended Leap-frog Methods ... - Purdue Math
Explicit Multi-Symplectic Extended Leap-frog Methods for Hamiltonian Wave Equations Wei Shia , Xinyuan Wua , Jianlin Xiab,∗ a Department
of Mathematics, Nanjing University; State Key Laboratory for Novel Software Technology at Nanjing University, Nanjing 210093, P.R.China b Department of Mathematics, Purdue University, West Lafayette, IN 47907, USA
Abstract In this paper, we study the integration of Hamiltonian wave equations whose solutions have oscillatory behaviors in time and/or space. We are mainly concerned with the research for multi-symplectic extended Runge-Kutta-Nystr¨om (ERKN) discretizations and the corresponding discrete conservation laws. We first show that the discretizations to the Hamiltonian wave equations using two symplectic ERKN methods in space and time respectively lead to an explicit multi-symplectic integrator (Eleap-frogI). Then we derive another multi-symplectic discretization using a symplectic ERKN method in time and a symplectic Partitioned Runge-Kutta method, which is equivalent to the well-known St¨ormer-Verlet method in space (Eleap-frogII). These two new multi-symplectic schemes are extensions of the leap-frog method. The numerical stability and dispersive properties of the new schemes are analyzed. Numerical experiments with comparisons are presented, where the two new explicit multi-symplectic methods and the leap-frog method are applied to the linear wave equation and the sine-Gordon equation. The numerical results confirm the superior performance and some significant advantages of our new integrators in the sense of structure preservation. Keywords: Multi-symplectic integrators, Hamiltonian wave equations, extended leap-frog methods, conservation laws, dispersive properties
Mathematics Subject Classification (2000): 35L05, 37M15, 65L06, 65M12, 65P10, 70S10
1. Introduction It is well-known that symplectic integrators are robust, efficient, and accurate in preserving the longtime behavior of the solutions of Hamiltonian ordinary differential equations (ODEs) [1]. The basic idea of a symplectic integrator is that the numerical scheme is designed to preserve the symplectic form at each time step. Recently, it is shown that many conservative partial differential equations (PDEs), such as wave equations, generalized KdV and nonlinear Schr¨odinger equations, etc., allow for a structure similar to the symplectic structure of Hamiltonian ODEs, called the multi-symplectic formulation (see, e.g., [2, 3, 4]). For example, in [4], Marsden and Shkoller develop the multi-symplectic structure of Hamiltonian PDEs from a Lagrangian formulation using a variational principle. It is natural to design integrators for Hamiltonian PDEs which preserve a semi-discretization of the symplectic form associated with the infinite-dimensional evolution equation. In the past decade, rapid progress I The research of Wei Shi and Xinyuan Wu was supported in part by the Natural Science Foundation of China under Grant 10771099 and by the Specialized Research Foundation for the Doctoral Program of Higher Education under Grant 20100091110033, by the 985 Project at Nanjing University under Grant 9112020301, by A Project Founded by the Priority Academic Program Development of Jiangsu Higher Education Institutions, and by the University Postgraduate Research and Innovation Project of Jiangsu Province 2012 under Grant CXLX12− 0033. II The research of Jianlin Xia was supported in part by NSF grants DMS-1115572 and CHE-0957024. ∗ Corresponding author. Email address: [email protected], [email protected], [email protected] (Jianlin Xia)
1
W. Shi, X. Wu, and J. Xia
2
has been made on the development of multi-symplectic integrators of Runge-Kutta type for various Hamiltonian PDEs (see, e.g., [5, 6, 7, 8, 9, 10, 11, 12]). Unfortunately, most of the existing multi-symplectic integrators fail to take into account the oscillatory behavior of the problems of interest, which is important for the solution that has the intrinsic oscillatory behavior. It is known that most of integrators adapted to the oscillatory problems are frequency-dependent, such as ARKN methods [13, 14, 15, 16, 17], ERKN methods [18, 19, 20, 21], and trigonometrically/exponentiallyfitted methods [22, 23, 24, 25]. For nonlinear Schr¨odinger equations, taking the frequency of the problem into account, several useful trigonometrically/exponentially-fitted methods are proposed (see, e.g., [26, 27, 28, 29, 30, 31, 32]). Some authors, such as Gonz´alez, et al. [33] and Franco [13], are interested in numerical integrators of Runge-Kutta-Nystr¨om type adapted to initial value problems of the form ¡ ¢ ½ 00 y (t) + ω2 y(t) = f y(t) , t ∈ [t0 , T ], (1) y(t0 ) = y0 , y0 (t0 ) = y00 , where ω is the main frequency and ω > 0. More recently, a new family of extended Runge-Kutta-Nystr¨om (ERKN) methods are proposed by Yang, et al. [18] for one-dimensional perturbed oscillators and systems of these oscillators with ω2 being a diagonal matrix. Furthermore, with the methodology in [18], a standard form of multidimensional ERKN integrators for the general system (1) is formulated by Wu, et al. [19], where ω2 is a symmetric positive semi-definite matrix and implicitly contains the frequencies of the problem. The corresponding order conditions are also derived by the authors based on the B-series theory associated with the extended Nystr¨om trees (EN-trees) [19]. These methods preserve the oscillatory feature of the unperturbed oscillators, not only for the updates but also for the internal stages. Numerical results have shown that ERKN methods are superior to other methods for oscillatory systems. In this paper, we investigate the multi-symplectic ERKN (MSERKN) methods for wave equations as Hamiltonian PDEs. In Section 2, we discuss the conservation laws for Hamiltonian wave equations. In Section 3, we show that the discretization in time and space for Hamiltonian wave equations with two symplectic ERKN (SERKN) methods leads to a multi-symplectic integrator. In Section 4, explicit multi-symplectic schemes are constructed based on the SERKN methods and a symplectic Partitioned Runge-Kutta (SPRK) method. Section 5 is devoted to the numerical experiments, including the dispersive properties of the schemes proposed in this paper. Some conclusions are given in the last section. 2. Conservation laws and multi-symplectic structures of wave equations 2.1. Multi-symplectic conservation laws and multi-symplectic integrators A multi-Hamiltonian PDE in one time dimension and one space dimension is a PDE that can be formulated in the following way (see [34]). Let M and K be any skew-symmetric real n × n matrices (n ≥ 3) and let S : Rn → R be any smooth function of the state variable z. Then a system of the following form is called a Hamiltonian system on a multi-symplectic structure: Mzt + Kz x = ∇z S (z), z ∈ Rn , (x, t) ∈ R2 ,
(2)
where ∇z is the gradient operator corresponding to the standard inner product on Rn denoted by h·, ·i. The Hamiltonian system (2) is multi-symplectic in the sense that it is associated with two-forms: ζ=
1 1 dz ∧ Mdz, κ = dz ∧ Kdz, 2 2
which define a space-time symplectic structure governed by the multi-symplectic conservation law ∂t ζ + ∂ x κ = 0. Let L be the length of the spatial domain. Then integrating (3) over the spatial domain leads to d dt
Z
L/2
ζdx + κ| x=L/2 − κ| x=−L/2 = 0. −L/2
(3)
W. Shi, X. Wu, and J. Xia
3
This equation, under appropriate boundary conditions such as z(L/2, t) = z(−L/2, t), results in the conservad tion of global symplecticity: ζ(z) = 0, where dt Z L/2 ζ(z) = ζ(z)dx. (4) −L/2
ζ(z) is a reminder indicating that the two-form ζ is integrated over all space. A properly chosen semidiscretization of this system may lead to a system of Hamiltonian ODEs with a symplectic form given by a discrete version of (4). An overview of this viewpoint in terms of symplectic time-integration methods is given by McLachlan [35]. Two other conservation laws for the energy and the momentum respectively follow from (3). The local energy conservation law is ∂t E(z) + ∂ x F(z) = 0, (5) 1 1 where E(z) = S (z)− zT Kz x is the energy density and F(z) = zT Kzt is the energy flux. The local momentum 2 2 conservation law is ∂t I(z) + ∂ xG(z) = 0, (6) 1 1 where I(z) = zT Mz x and G(z) = S (z) − zT Mzt are the momentum density and the momentum flux, 2 2 respectively. Similar to the conservation of multi-symplecticity, integrating the local conservation laws (5) and (6) over the spatial domain in Rn with periodic boundary conditions leads to global conservation laws d d ε(z) = 0 and I(z) = 0, (7) dt dt R L/2 R L/2 where ε(z) = −L/2 E(z)dx and I(z) = −L/2 I(z)dx. We note that the global conservation of energy or momentum (7) is necessary but not sufficient for the local conservation of energy or momentum. When it comes to the numerical integration, the symplecticity of numerical integrators should be determined. Definition 2.1. [34] A numerical scheme is called multi-symplectic if it preserves a discrete conservation law of (3). An alternative definition of multi-symplectic integrators based on a discrete form of the symplectic conservation law is suggested by Bridges and Reich [3]. It has been shown that popular methods such as the centered Preissman scheme and the leap-frog method are multi-symplectic and that such schemes have remarkable properties of conserving local energy and momentum. 2.2. Conservation laws for wave equations Consider the following nonlinear wave equation (see [36]): utt − u xx + V 0 (u) = 0,
(8)
where V(u) is certain smooth function and u is a scalar function. Introducing two new variables v = ut and w = u x , we can write (8) as vt − w x + V 0 (u) = 0, v = ut , w = ux . By taking z = (u, v, w)T ∈ R3 , we can obtain a first order PDE system of the abstract form (2) corresponding to (8) (see, e.g., [11]), with 0 0 1 0 −1 0 0 0 , K = 0 0 0 , M= 1 −1 0 0 0 0 0
4
W. Shi, X. Wu, and J. Xia
1 and the Hamiltonian is S (z) = (v2 − w2 ) + V(u). As mentioned in Subsection 2.1, this system has the 2 multi-symplectic conservation law (3), where ζ and κ are the pre-symplectic forms 1 dz ∧ Mdz = dv ∧ du, 2 1 κ = dz ∧ Kdz = du ∧ dw. 2
ζ=
(9) (10)
By applying the general approach of deriving the local energy and momentum conservation law to the nonlinear wave equation (8), as discussed in the previous subsection, it follows that the local energy conserva1 tion law is (5) with E(z) = (v2 + w2 ) + V(u) and F(z) = −vw. Meanwhile, the local momentum conservation 2 1 law is (6) with I(z) = vw and G(z) = − (v2 + w2 ) + V(u) (see [11]). 2 3. Extended RKN discretization of wave equations 3.1. Extended RKN methods for ODEs In [18], Yang, et al. investigate ERKN methods for the IVP (1) of the following form : s P Yi = φ0 (ci ν)yn + ci φ1 (ci ν)hy0n + h2 a¯ i j (ν) f (tn + c j h, Y j ), i = 1, . . . , s, j=1 s P yn+1 = φ0 (ν)yn + φ1 (ν)hy0n + h2 b¯ i (ν) f (tn + ci h, Yi ), i=1 s 0 0 y = −νwφ1 (ν)yn + φ0 (ν)y + h P bi (ν) f (tn + ci h, Yi ), n+1
n
(11)
i=1
where a¯ i j (ν), bi (ν) and b¯ i (ν), i, j = 1, . . . , s, are assumed to be even functions of ν = hω, φ0 (ν) = cos(ν), and sin(ν) φ1 (ν) = . We denote (11) by {ci , bi , b¯ i , a¯ i j , ω} in this paper. ν Very recently, ERKN methods are considered for Hamiltonian ODEs, whose solutions exhibit a periodic or oscillatory character. The coefficients of these methods depend on the product of the step size h and the fitted frequency ω, and the methods can be expressed in Butcher tableau as :
c
¯ A(ν) b¯ T (ν) bT (ν)
=
c1 .. .
a¯ 11 (ν) .. .
cs
a¯ 1s (ν) .. .
a¯ s1 (ν)
... .. . ···
b¯ 1 (ν)
···
b¯ s (ν)
b1 (ν)
···
b s (ν)
a¯ ss (ν)
The symplecticity conditions for ERKN methods are given by the next theorem. Theorem 3.1. (Wu, et al. [37, 20]) If the coefficients of an s-stage ERKN method satisfy the following conditions ¶ µ ¶ µ 2 b¯ i (ν) φ0 (ν) + ci ν φ1 (ν)φ1 (ci ν) = bi (ν) φ1 (ν) − ci φ0 (ν)φ1 (ci ν) , φ0 (ci ν) φ0 (ci ν) iµ = 1, . . . , s, ¶ a¯ i j (ν)φ0 (ν) ν2 φ1 (ν) (12) ¯ − b¯ i (ν)¯ai j (ν) b (ν) b (ν) − i j φ (c ν) φ (c ν) 0 i 0 i ¶ µ a¯ ji (ν)φ0 (ν) ν2 φ1 (ν) − b¯ j (ν)¯a ji (ν) , i, j = 1, . . . , s, = b j (ν) b¯ i (ν) − φ0 (c j ν) φ0 (c j ν) then the method is symplectic.
5
W. Shi, X. Wu, and J. Xia
R 3.1. When ω → 0, the scheme (11) reduces to a classical RKN method in Hairer, et al. [38]. In this case, the symplecticity conditions of (12) become ½ b¯ i = (1 − ci )bi , i = 1, . . . , s, bi (b¯ j − a¯ i j ) = b j (b¯ i − a¯ ji ), i ≤ j = 1, . . . , s, which are obtained by Suris [39] for the symplecticity of an s-stage RKN method. 3.2. Multi-symplectic extended RKN discretization In order to derive the new MSERKN methods, we first consider the discretization of the equation (8) by employing the following two different SERKN methods in temporal and spatial directions, respectively: c1 .. .
a¯ 11 (νt ) .. .
a¯ 1r (νt ) .. .
c˜ 1 .. .
aˆ 11 (ν x ) .. .
a¯ r1 (νt )
... .. . ···
cr
a¯ rr (νt )
c˜ s
aˆ s1 (ν x )
b¯ 1 (νt )
···
b¯ r (νt )
b1 (νt )
···
br (νt )
bˆ 1 (ν x ) b˜ 1 (ν x )
... .. . ···
aˆ 1s (ν x ) .. .
···
bˆ s (ν x ) b˜ s (ν x )
···
aˆ ss (ν x )
Theorem 3.2. Assume that a method results from the temporal discretization of the wave equation (8) with an r-stage SERKN scheme {ci , bi , b¯ i , a¯ i j , ωt } and spatial discretization with an s-stage SERKN scheme {˜ci , b˜ i , bˆ i , aˆ i j , ω x }, where the coefficients satisfy the symplectic conditions ¶ µ ¶ µ ci φ0 (νt )φ1 (ci νt ) ci νt 2 φ1 (νt )φ1 (ci νt ) ¯ = bi (νt ) φ1 (νt ) − , bi (νt ) φ0 (νt ) + φ0 (ci νt ) φ0 (ci νt ) iµ= 1, . . . , r, ¶ a¯ i j (νt )φ0 (νt ) νt 2 φ1 (νt ) ¯ (13) bi (νt ) b j (νt ) − − b¯ i (νt )¯ai j (νt ) φ0 (ci νt ) µ φ0 (ci νt ) ¶ a¯ ji (νt )φ0 (νt ) νt 2 φ1 (νt ) = b j (νt ) b¯ i (νt ) − − b¯ j (νt )¯a ji (νt ) , φ0 (c j νt ) φ0 (c j νt ) i, j = 1, . . . , r, and
¶ µ ¶ µ 2 ˆ i (ν x ) φ0 (ν x ) + c˜ i ν x φ1 (ν x )φ1 (˜ci ν x ) = b˜ i (ν x ) φ1 (ν x ) − c˜ i φ0 (ν x )φ1 (˜ci ν x ) , b φ0 (˜ci ν x ) φ0 (˜ci ν x ) iµ = 1, . . . , s, ¶ aˆ i j (ν x )φ0 (ν x ) ν x 2 φ1 (ν x ) b˜ i (ν x ) bˆ j (ν x ) − − bˆ i (ν x )ˆai j (ν x ) φ0 (˜ci ν x ) µ φ0 (˜ci ν x ) ¶ a ˆ (ν )φ (ν ) ν x 2 φ1 (ν x ) ji x 0 x = b˜ j (ν x ) bˆ i (ν x ) − − bˆ j (ν x )ˆa ji (ν x ) , φ0 (˜c j ν x ) φ0 (˜c j ν x ) i, j = 1, . . . , s,
(14)
respectively. The method is said to be multi-symplectic in the sense that it satisfies the discrete multisymplectic conservation law µ ¶ s ¢ P 1 ¡˜ 2 n+1 n n ˆ bi φ0 (ν x ) + bi ν x φ1 (ν x ) (dvn+1 m,i ∧ dum,i − dvm,i ∧ dum,i )h ci ν x ) i=1 µ φ0 (˜ ¶ (15) r ¢ P 1 ¡ n,k n,k n,k + bk φ0 (νt ) + b¯ k νt2 φ1 (νt ) (dun,k ∧ dw − du ∧ dw )τ = 0, m m m+1 m+1 k=1 φ0 (ck νt ) where τ and h are the step-size of time and space, respectively.
6
W. Shi, X. Wu, and J. Xia
Proof. (8) as
We first apply an s-stage SERKN method with the spatial discretization and rewrite equation s P Ui = φ0 (˜ci ν x )um + c˜ i φ1 (˜ci ν x )hwm + h2 aˆ i j (ν x ) f j , j=1 s P um+1 = φ0 (ν x )um + φ1 (ν x )hwm + h2 bˆ i (ν x ) fi , i=1 s P wm+1 = −ν x ω x φ1 (ν x )um + φ0 (ν x )wm + h b˜ i (ν x ) fi ,
i = 1, . . . , s,
i=1
where fi = f (Ui ) = ω2x Ui + ∂ xx Ui . Having φ0 (ν) = cos(ν) and φ1 (ν) = dum+1 ∧ dwm+1 = dum ∧ dwm + h s ¡ P
sin(ν) in mind, we proceed with ν
s ¡ ¢ P b˜ i (ν x )φ0 (ν x ) + bˆ i (ν x )ν2x φ1 (ν x ) dum ∧ d fi i=1
s ¢ P b˜ i (ν x )φ1 (ν x ) − bˆ i (ν x )φ0 (ν x ) dwm ∧ d fi + h3 bˆ i (ν x )b˜ j (ν x )d fi ∧ d f j i=1 i, j=1 ¶ µ s ¢ P 1 ¡˜ 2 ˆ bi (ν x )φ0 (ν x ) + bi (ν x )ν x φ1 (ν x ) dUi ∧ d fi = dum ∧ dwm + h φ0 (˜ci ν x ) i=1 µ ¶ s ¢ c˜ i φ1 (˜ci ν x ) P ¡ 2 2 ˜ ˆ ˜ ˆ −h − bi (ν x )φ1 (ν x )+ bi (ν x )φ0 (ν x ) dwm ∧d fi bi (ν x )φ0 (ν x )+ bi (ν x )ν x φ1 (ν x ) φ0 (˜ci ν x ) i=1 µ ¶ 2 s P b˜ i (ν x )φ0 (ν x ) + bˆ i (ν x )ν x φ1 (ν x ) 3 ˆ ˜ +h aˆ i j + bi (ν x )b j (ν x ) d fi ∧ d f j . φ0 (˜ci ν x ) i, j=1 2
+h
Together with (14), this is further transformed to ¶ µ s ¢ P 1 ¡˜ dum+1 ∧ dwm+1 = dum ∧ dwm + h bi (ν x )φ0 (ν x ) + bˆ i (ν x )ν2x φ1 (ν x ) dUi ∧ d fi ci ν x ) i=1 µ φ0 (˜ ¶ s ¢ P 1 ¡˜ = dum ∧ dwm + h bi (ν x )φ0 (ν x ) + bˆ i (ν x )ν2x φ1 (ν x ) dUi ∧ (ω2x dUi + ∂ xx dUi ). φ0 (˜ci ν x ) i=1 Then, we obtain the following semi-symplectic conservation law in space: ¶ µ s ¢ P 1 ¡˜ 2 ˆ bi (ν x )φ0 (ν x ) + bi (ν x )ν x φ1 (ν x ) dUi ∧ ∂ x dWi . dum+1 ∧ dwm+1 = dum ∧ dwm + h φ0 (˜ci ν x ) i=1
(16)
The next step is the discretization in time over a time interval [0, τ]. For simplicity, we set m = 0 and assume that xm = 0. We use the r-stage SERKN method and obtain r k 0 0 2P U = φ (c ν )u + c φ (c ν )τv + τ a¯ kl (νt )gli , k = 1, . . . , r, 0 k t k 1 k t i i i l=1 r P u1i = φ0 (νt )u0i + φ1 (vt )τv0i + τ2 b¯ k (νt )gki , k=1 r P v1i = −νt ωt φ1 (νt )u0i + φ0 (νt )v0i + τ bk (νt )gki . k=1
Here, we introduce the notation Uik ≈ u(˜ci h, ck τ), Vik ≈ v(˜ci h, ck τ), u1i ≈ u(˜ci h, τ), v1i ≈ v(˜ci h, τ), gki = g(Uik ) = ω2t Uik + ∂tt Uik . Likewise, we can obtain the identity du1i ∧ dv1i = du0i ∧ v0i + τ
r P k=1
µ
¶ ¢ 1 ¡ bk (νt )φ0 (νt ) + b¯ k (νt )νt2 φ1 (νt ) dUik ∧ ∂t dVik . φ0 (ck νt )
7
W. Shi, X. Wu, and J. Xia
This is the symplectic conservation law in time direction. For any x = xm,i = xm + c˜ i h, t = tn,k = tn + ck τ, the above identity can be transformed to ¶ r µ X ¢ 1 ¡ n,k n,k n+1 n+1 n n 2 ¯ dum,i ∧ dvm,i − dum,i ∧ dvm,i = τ bk (νt )φ0 (νt ) + bk (νt )νt φ1 (νt ) (dUm,i ∧ ∂t dVm,i ). (17) φ (c ν ) 0 k t k=1 Next, we rewrite (16) for t = tn,k as n,k n,k n,k dun,k m+1 ∧ dwm+1 − dum ∧ wm = h
s µ X i=1
¶ ¢ 1 ¡˜ n,k n,k bi (ν x )φ0 (ν x ) + bˆ i (ν x )ν2x φ1 (ν x ) (dUm,i ∧ ∂ x dWm,i ). (18) φ0 (˜ci ν x )
It finally follows from (17) and (18) that the discrete multi-symplectic conservation law (15) holds. Meanwhile, the fully-discrete scheme for (8) is given by s n,k n,k n,k n,k 2P U = φ (˜ c ν )u + c ˜ φ (˜ c ν )hw + h aˆ i j (ν x ) fm, 0 i x i 1 i x m m m,i j, j=1 s P n,k n,k um+1 = φ0 (ν x )un,k + φ1 (ν x )hwn,k + h2 bˆ i (ν x ) fm,i , m m i=1 s P n,k n,k n,k = −ν x ω x φ1 (ν x )un,k b˜ i (ν x ) fm,i , w m + φ0 (ν x )wm + h m+1
i = 1, . . . , s,
i=1
r P n,k Um,i = φ0 (ck νt )unm,i + ck φ1 (ck νt )τvnm,i + τ2 a¯ kl (νt )gn,l m,i , l=1 r n n 2P ¯ un+1 bk (νt )gn,k m,i = φ0 (νt )um,i + φ1 (νt )τvm,i + τ m,i , k=1 r P n n bk (νt )gn,k vn+1 m,i = −νt ωt φ1 (νt )um,i + φ0 (νt )vm,i + τ m,i ,
¤
(19) k = 1, . . . , r,
k=1
where we have used the following notation: n,k Um,i ≈ u(xm + c˜ i h, tn + ck τ),
un+1 ˜ i h, tn + τ), m,i ≈ u(xm + c
vn+1 ˜ i h, tn + τ), m,i ≈ v(xm + c
un,k m+1 ≈ u(xm + h, tn + ck τ),
wn,k m+1 ≈ w(xm + h, tn + ck τ),
n,k n,k n,k n,k 2 n,k 2 n,k 0 gn,k m,i = g(U m,i ) = ωt U m,i + ∂tt U m,i = ωt U m,i + ∂ xx U m,i − V (U m,i ), n,k n,k n,k n,k fm,i = f (Um,i ) = ω2x Um,i + ∂ xx Um,i .
The formula (15) can be understood as the approximation of the integral of (3) over the domain [xm , xm+1 ]× [tn , tn+1 ]. That is, (15) approximates Z xm+1 ¡ ¢ dv(x, tn+1 ) ∧ du(x, tn+1 ) − dv(x, tn ) ∧ du(x, tn ) dx Z
xm tn+1
+
¡
¢ du(xm+1 , t) ∧ dw(xm+1 , t) − du(xm , t) ∧ dw(xm , t) dt = 0,
(20)
tn
with SERKN methods for the evaluation of the two integrals. When the problems with periodic boundary conditions are considered, (15) has another interesting consequence. Let us take the sum of the identity (15) over all spatial grid points m = 1, . . . , M as: M nX s µ X
¶ ¢ 1 ¡˜ 2 n+1 n n ˆ bi (ν x )φ0 (ν x ) + bi (ν x )ν x φ1 (ν x ) (dvn+1 m,i ∧ dum,i − dvm,i ∧ dum,i )h φ (˜ c ν ) 0 i x m=1 i=1 ¶ r µ o X ¢ 1 ¡ n,k 2 n,k n,k ¯ + bk (νt )φ0 (νt ) + bk (νt )νt φ1 (νt ) (dun,k ∧ dw − du ∧ dw )τ = 0. m m m+1 m+1 φ0 (ck νt ) k=1
(21)
8
W. Shi, X. Wu, and J. Xia
Periodic boundary conditions in space imply M X r µ X
¶ ¢ 1 ¡ n,k n,k n,k bk (νt )φ0 (νt ) + b¯ k (νt )νt2 φ1 (νt ) (dun,k m+1 ∧ dwm+1 − dum ∧ dwm )τ φ (c ν ) 0 k t m=1 k=1 ¶ r µ X ¢ 1 ¡ n,k n,k n,k 2 ¯ = bk (νt )φ0 (νt ) + bk (νt )νt φ1 (νt ) (dun,k M+1 ∧ dw M+1 − du1 ∧ dw1 )τ = 0, φ (c ν ) 0 k t k=1 which in turn transforms (21) to M X s µ X
¶ ¢ 1 ¡˜ n+1 bi (ν x )φ0 (ν x ) + bˆ i (ν x )ν2x φ1 (ν x ) dvn+1 m,i ∧ dum,i φ (˜ c ν ) 0 i x m=1 i=1 ¶ M X s µ X ¢ 1 ¡˜ 2 ˆ = bi (ν x )φ0 (ν x ) + bi (ν x )ν x φ1 (ν x ) dvnm,i ∧ dunm,i . φ (˜ c ν ) 0 i x m=1 i=1 This is precisely the conservation of symplecticity in time. It is also the discretization of the total symplectic conservation (4) with (9). R 3.2. By abuse of notation, we denote the fitted temporal frequency and spatial frequency in the two SERKN methods {ci , bi , b¯ i , a¯ i j , ωt } and {˜ci , b˜ i , bˆ i , aˆ i j , ω x } by ωt and ω x , respectively. Obviously, when ωt → 0 and ω x → 0, the MSERKN methods become multi-symplectic RKN (MSRKN) methods. R 3.3. Theorem 3.2 shows that concatenation of SERKN methods can produce a multi-symplectic integrator for Hamiltonian wave equations. It is remarked that the above idea may apply equally to some other Hamiltonian PDEs. For Hamiltonian PDEs with oscillations in both space and time directions, two SERKN methods are used in space and time respectively. For the case that oscillation exists only in time direction, we propose the following multi-symplectic discretizations. Theorem 3.3. The wave equation (8) discretized by an SPRK method in space and an SERKN method in time has a discrete multi-symplectic conservation law. The proof is similar to that of Theorem 3.2. As an example, we use SERKN discretization in time and the following SPRK discretization in space (which can be regarded as the St¨ormer-Verlet formula):
|
0
0
0
1
1 2
1 2
1 2
1 2
{z u
Proof. obtain
1 2 1 2
}
|
1 2 1 2
0
1 2
1 2
{z
0
(22) }
w
We first apply the SPRK method (22) to the spatial discretization for u and w respectively and Um,1 = um , U m,2 = um+1 , h Wm,1 = wm + ∂ xx um , 2 h Wm,2 = wm + ∂ xx um , 2 h um+1 = um + (Wm,1 + Wm,2 ), 2 h w m+1 = wm + (∂ xx U m,1 + ∂ xx U m,2 ). 2
(23)
W. Shi, X. Wu, and J. Xia
9
Then proceed with dum+1 ∧ dwm+1 =
d(um + hwm +
=
dum ∧ dwm +
=
dum ∧ dwm +
=
dum ∧ dwm +
=
dum ∧ dwm +
h2 h h ∂ xx um ) ∧ d(wm + ∂ xx um + ∂ xx um+1 ) 2 2 2 h h dum ∧ ∂ xx dum + dum+1 ∧ ∂ xx dum+1 2 2 h h dUm,1 ∧ ∂ xx dUm,1 + dUm,2 ∧ ∂ xx dUm,2 2 2 h h dUm,1 ∧ ∂ x dWm,1 + dUm,2 ∧ ∂ x dWm,2 2 2 2 hP (dUm,i ∧ ∂ x dWm,i ). 2 i=1
That is, dum+1 ∧ dwm+1 − dum ∧ dwm =
2 hP (dUm,i ∧ ∂ x dWm,i ). 2 i=1
(24)
For t = tn,k , we rewrite (24) as n,k n,k n,k dun,k m+1 ∧ dwm+1 − dum ∧ dwm =
2 hP n,k n,k (dUm,i ∧ ∂ x dWm,i ). 2 i=1
(25)
It follows from (17) and (25) that the multi-symplectic conservation law is
+
2 1 P n n (dvn+1 ∧ dun+1 m,i − dvm,i ∧ dum,i )h 2 m,i i=1 µ ¶ r ¡ ¢ P 1 n,k n,k n,k bk (νt )φ0 (νt ) + b¯ k (νt )νt2 φ1 (νt ) (dun,k m+1 ∧ dwm+1 − dum ∧ dwm )τ = 0. k=1 φ0 (ck νt )
(26)
¤ Similarly, we can also obtain the following fully-discrete scheme for (8) h2 n,k n,k un,k ∂ xx un,k = u + hw + m , m m m+1 2 h n,k n,k n,k wm+1 = wn,k m + (∂ xx um + ∂ xx um+1 ), 2 r P k = 1, . . . , r, Umn,k = φ0 (ck νt )unm + ck φ1 (ck νt )τvnm + τ2 a¯ kl (νt )gn,l m, l=1 r n n 2P ¯ un+1 bk (νt )gn,k m = φ0 (νt )um + φ1 (νt )τvm + τ m , k=1 r vn+1 = −νt ωt φ1 (νt )un + φ0 (νt )vn + τ P bk (νt )gn,k . m m m m
(27)
k=1
4. Construction of explicit multi-symplectic schemes Based on the general framework presented in Section 3, a large number of novel multi-symplectic schemes can be constructed. Here, an undoubted fact is that, symplectic Runge-Kutta schemes for the integration of general Hamiltonian systems are usually implicit, and multi-symplectic integrators constructed by concatenating symplectic Runge-Kutta methods and/or symplectic partitioned Runge-Kutta methods are usually implicit. Thus, in practice, one has to solve implicit algebraic equations, often by using some iterative approximation methods. In this case, the resulting integration scheme may no longer be symplectic. The advantage of explicit symplectic integrators, when compared with fully implicit or partly implicit methods, is that they do not require the solution of large and complicated systems of nonlinear algebraic or transcendental equations. Therefore, it is clear that the computational cost of an implicit multi-symplectic method is usually higher than that of an explicit one. Consequently, in this section, we construct two explicit multi-symplectic schemes, which are actually two multi-symplectic extended leap-frog methods. It can be observed that when ω x → 0, ωt → 0, these two new schemes Eleap-frogI and Eleap-frogII reduce to classical leap-frog schemes.
W. Shi, X. Wu, and J. Xia
10
4.1. Extended leap-frog scheme I: An explicit multi-symplectic ERKN scheme In this subsection, we pay attention to constructing an explicit 2-stage MSERKN scheme of order 2. n,k n+1 First of all, since ∂ xx Um,i appear in (19), we observe that we cannot obtain un+1 m,i and vm,i explicitly provided n,k that ∂ xx Um,i could not be expressed explicitly. This means that, even if we use explicit SERKN methods for both spatial and temporal discretizations, the resulting multi-symplectic scheme may not be necessarily explicit. We take the spatial discretization as an example. If we chose c˜ = [0, 1]T , bˆ 1 (ν x ) = aˆ 21 (ν x ), and n,k n,k n,k bˆ 2 (ν x ) = aˆ 22 (ν x ) = 0, then Um,1 = un,k m , U m,2 = um+1 , and the second and third equations in (19) become n,k 2 n,k n,k n,k n,k 2ˆ um+1 =φ0 (ν x )um + φ1 (ν x )hwm + h b1 (ν x )(ω x um + ∂ xx um ), n,k n,k 2 n,k n,k ˜ wm+1 = − ν x ω x φ1 (ν x )un,k m + φ0 (ν x )wm + hb1 (ν x )(ω x um + ∂ xx um ) n,k + hb˜ 2 (ν x )(ω2x un,k m+1 + ∂ xx um+1 ). From this we derive that ∂ xx un,k m = and then ∂ xx un,k m+1 =
n,k n,k un,k m+1 − φ0 (ν x )um − φ1 (ν x )hwm − ω2x un,k m , h2 bˆ 1 (ν x )
n,k n,k un,k m+2 − φ0 (ν x )um+1 − φ1 (ν x )hwm+1 − ω2x un,k m+1 . h2 bˆ 1 (ν x )
Then it follows that n,k n,k 2 n,k n,k ˜ wn,k m+1 = − ν x ω x φ1 (ν x )um + φ0 (ν x )wm + hb1 (ν x )(ω x um + ∂ xx um ) + hb˜ 2 (ν x )(ω2x un,k + ∂ xx un,k ) m+1
m+1
b˜ 1 (ν x ) n,k b˜ 1 (ν x ) b˜ 1 (ν x ) um+1 − φ0 (ν x )un,k φ1 (ν x )wn,k m − ˆ m ˆ ˆ hb1 (ν x ) hb1 (ν x ) b1 (ν x ) b˜ 2 (ν x ) n,k b˜ 2 (ν x ) b˜ 2 (ν x ) + um+2 − φ0 (ν x )un,k − φ1 (ν x )wn,k m+1 m+1 . hbˆ 1 (ν x ) hbˆ 1 (ν x ) bˆ 1 (ν x )
n,k = − ν x ω x φ1 (ν x )un,k m + φ0 (ν x )wm +
Assuming
b˜ 2 (ν x )un,k m+2
bˆ 1 φ0 (ν x ) b˜ 1 (ν x ) = , φ1 (ν x ) ¡ ¢ ¡ ¢ n,k ˜ ˜ − bˆ 1 (ν x )ν2x φ1 (ν x ) + b˜ 1 (ν x )φ0 (ν x ) un,k m + b1 (ν x ) − b2 (ν x )φ0 (ν x ) u
(28)
m+1 . hbˆ 1 (ν x ) + hb˜ 2 (ν x )φ1 (ν x ) Therefore, from (28), together with the symplectic conditions (13), (14) and order conditions of ERKN methods proposed in [18], we can construct several 2-stage SERKN methods of order 2. In this paper, we consider the following 2-stage SERKN methods of order 2 for the two directional discretizations, which can be expressed in Butcher tableaus:
then we have wn,k m+1 =
0 1
0 φ2 (νt )
0 0
¯ t) b(ν
φ2 (νt ) φ0 (νt )φ2 (νt ) φ1 (νt )
0 φ0 (νt )φ2 (νt ) φ1 (νt ) − φ1 (νt )
0 1
0 φ2 (ν x )
0 0
ˆ x) b(ν ˜ x) b(ν
φ2 (ν x ) φ0 (ν x )φ2 (ν x ) φ1 (ν x )
0 φ0 (ν x )φ2 (ν x ) φ1 (ν x ) − φ1 (ν x )
b(νt )
(29)
(30)
11
W. Shi, X. Wu, and J. Xia
1 − φ0 (ν) where φ2 (ν) = . ν2 Note that as ωt → 0, ω x → 0, the above two SERKN methods (29) and (30) for time and space discretizations, respectively, reduce to the same classical 2-stage symplectic RKN (SRKN) method of order 2: 0
0
0
1
1 2
0
1 2
0
1 2
1 2
This coincides with the leap-frog method qn+1 − 2qn + qn−1 = f (qn ) h2 for problem p˙ = f (q), q˙ = p. The leap-frog method can be translated from the following 1-stage SRKN method of order 2 by exchanging the roles of p and q: 1 2
0 (31)
1 2
1 With the discretizations of the SERKN methods (29) and (30), (19) becomes the following scheme Umn,1 = unm , n,1 n,1 Um+1 − 2Umn,1 + Um−1 n,1 , (∂ U) = xx m 2h2 φ2 (ν x ) ¡ 2 n,1 ¢ 0 n,2 n n 2 n,1 n,1 U = φ (ν )u + φ (ν )τv + τ φ (ν ) ω U + ∂ U − V (U ) , 0 t 1 t 2 t xx m m m t m m m n,2 U n,2 − 2U n,2 + Um−1 n,2 (32) (∂ xx U)m = m+1 2 m , 2h φ2 (ν x ) ¡ 2 n,1 ¢ n+1 0 n n 2 n,1 n,1 um = φ0 (νt )um + φ1 (νt )τvm + τ φ2 (νt ) ωt Um + ∂ xx Um − V (Um ) , n+1 vm = − νt¡ωt φ1 (νt¡)unm + φ0 (νt )vnm ¢ ¡ ¢¢ 0 0 + τ b1 (νt ) ω2t Umn,1 + ∂ xx Umn,1 − V (Umn,1 ) + b2 (νt ) ω2t Umn,2 + ∂ xx Umn,2 − V (Umn,2 ) . In fact, the scheme (32) gives the discretizations for (∂ xx u)nm and (∂tt u)nm by (∂ xx u)nm = and (∂tt u)nm =
unm+1 − 2unm + unm−1 2h2 φ2 (ν x )
n n−1 un+1 m − 2um + um , respectively. This also deduces an extended leap-frog scheme: 2 2τ φ2 (νt ) n n−1 un+1 un − 2un + unm−1 0 m − 2um + um = m+1 2 m − V (unm ). 2 2τ φ2 (νt ) 2h φ2 (ν x )
(33)
We refer to the scheme (33) as Eleap-frogI and use this scheme to solve a linear wave equation (Problem 1) in the next section. Meanwhile, the discretizations of w and v are given by wnm+1 =
n un+2 unm+2 − unm m − um , vn+1 . m = 2hφ1 (ν x ) 2τφ1 (νt )
Furthermore, we have the following theorem.
W. Shi, X. Wu, and J. Xia
12
Theorem 4.1. By simply exchanging the roles of u and w in space and u and v in time, the MSERKN scheme given by the following 1-stage SERKN methods of order 2 for the two directional discretizations can be transformed into the MSERKN scheme (33): 1 2
¯ t) b(ν
φ21 ( ν2t )
0 ± 2φ0 ( ν2t )
b(νt )
φ1 ( ν2t )
1 2
0 ± φ21 ( ν2x ) 2φ0 ( ν2x )
ˆ x) b(ν ˜ x) b(ν
φ1 ( ν2x )
(34)
(35)
Proof. Firstly, the symplecticity of the schemes (34) and (35) can be directly obtained from Theorem 3.1. Then, we only show the details for the space discretization, and those for the time discretization are similar. Applying the 1-stage SERKN method (35) with the spatial discretization yields 1 νx νx U1 = φ0 ( 2 )um + 2 φ1 ( 2 )hwm , (36) ˆ x ) f1 , um+1 = φ0 (ν x )um + φ1 (ν x )hwm + h2 b(ν w ˜ m+1 = −ν x ω x φ1 (ν x )um + φ0 (ν x )wm + hb(ν x ) f1 , where f1 = f (U1 ) = ω2x U1 + ∂ xx U1 . From the last two equations in (36), we get ˆ x ) f1 um+1 =φ0 (ν x )um + φ1 (ν x )hwm + h2 b(ν ˆ x ) wm+1 + ν x ω x φ1 (ν x )um − φ0 (ν x )wm =φ0 (ν x )um + φ1 (ν x )hwm + h2 b(ν ˜ x) hb(ν ˆ x) ˆ x) ˆ x) b(ν b(ν b(ν =φ0 (ν x )um + φ1 (ν x )hwm + h wm+1 + ν2x φ1 (ν x )um − h φ (ν )w ˜ x) ˜ x) ˜ x) 0 x m b(ν b(ν b(ν φ1 (ν x ) =um + h (wm + wm+1 ). 1 + φ0 (ν x ) This shows um+1 − um = h
φ1 (ν x ) (wm + wm+1 ). 1 + φ0 (ν x )
(37)
The first and last equations of (36) imply ˜ x ) f (U1 ) wm+1 = − ν x ω x φ1 (ν x )um + φ0 (ν x )wm + hb(ν ˜ x )(ω2x U1 + ∂ xx U1 ) = − ν x ω x φ1 (ν x )um + φ0 (ν x )wm + hb(ν µ ¶ νx νx (38) 2 21 ˜ = − ν x ω x φ1 (ν x )um + φ0 (ν x )wm + hb(ν x ) ω x φ0 ( )um + ω x φ1 ( )hwm + ∂ xx U1 2 2 2 νx =wm + hφ1 ( )∂ xx U1 . 2 Following the approach to obtaining the leap-frog method by exchanging the roles of p and q in the SRKN method (31), we also exchange the roles of u and w. Then it can be deduced from the first equation of (36), (37) and (38) that 1 νx νx W1 = φ0 ( )wm + φ1 ( )h∂ xx um , 2 2 2 φ1 (ν x ) (39) wm+1 = wm + h (∂ xx um + ∂ xx um+1 ), 1 + φ0 (ν x ) νx u )W1 . m+1 = um + hφ1 ( 2
13
W. Shi, X. Wu, and J. Xia
The first equation of (39) implies wm+1/2
µ ¶ νx 1 νx νx φ1 (ν x ) = φ0 ( )wm + φ1 ( )h∂ xx um = φ0 ( ) wm + h ∂ xx um , 2 2 2 2 1 + φ0 (ν x )
(40)
1 φ1 (ν x ) ∂ xx um . ν x wm+1/2 − h φ0 ( 2 ) 1 + φ0 (ν x )
(41)
and then wm =
Inserting (41) into the second equation of (39) gives wm+1 = and furthermore, wm = That is, wm−1/2
1 φ1 (ν x ) ∂ xx um+1 , wm+1/2 + h φ0 ( ν2x ) 1 + φ0 (ν x ) φ1 (ν x ) 1 wm−1/2 + h ∂ xx um . φ0 ( ν2x ) 1 + φ0 (ν x )
µ ¶ φ1 (ν x ) νx 1 νx νx ∂ xx um = φ0 ( )wm − φ1 ( )h∂ xx um . = φ0 ( ) wm − h 2 1 + φ0 (ν x ) 2 2 2
Therefore, it follows from (40) and (42) that νx wm+1/2 + wm−1/2 = 2φ0 ( )wm , 2 ν wm+1/2 − wm−1/2 = φ1 ( x )h∂ xx um . 2 On the other hand, the last equation of (39) means wm+1/2 =
1 2
νx u − um − (um − um−1 ) wm+1/2 − wm−1/2 = m+1 = φ1 ( )h∂ xx um . hφ1 ( ν2x ) 2
(44)
um+1 − um−1 2hφ1 (ν x ) u m+1 − 2um + um−1 . (∂ xx u)m = 2h2 φ2 (ν x ) wm =
For the time discretization, it also can be proved that (∂tt u)n = holds.
(43)
um+1 − um , which leads to hφ1 ( ν2x )
νx um+1 − um + (um − um−1 ) = 2φ0 ( )wm , wm+1/2 + wm−1/2 = hφ ( νx ) 2
It follows from (44) that
(42)
n un+1 − 2un + un−1 un+2 m − um n+1 and v = m 2τ2 φ2 (νt ) 2τφ1 (νt ) ¤
4.2. Extended leap-frog scheme II: An explicit multi-symplectic ERKN-PRK scheme In this subsection, we construct another explicit multi-symplectic scheme based on Theorem 3.3. We apply the SPRK method (22) to the space discretization, which is equivalent to the St¨ormer-Verlet that gives the second order central difference of ∂ xx Umn,k , and apply the explicit 2-stage SERKN method (29) to the time
14
W. Shi, X. Wu, and J. Xia
discretization. Then the multi-symplectic scheme (27) becomes n,1 U = unm , m n,1 n,1 Um+1 − 2Umn,1 + Um−1 n,1 (∂ U) = , xx m h2 ¡ ¢ 0 Umn,2 = φ0 (νt )unm + φ1 (νt )τvnm + τ2 φ2 (νt ) ω2t Umn,1 + ∂ xx Umn,1 − V (Umn,1 ) , n,2 U n,2 − 2Umn,2 + Um−1 n,2 (∂ xx U)m = m+1 , h2 ¡ 2 n,1 ¢ 0 n+1 n n 2 n,1 n,1 um = φ0 (νt )um + φ1 (νt )τvm + τ φ2 (νt ) ωt Um + ∂ xx Um − V (Um ) , vn+1 = − νt¡ωt φ1 (νt¡)unm + φ0 (νt )vnm ¢ ¡ ¢¢ m 0 0 + τ b1 (νt ) ω2t Umn,1 + ∂ xx Umn,1 − V (Umn,1 ) + b2 (νt ) ω2t Umn,2 + ∂ xx Umn,2 − V (Umn,2 ) .
(45)
Just like the situation of the scheme Eleap-frogI, it can be observed that this scheme is equivalent to the following form n n−1 un+1 unm+1 − 2unm + unm−1 0 m − 2um + um = − V (unm ), (46) 2 2 2τ φ2 (νt ) h and the discretizations of w and v are wnm+1 =
n unm+2 − unm un+2 m − um , vn+1 . m = 2h 2τφ1 (νt )
It is also equivalent to the concatenation of the 1-stage SERKN method (34) of order 2 and SPRK method (22) applied to the temporal and spatial discretizations, respectively. We refer to the above scheme (46) as Eleap-frogII, and in the next section, we use this scheme to solve a nonlinear wave equation (Problem 2). 4.3. Analysis of linear stability In this subsection, we pay attention to another significant aspect for those schemes proposed in this paper, namely, the stability. It is known that the attractiveness of explicit schemes relies on the computational efficiency, but such efficiency is always at the cost of stability. Although an application of the linear stability analysis to nonlinear equations cannot be justified, it is found to be effective in numerical experiments. Using the well-known Fourier method [40], we give some linear stability analysis for the explicit multi-symplectic Eleap-frogI scheme, and the stability for scheme Eleap-frogII can be analyzed in the similar way. Set unm = uˆ n exp{iξmh}, vnm = vˆ n exp{iξmh}, Umn,1 = Uˆ n,1 exp{iξmh}, Umn,2 = Uˆ n,2 exp{iξmh}, etc., √ 0 where ξ ∈ R and i = −1. Then, substituting the above expressions into scheme (32) with V (Umn,i ) = Umn,i , we get n,1 n Uˆ = uˆ , ¶ µ iξh −iξh Uˆ n,2 = φ (ν )ˆun + φ (ν )τˆvn + τ2 φ (ν ) ω2 + (e − 2 + e ) − 1 Uˆ n,1 , 0 t 1 t 2 t t 2h2 φ2 (ν x ) n+1 n,2 uˆ = Uˆ , n+1 = −νt ωt φ1 (νt )ˆun + φ0 (νt )ˆvn vˆ µ ¶ ³ ´ ³ ´ (eiξh − 2 + e−iξh ) (eiξh − 2 + e−iξh ) 2 n,1 2 n,2 ˆ ˆ − 1 U + b2 (νt ) ωt + −1 U . +τ b1 (νt ) ωt + 2h2 φ2 (ν x ) 2h2 φ2 (ν x ) Eliminating Uˆ n,1 , Uˆ n,2 , we obtain
µ
uˆ n+1 vˆ n+1
¶
µ =G
uˆ n vˆ n
¶ ,
15
W. Shi, X. Wu, and J. Xia
where
τ2 φ2 (νt )θ 2 1 − τ φ2 (νt ) + h2 φ2 (ν x ) G= G21
φ1 (νt )τ
, 2 τ φ (ν )θ 2 t 1 − τ2 φ2 (νt ) + 2 h φ2 (ν x )
with θ = (cos ξh − 1), and G21 =
τφ1 (νt )θ + τ3 b2 (νt )νt2 φ2 (νt )θ(ω2t − 2) τ3 b2 φ2 (νt )θ2 − φ1 (νt )τ + b2 (νt )φ2 (νt )τ3 (1 − ω2t ). + h2 φ2 (ν x ) h4 φ22 (ν x )
We can easily verify that G11G22 − G12G21 = 1. Hence, the eigenvalues of G = G(ξ, θ, ν x , νt ) satisfy λ2 − (G11 + G22 )λ + (G11G22 − G12G21 ) = 0. That is, λ2 − (G11 + G22 )λ + 1 = 0.
(47)
Obviously, the roots λi (i = 1, 2) of equation (47) satisfy |λi | ≤ 1 if and only if |G11 + G22 | ≤ 2. Then the method satisfies the von Neumann condition, which is necessary for the stability [40]. Assuming |G11 + G22 | ≤ 2 that ¯ ¯ 2 ¯ ¯ ¯1 − τ2 φ2 (νt ) + τ φ2 (νt )θ ¯ ≤ 1, ¯ 2 h φ2 (ν x ) ¯ which is equivalent to
τ2 φ2 (νt )θ ≤ 1. (48) h2 φ2 (ν x ) τ π Denote r = , and suppose 0 < r ≤ α and 0 < νt ≤ ν x ≤ . Since φ2 (ν) is a monotonically decreasing h 2 π 1 function on the interval [0, ] with φ2 (0) = and −2 ≤ θ ≤ 0, we obtain 2 2 −1 ≤ 1 − τ2 φ2 (νt ) +
1 1/2 τ2 φ2 (νt )θ − 2α2 < 1 − τ2 φ2 (νt ) + 2 < 1. 2 φ2 (π/2) h φ2 (ν x ) √ 6 1 1/2 Obviously, if 0 < α ≤ , then − 2α2 ≥ −1, and the inequality (48) is satisfied. π 2 φ2 (π/2) Furthermore, it also gives a sufficient condition for the stability of Eleap-frogI. In fact, conditions 0 < √ 6 π τ2 φ2 (νt )θ α≤ and νt ≤ ν x ≤ make −1 < 1 − τ2 φ2 (νt ) + 2 < 1, which implies that the equation (47) has π 2 h φ2 (ν x ) two distinct conjugate roots λi with |λi | = 1, i = 1, 2. Therefore, the method is linearly stable if √ √ 6 6 π 0 0 is just the initial condition translated through space: u(x, t) = u(x − c0 t, 0). ∂D ± ∂D 0 . An important quantity of dispersive systems is the group velocity defined by, vg (k) = ω (k) = − ∂k ∂ω The group velocity is usually a function of the wave number k and characterizes the speed of energy transport 00 in wave packets. In this case, ω (k) , 0, and waves with different wave numbers travel with different phase velocities v p , which brings the group velocity dispersion (GVD) and ultimately causes the spatial spreading of the wave packet. The dispersive properties of the discretizations can be studied by calculating their numerical dispersion relations. To do so, we start with the discretizations of the linearized sine-Gordon equation (55) associated with the Eleap-frogI and Eleap-frogII schemes in their single equation format D2 un−1 D2 un Eleap − frogI : t m − x m−1 + χunm = 0, 2φ2 (νt ) 2φ2 (ν x ) (57) D2 un−1 Eleap − frogII : t m − D2x unm−1 + χunm = 0, 2φ2 (νt )
22
W. Shi, X. Wu, and J. Xia
n un+1 un − unm m − um and D x unm = m+1 are the finite difference operators. τ h P i(Km−Ωn) Assume a discrete general solution to a linear PDE of the form Ae , K = kh, Ω = ωτ, where
where Dt unm =
K
each mode is a solution provided K and Ω satisfy a numerical dispersion relation DN (Ω, K) = 0. For the equation (57), the corresponding numerical dispersion relations are given by µ√ ¶2 µ √ ¶2 2φ2 (νt ) 2φ2 (ν x ) Ω K sin − sin − χ = 0, Eleap − frogI : φ (ν )τ 2 φ (ν )h 2 µ √2 t ¶2 µ 2 x ¶2 2φ2 (νt ) Ω 2 K sin − sin − χ = 0. Eleap − frogII : φ2 (νt )τ 2 h 2
(58)
Both these two schemes preserve the form of the analytical dispersion relation (56). The exact relationship is given by the following theorem. Theorem 5.1. The multi-symplectic schemes (58) qualitatively preserves the dispersion relation of the linear PDE. Specifically, there exist diffeomorphisms ψ1 and ψ2 satisfying the exact dispersion relationship ¡ ¢ DN (Ω, K) = D ψ1 (Ω), ψ2 (K) = ψ21 − ψ22 − χ = 0, (59) where
¶ µ√ √ ¡ ¢ 2φ2 (νt ) Ω 2φ2 (ν x ) K sin , sin , Eleap − frogI : ψ1 (Ω), ψ2 (K) = φ2 (νt )τ 2 φ2 (ν x )h 2 ¶ µ√ ¡ ¢ 2φ2 (νt ) Ω 2 K sin , sin , Eleap − frogII : ψ1 (Ω), ψ2 (K) = φ2 (νt )τ 2 h 2
(60)
for −π < K, Ω < π. Solving the numerical dispersion relations (58) gives s φ2 (νt )r2 K φ2 (νt )hr Eleap − frogI : Ω(K) = ±2 arcsin sin2 + χ, φ2 (ν x ) 2 2 r K φ (ν )hr Eleap − frogII : Ω(K) = ±2 arcsin 2φ2 (νt )r2 sin2 + 2 t χ. 2 2 The significance of Theorem 5.1 for the non-dispersive equation (55) (with χ = 0) can be seen in Figures 0 00 9–11, which show the dispersion curves Ω(K), the group velocity Ω (K), the group velocity dispersion Ω (K) for the exact (56), and numerical (58). In the experiment, two different values of the ratio r = τ/h, with h = 0.01 are used. The exact relation is given by Ω(K) = ±rK. The significance of Theorem 5.1 for the linearized sine-Gordon equation (55) (with χ = 1) can be seen in √ Figures 12–14. The exact relation is given by Ω(K) = ±r K 2 + h2 . From Figures 9–14, we can see that the existence of diffeomorphisms (60) is not enough to preserve the qualitative features of the analytical solution. All of the schemes introduce numerical dispersions. From Figures 10 and 13, because Ω(K) and K both have direction, we can see that the dispersion relation curves are all monotonically increasing at rates given by their numerical group velocities. Similarly, Figures 11 and 14 display that the group velocity curves (see Figures 10 and 13 ) are monotonically decreasing. Therefore, considering the numerical solutions of the linear wave equations given by Eleap-frogI and Eleap-frogII schemes, we further obtain the following corollaries. ∂DN ± ∂DN 0 at (K, Ω) of the two schemes Corollary 5.1. The numerical group velocities Vg (k) = Ω (K) = − ∂K ∂Ω ∂D ± ∂D 0 Eleap-frogI and Eleap-frogII have the same sign as the group velocity vg (k) = ω (k) = − for the ∂k ∂ω associated pair (k, ω). Corollary 5.2. Both Eleap-frogI and Eleap-frogII schemes introduce negative dispersion, and when the wave number is increasing, the numerical group velocities are monotonically decreasing.
23
W. Shi, X. Wu, and J. Xia
r=0.7
r=0.4 1.5 Eleap−frogI Ω⋅K>0 Eleap−frogI Ω⋅K0 exact Ω⋅K0 Eleap−frogI Ω⋅K0 exact Ω⋅K0 Eleap−frogI Ω⋅K0 exact Ω⋅K0 Eleap−frogI Ω⋅K0 exact Ω⋅K0
’
0.3
Eleap−frogI Ω (K)⋅K>0
0.2
Eleap−frogI Ω (K)⋅K0 Eleap−frogII Ω⋅K0 exact Ω⋅K0 Eleap−frogII Ω⋅K0 exact Ω⋅K0 Eleap−frogII Ω⋅K0 exact Ω⋅K0 Eleap−frogII Ω⋅K0 exact Ω⋅K0
0.4
’
Eleap−frogII Ω (K)⋅K>0
Eleap−frogII Ω’(K)⋅K0
0.2
exact Ω (K)⋅K0
0.4
’
Ω (K)
’
Eleap−frogII Ω (K)⋅K
of Mathematics, Nanjing University; State Key Laboratory for Novel Software Technology at Nanjing University, Nanjing 210093, P.R.China b Department of Mathematics, Purdue University, West Lafayette, IN 47907, USA
Abstract In this paper, we study the integration of Hamiltonian wave equations whose solutions have oscillatory behaviors in time and/or space. We are mainly concerned with the research for multi-symplectic extended Runge-Kutta-Nystr¨om (ERKN) discretizations and the corresponding discrete conservation laws. We first show that the discretizations to the Hamiltonian wave equations using two symplectic ERKN methods in space and time respectively lead to an explicit multi-symplectic integrator (Eleap-frogI). Then we derive another multi-symplectic discretization using a symplectic ERKN method in time and a symplectic Partitioned Runge-Kutta method, which is equivalent to the well-known St¨ormer-Verlet method in space (Eleap-frogII). These two new multi-symplectic schemes are extensions of the leap-frog method. The numerical stability and dispersive properties of the new schemes are analyzed. Numerical experiments with comparisons are presented, where the two new explicit multi-symplectic methods and the leap-frog method are applied to the linear wave equation and the sine-Gordon equation. The numerical results confirm the superior performance and some significant advantages of our new integrators in the sense of structure preservation. Keywords: Multi-symplectic integrators, Hamiltonian wave equations, extended leap-frog methods, conservation laws, dispersive properties
Mathematics Subject Classification (2000): 35L05, 37M15, 65L06, 65M12, 65P10, 70S10
1. Introduction It is well-known that symplectic integrators are robust, efficient, and accurate in preserving the longtime behavior of the solutions of Hamiltonian ordinary differential equations (ODEs) [1]. The basic idea of a symplectic integrator is that the numerical scheme is designed to preserve the symplectic form at each time step. Recently, it is shown that many conservative partial differential equations (PDEs), such as wave equations, generalized KdV and nonlinear Schr¨odinger equations, etc., allow for a structure similar to the symplectic structure of Hamiltonian ODEs, called the multi-symplectic formulation (see, e.g., [2, 3, 4]). For example, in [4], Marsden and Shkoller develop the multi-symplectic structure of Hamiltonian PDEs from a Lagrangian formulation using a variational principle. It is natural to design integrators for Hamiltonian PDEs which preserve a semi-discretization of the symplectic form associated with the infinite-dimensional evolution equation. In the past decade, rapid progress I The research of Wei Shi and Xinyuan Wu was supported in part by the Natural Science Foundation of China under Grant 10771099 and by the Specialized Research Foundation for the Doctoral Program of Higher Education under Grant 20100091110033, by the 985 Project at Nanjing University under Grant 9112020301, by A Project Founded by the Priority Academic Program Development of Jiangsu Higher Education Institutions, and by the University Postgraduate Research and Innovation Project of Jiangsu Province 2012 under Grant CXLX12− 0033. II The research of Jianlin Xia was supported in part by NSF grants DMS-1115572 and CHE-0957024. ∗ Corresponding author. Email address: [email protected], [email protected], [email protected] (Jianlin Xia)
1
W. Shi, X. Wu, and J. Xia
2
has been made on the development of multi-symplectic integrators of Runge-Kutta type for various Hamiltonian PDEs (see, e.g., [5, 6, 7, 8, 9, 10, 11, 12]). Unfortunately, most of the existing multi-symplectic integrators fail to take into account the oscillatory behavior of the problems of interest, which is important for the solution that has the intrinsic oscillatory behavior. It is known that most of integrators adapted to the oscillatory problems are frequency-dependent, such as ARKN methods [13, 14, 15, 16, 17], ERKN methods [18, 19, 20, 21], and trigonometrically/exponentiallyfitted methods [22, 23, 24, 25]. For nonlinear Schr¨odinger equations, taking the frequency of the problem into account, several useful trigonometrically/exponentially-fitted methods are proposed (see, e.g., [26, 27, 28, 29, 30, 31, 32]). Some authors, such as Gonz´alez, et al. [33] and Franco [13], are interested in numerical integrators of Runge-Kutta-Nystr¨om type adapted to initial value problems of the form ¡ ¢ ½ 00 y (t) + ω2 y(t) = f y(t) , t ∈ [t0 , T ], (1) y(t0 ) = y0 , y0 (t0 ) = y00 , where ω is the main frequency and ω > 0. More recently, a new family of extended Runge-Kutta-Nystr¨om (ERKN) methods are proposed by Yang, et al. [18] for one-dimensional perturbed oscillators and systems of these oscillators with ω2 being a diagonal matrix. Furthermore, with the methodology in [18], a standard form of multidimensional ERKN integrators for the general system (1) is formulated by Wu, et al. [19], where ω2 is a symmetric positive semi-definite matrix and implicitly contains the frequencies of the problem. The corresponding order conditions are also derived by the authors based on the B-series theory associated with the extended Nystr¨om trees (EN-trees) [19]. These methods preserve the oscillatory feature of the unperturbed oscillators, not only for the updates but also for the internal stages. Numerical results have shown that ERKN methods are superior to other methods for oscillatory systems. In this paper, we investigate the multi-symplectic ERKN (MSERKN) methods for wave equations as Hamiltonian PDEs. In Section 2, we discuss the conservation laws for Hamiltonian wave equations. In Section 3, we show that the discretization in time and space for Hamiltonian wave equations with two symplectic ERKN (SERKN) methods leads to a multi-symplectic integrator. In Section 4, explicit multi-symplectic schemes are constructed based on the SERKN methods and a symplectic Partitioned Runge-Kutta (SPRK) method. Section 5 is devoted to the numerical experiments, including the dispersive properties of the schemes proposed in this paper. Some conclusions are given in the last section. 2. Conservation laws and multi-symplectic structures of wave equations 2.1. Multi-symplectic conservation laws and multi-symplectic integrators A multi-Hamiltonian PDE in one time dimension and one space dimension is a PDE that can be formulated in the following way (see [34]). Let M and K be any skew-symmetric real n × n matrices (n ≥ 3) and let S : Rn → R be any smooth function of the state variable z. Then a system of the following form is called a Hamiltonian system on a multi-symplectic structure: Mzt + Kz x = ∇z S (z), z ∈ Rn , (x, t) ∈ R2 ,
(2)
where ∇z is the gradient operator corresponding to the standard inner product on Rn denoted by h·, ·i. The Hamiltonian system (2) is multi-symplectic in the sense that it is associated with two-forms: ζ=
1 1 dz ∧ Mdz, κ = dz ∧ Kdz, 2 2
which define a space-time symplectic structure governed by the multi-symplectic conservation law ∂t ζ + ∂ x κ = 0. Let L be the length of the spatial domain. Then integrating (3) over the spatial domain leads to d dt
Z
L/2
ζdx + κ| x=L/2 − κ| x=−L/2 = 0. −L/2
(3)
W. Shi, X. Wu, and J. Xia
3
This equation, under appropriate boundary conditions such as z(L/2, t) = z(−L/2, t), results in the conservad tion of global symplecticity: ζ(z) = 0, where dt Z L/2 ζ(z) = ζ(z)dx. (4) −L/2
ζ(z) is a reminder indicating that the two-form ζ is integrated over all space. A properly chosen semidiscretization of this system may lead to a system of Hamiltonian ODEs with a symplectic form given by a discrete version of (4). An overview of this viewpoint in terms of symplectic time-integration methods is given by McLachlan [35]. Two other conservation laws for the energy and the momentum respectively follow from (3). The local energy conservation law is ∂t E(z) + ∂ x F(z) = 0, (5) 1 1 where E(z) = S (z)− zT Kz x is the energy density and F(z) = zT Kzt is the energy flux. The local momentum 2 2 conservation law is ∂t I(z) + ∂ xG(z) = 0, (6) 1 1 where I(z) = zT Mz x and G(z) = S (z) − zT Mzt are the momentum density and the momentum flux, 2 2 respectively. Similar to the conservation of multi-symplecticity, integrating the local conservation laws (5) and (6) over the spatial domain in Rn with periodic boundary conditions leads to global conservation laws d d ε(z) = 0 and I(z) = 0, (7) dt dt R L/2 R L/2 where ε(z) = −L/2 E(z)dx and I(z) = −L/2 I(z)dx. We note that the global conservation of energy or momentum (7) is necessary but not sufficient for the local conservation of energy or momentum. When it comes to the numerical integration, the symplecticity of numerical integrators should be determined. Definition 2.1. [34] A numerical scheme is called multi-symplectic if it preserves a discrete conservation law of (3). An alternative definition of multi-symplectic integrators based on a discrete form of the symplectic conservation law is suggested by Bridges and Reich [3]. It has been shown that popular methods such as the centered Preissman scheme and the leap-frog method are multi-symplectic and that such schemes have remarkable properties of conserving local energy and momentum. 2.2. Conservation laws for wave equations Consider the following nonlinear wave equation (see [36]): utt − u xx + V 0 (u) = 0,
(8)
where V(u) is certain smooth function and u is a scalar function. Introducing two new variables v = ut and w = u x , we can write (8) as vt − w x + V 0 (u) = 0, v = ut , w = ux . By taking z = (u, v, w)T ∈ R3 , we can obtain a first order PDE system of the abstract form (2) corresponding to (8) (see, e.g., [11]), with 0 0 1 0 −1 0 0 0 , K = 0 0 0 , M= 1 −1 0 0 0 0 0
4
W. Shi, X. Wu, and J. Xia
1 and the Hamiltonian is S (z) = (v2 − w2 ) + V(u). As mentioned in Subsection 2.1, this system has the 2 multi-symplectic conservation law (3), where ζ and κ are the pre-symplectic forms 1 dz ∧ Mdz = dv ∧ du, 2 1 κ = dz ∧ Kdz = du ∧ dw. 2
ζ=
(9) (10)
By applying the general approach of deriving the local energy and momentum conservation law to the nonlinear wave equation (8), as discussed in the previous subsection, it follows that the local energy conserva1 tion law is (5) with E(z) = (v2 + w2 ) + V(u) and F(z) = −vw. Meanwhile, the local momentum conservation 2 1 law is (6) with I(z) = vw and G(z) = − (v2 + w2 ) + V(u) (see [11]). 2 3. Extended RKN discretization of wave equations 3.1. Extended RKN methods for ODEs In [18], Yang, et al. investigate ERKN methods for the IVP (1) of the following form : s P Yi = φ0 (ci ν)yn + ci φ1 (ci ν)hy0n + h2 a¯ i j (ν) f (tn + c j h, Y j ), i = 1, . . . , s, j=1 s P yn+1 = φ0 (ν)yn + φ1 (ν)hy0n + h2 b¯ i (ν) f (tn + ci h, Yi ), i=1 s 0 0 y = −νwφ1 (ν)yn + φ0 (ν)y + h P bi (ν) f (tn + ci h, Yi ), n+1
n
(11)
i=1
where a¯ i j (ν), bi (ν) and b¯ i (ν), i, j = 1, . . . , s, are assumed to be even functions of ν = hω, φ0 (ν) = cos(ν), and sin(ν) φ1 (ν) = . We denote (11) by {ci , bi , b¯ i , a¯ i j , ω} in this paper. ν Very recently, ERKN methods are considered for Hamiltonian ODEs, whose solutions exhibit a periodic or oscillatory character. The coefficients of these methods depend on the product of the step size h and the fitted frequency ω, and the methods can be expressed in Butcher tableau as :
c
¯ A(ν) b¯ T (ν) bT (ν)
=
c1 .. .
a¯ 11 (ν) .. .
cs
a¯ 1s (ν) .. .
a¯ s1 (ν)
... .. . ···
b¯ 1 (ν)
···
b¯ s (ν)
b1 (ν)
···
b s (ν)
a¯ ss (ν)
The symplecticity conditions for ERKN methods are given by the next theorem. Theorem 3.1. (Wu, et al. [37, 20]) If the coefficients of an s-stage ERKN method satisfy the following conditions ¶ µ ¶ µ 2 b¯ i (ν) φ0 (ν) + ci ν φ1 (ν)φ1 (ci ν) = bi (ν) φ1 (ν) − ci φ0 (ν)φ1 (ci ν) , φ0 (ci ν) φ0 (ci ν) iµ = 1, . . . , s, ¶ a¯ i j (ν)φ0 (ν) ν2 φ1 (ν) (12) ¯ − b¯ i (ν)¯ai j (ν) b (ν) b (ν) − i j φ (c ν) φ (c ν) 0 i 0 i ¶ µ a¯ ji (ν)φ0 (ν) ν2 φ1 (ν) − b¯ j (ν)¯a ji (ν) , i, j = 1, . . . , s, = b j (ν) b¯ i (ν) − φ0 (c j ν) φ0 (c j ν) then the method is symplectic.
5
W. Shi, X. Wu, and J. Xia
R 3.1. When ω → 0, the scheme (11) reduces to a classical RKN method in Hairer, et al. [38]. In this case, the symplecticity conditions of (12) become ½ b¯ i = (1 − ci )bi , i = 1, . . . , s, bi (b¯ j − a¯ i j ) = b j (b¯ i − a¯ ji ), i ≤ j = 1, . . . , s, which are obtained by Suris [39] for the symplecticity of an s-stage RKN method. 3.2. Multi-symplectic extended RKN discretization In order to derive the new MSERKN methods, we first consider the discretization of the equation (8) by employing the following two different SERKN methods in temporal and spatial directions, respectively: c1 .. .
a¯ 11 (νt ) .. .
a¯ 1r (νt ) .. .
c˜ 1 .. .
aˆ 11 (ν x ) .. .
a¯ r1 (νt )
... .. . ···
cr
a¯ rr (νt )
c˜ s
aˆ s1 (ν x )
b¯ 1 (νt )
···
b¯ r (νt )
b1 (νt )
···
br (νt )
bˆ 1 (ν x ) b˜ 1 (ν x )
... .. . ···
aˆ 1s (ν x ) .. .
···
bˆ s (ν x ) b˜ s (ν x )
···
aˆ ss (ν x )
Theorem 3.2. Assume that a method results from the temporal discretization of the wave equation (8) with an r-stage SERKN scheme {ci , bi , b¯ i , a¯ i j , ωt } and spatial discretization with an s-stage SERKN scheme {˜ci , b˜ i , bˆ i , aˆ i j , ω x }, where the coefficients satisfy the symplectic conditions ¶ µ ¶ µ ci φ0 (νt )φ1 (ci νt ) ci νt 2 φ1 (νt )φ1 (ci νt ) ¯ = bi (νt ) φ1 (νt ) − , bi (νt ) φ0 (νt ) + φ0 (ci νt ) φ0 (ci νt ) iµ= 1, . . . , r, ¶ a¯ i j (νt )φ0 (νt ) νt 2 φ1 (νt ) ¯ (13) bi (νt ) b j (νt ) − − b¯ i (νt )¯ai j (νt ) φ0 (ci νt ) µ φ0 (ci νt ) ¶ a¯ ji (νt )φ0 (νt ) νt 2 φ1 (νt ) = b j (νt ) b¯ i (νt ) − − b¯ j (νt )¯a ji (νt ) , φ0 (c j νt ) φ0 (c j νt ) i, j = 1, . . . , r, and
¶ µ ¶ µ 2 ˆ i (ν x ) φ0 (ν x ) + c˜ i ν x φ1 (ν x )φ1 (˜ci ν x ) = b˜ i (ν x ) φ1 (ν x ) − c˜ i φ0 (ν x )φ1 (˜ci ν x ) , b φ0 (˜ci ν x ) φ0 (˜ci ν x ) iµ = 1, . . . , s, ¶ aˆ i j (ν x )φ0 (ν x ) ν x 2 φ1 (ν x ) b˜ i (ν x ) bˆ j (ν x ) − − bˆ i (ν x )ˆai j (ν x ) φ0 (˜ci ν x ) µ φ0 (˜ci ν x ) ¶ a ˆ (ν )φ (ν ) ν x 2 φ1 (ν x ) ji x 0 x = b˜ j (ν x ) bˆ i (ν x ) − − bˆ j (ν x )ˆa ji (ν x ) , φ0 (˜c j ν x ) φ0 (˜c j ν x ) i, j = 1, . . . , s,
(14)
respectively. The method is said to be multi-symplectic in the sense that it satisfies the discrete multisymplectic conservation law µ ¶ s ¢ P 1 ¡˜ 2 n+1 n n ˆ bi φ0 (ν x ) + bi ν x φ1 (ν x ) (dvn+1 m,i ∧ dum,i − dvm,i ∧ dum,i )h ci ν x ) i=1 µ φ0 (˜ ¶ (15) r ¢ P 1 ¡ n,k n,k n,k + bk φ0 (νt ) + b¯ k νt2 φ1 (νt ) (dun,k ∧ dw − du ∧ dw )τ = 0, m m m+1 m+1 k=1 φ0 (ck νt ) where τ and h are the step-size of time and space, respectively.
6
W. Shi, X. Wu, and J. Xia
Proof. (8) as
We first apply an s-stage SERKN method with the spatial discretization and rewrite equation s P Ui = φ0 (˜ci ν x )um + c˜ i φ1 (˜ci ν x )hwm + h2 aˆ i j (ν x ) f j , j=1 s P um+1 = φ0 (ν x )um + φ1 (ν x )hwm + h2 bˆ i (ν x ) fi , i=1 s P wm+1 = −ν x ω x φ1 (ν x )um + φ0 (ν x )wm + h b˜ i (ν x ) fi ,
i = 1, . . . , s,
i=1
where fi = f (Ui ) = ω2x Ui + ∂ xx Ui . Having φ0 (ν) = cos(ν) and φ1 (ν) = dum+1 ∧ dwm+1 = dum ∧ dwm + h s ¡ P
sin(ν) in mind, we proceed with ν
s ¡ ¢ P b˜ i (ν x )φ0 (ν x ) + bˆ i (ν x )ν2x φ1 (ν x ) dum ∧ d fi i=1
s ¢ P b˜ i (ν x )φ1 (ν x ) − bˆ i (ν x )φ0 (ν x ) dwm ∧ d fi + h3 bˆ i (ν x )b˜ j (ν x )d fi ∧ d f j i=1 i, j=1 ¶ µ s ¢ P 1 ¡˜ 2 ˆ bi (ν x )φ0 (ν x ) + bi (ν x )ν x φ1 (ν x ) dUi ∧ d fi = dum ∧ dwm + h φ0 (˜ci ν x ) i=1 µ ¶ s ¢ c˜ i φ1 (˜ci ν x ) P ¡ 2 2 ˜ ˆ ˜ ˆ −h − bi (ν x )φ1 (ν x )+ bi (ν x )φ0 (ν x ) dwm ∧d fi bi (ν x )φ0 (ν x )+ bi (ν x )ν x φ1 (ν x ) φ0 (˜ci ν x ) i=1 µ ¶ 2 s P b˜ i (ν x )φ0 (ν x ) + bˆ i (ν x )ν x φ1 (ν x ) 3 ˆ ˜ +h aˆ i j + bi (ν x )b j (ν x ) d fi ∧ d f j . φ0 (˜ci ν x ) i, j=1 2
+h
Together with (14), this is further transformed to ¶ µ s ¢ P 1 ¡˜ dum+1 ∧ dwm+1 = dum ∧ dwm + h bi (ν x )φ0 (ν x ) + bˆ i (ν x )ν2x φ1 (ν x ) dUi ∧ d fi ci ν x ) i=1 µ φ0 (˜ ¶ s ¢ P 1 ¡˜ = dum ∧ dwm + h bi (ν x )φ0 (ν x ) + bˆ i (ν x )ν2x φ1 (ν x ) dUi ∧ (ω2x dUi + ∂ xx dUi ). φ0 (˜ci ν x ) i=1 Then, we obtain the following semi-symplectic conservation law in space: ¶ µ s ¢ P 1 ¡˜ 2 ˆ bi (ν x )φ0 (ν x ) + bi (ν x )ν x φ1 (ν x ) dUi ∧ ∂ x dWi . dum+1 ∧ dwm+1 = dum ∧ dwm + h φ0 (˜ci ν x ) i=1
(16)
The next step is the discretization in time over a time interval [0, τ]. For simplicity, we set m = 0 and assume that xm = 0. We use the r-stage SERKN method and obtain r k 0 0 2P U = φ (c ν )u + c φ (c ν )τv + τ a¯ kl (νt )gli , k = 1, . . . , r, 0 k t k 1 k t i i i l=1 r P u1i = φ0 (νt )u0i + φ1 (vt )τv0i + τ2 b¯ k (νt )gki , k=1 r P v1i = −νt ωt φ1 (νt )u0i + φ0 (νt )v0i + τ bk (νt )gki . k=1
Here, we introduce the notation Uik ≈ u(˜ci h, ck τ), Vik ≈ v(˜ci h, ck τ), u1i ≈ u(˜ci h, τ), v1i ≈ v(˜ci h, τ), gki = g(Uik ) = ω2t Uik + ∂tt Uik . Likewise, we can obtain the identity du1i ∧ dv1i = du0i ∧ v0i + τ
r P k=1
µ
¶ ¢ 1 ¡ bk (νt )φ0 (νt ) + b¯ k (νt )νt2 φ1 (νt ) dUik ∧ ∂t dVik . φ0 (ck νt )
7
W. Shi, X. Wu, and J. Xia
This is the symplectic conservation law in time direction. For any x = xm,i = xm + c˜ i h, t = tn,k = tn + ck τ, the above identity can be transformed to ¶ r µ X ¢ 1 ¡ n,k n,k n+1 n+1 n n 2 ¯ dum,i ∧ dvm,i − dum,i ∧ dvm,i = τ bk (νt )φ0 (νt ) + bk (νt )νt φ1 (νt ) (dUm,i ∧ ∂t dVm,i ). (17) φ (c ν ) 0 k t k=1 Next, we rewrite (16) for t = tn,k as n,k n,k n,k dun,k m+1 ∧ dwm+1 − dum ∧ wm = h
s µ X i=1
¶ ¢ 1 ¡˜ n,k n,k bi (ν x )φ0 (ν x ) + bˆ i (ν x )ν2x φ1 (ν x ) (dUm,i ∧ ∂ x dWm,i ). (18) φ0 (˜ci ν x )
It finally follows from (17) and (18) that the discrete multi-symplectic conservation law (15) holds. Meanwhile, the fully-discrete scheme for (8) is given by s n,k n,k n,k n,k 2P U = φ (˜ c ν )u + c ˜ φ (˜ c ν )hw + h aˆ i j (ν x ) fm, 0 i x i 1 i x m m m,i j, j=1 s P n,k n,k um+1 = φ0 (ν x )un,k + φ1 (ν x )hwn,k + h2 bˆ i (ν x ) fm,i , m m i=1 s P n,k n,k n,k = −ν x ω x φ1 (ν x )un,k b˜ i (ν x ) fm,i , w m + φ0 (ν x )wm + h m+1
i = 1, . . . , s,
i=1
r P n,k Um,i = φ0 (ck νt )unm,i + ck φ1 (ck νt )τvnm,i + τ2 a¯ kl (νt )gn,l m,i , l=1 r n n 2P ¯ un+1 bk (νt )gn,k m,i = φ0 (νt )um,i + φ1 (νt )τvm,i + τ m,i , k=1 r P n n bk (νt )gn,k vn+1 m,i = −νt ωt φ1 (νt )um,i + φ0 (νt )vm,i + τ m,i ,
¤
(19) k = 1, . . . , r,
k=1
where we have used the following notation: n,k Um,i ≈ u(xm + c˜ i h, tn + ck τ),
un+1 ˜ i h, tn + τ), m,i ≈ u(xm + c
vn+1 ˜ i h, tn + τ), m,i ≈ v(xm + c
un,k m+1 ≈ u(xm + h, tn + ck τ),
wn,k m+1 ≈ w(xm + h, tn + ck τ),
n,k n,k n,k n,k 2 n,k 2 n,k 0 gn,k m,i = g(U m,i ) = ωt U m,i + ∂tt U m,i = ωt U m,i + ∂ xx U m,i − V (U m,i ), n,k n,k n,k n,k fm,i = f (Um,i ) = ω2x Um,i + ∂ xx Um,i .
The formula (15) can be understood as the approximation of the integral of (3) over the domain [xm , xm+1 ]× [tn , tn+1 ]. That is, (15) approximates Z xm+1 ¡ ¢ dv(x, tn+1 ) ∧ du(x, tn+1 ) − dv(x, tn ) ∧ du(x, tn ) dx Z
xm tn+1
+
¡
¢ du(xm+1 , t) ∧ dw(xm+1 , t) − du(xm , t) ∧ dw(xm , t) dt = 0,
(20)
tn
with SERKN methods for the evaluation of the two integrals. When the problems with periodic boundary conditions are considered, (15) has another interesting consequence. Let us take the sum of the identity (15) over all spatial grid points m = 1, . . . , M as: M nX s µ X
¶ ¢ 1 ¡˜ 2 n+1 n n ˆ bi (ν x )φ0 (ν x ) + bi (ν x )ν x φ1 (ν x ) (dvn+1 m,i ∧ dum,i − dvm,i ∧ dum,i )h φ (˜ c ν ) 0 i x m=1 i=1 ¶ r µ o X ¢ 1 ¡ n,k 2 n,k n,k ¯ + bk (νt )φ0 (νt ) + bk (νt )νt φ1 (νt ) (dun,k ∧ dw − du ∧ dw )τ = 0. m m m+1 m+1 φ0 (ck νt ) k=1
(21)
8
W. Shi, X. Wu, and J. Xia
Periodic boundary conditions in space imply M X r µ X
¶ ¢ 1 ¡ n,k n,k n,k bk (νt )φ0 (νt ) + b¯ k (νt )νt2 φ1 (νt ) (dun,k m+1 ∧ dwm+1 − dum ∧ dwm )τ φ (c ν ) 0 k t m=1 k=1 ¶ r µ X ¢ 1 ¡ n,k n,k n,k 2 ¯ = bk (νt )φ0 (νt ) + bk (νt )νt φ1 (νt ) (dun,k M+1 ∧ dw M+1 − du1 ∧ dw1 )τ = 0, φ (c ν ) 0 k t k=1 which in turn transforms (21) to M X s µ X
¶ ¢ 1 ¡˜ n+1 bi (ν x )φ0 (ν x ) + bˆ i (ν x )ν2x φ1 (ν x ) dvn+1 m,i ∧ dum,i φ (˜ c ν ) 0 i x m=1 i=1 ¶ M X s µ X ¢ 1 ¡˜ 2 ˆ = bi (ν x )φ0 (ν x ) + bi (ν x )ν x φ1 (ν x ) dvnm,i ∧ dunm,i . φ (˜ c ν ) 0 i x m=1 i=1 This is precisely the conservation of symplecticity in time. It is also the discretization of the total symplectic conservation (4) with (9). R 3.2. By abuse of notation, we denote the fitted temporal frequency and spatial frequency in the two SERKN methods {ci , bi , b¯ i , a¯ i j , ωt } and {˜ci , b˜ i , bˆ i , aˆ i j , ω x } by ωt and ω x , respectively. Obviously, when ωt → 0 and ω x → 0, the MSERKN methods become multi-symplectic RKN (MSRKN) methods. R 3.3. Theorem 3.2 shows that concatenation of SERKN methods can produce a multi-symplectic integrator for Hamiltonian wave equations. It is remarked that the above idea may apply equally to some other Hamiltonian PDEs. For Hamiltonian PDEs with oscillations in both space and time directions, two SERKN methods are used in space and time respectively. For the case that oscillation exists only in time direction, we propose the following multi-symplectic discretizations. Theorem 3.3. The wave equation (8) discretized by an SPRK method in space and an SERKN method in time has a discrete multi-symplectic conservation law. The proof is similar to that of Theorem 3.2. As an example, we use SERKN discretization in time and the following SPRK discretization in space (which can be regarded as the St¨ormer-Verlet formula):
|
0
0
0
1
1 2
1 2
1 2
1 2
{z u
Proof. obtain
1 2 1 2
}
|
1 2 1 2
0
1 2
1 2
{z
0
(22) }
w
We first apply the SPRK method (22) to the spatial discretization for u and w respectively and Um,1 = um , U m,2 = um+1 , h Wm,1 = wm + ∂ xx um , 2 h Wm,2 = wm + ∂ xx um , 2 h um+1 = um + (Wm,1 + Wm,2 ), 2 h w m+1 = wm + (∂ xx U m,1 + ∂ xx U m,2 ). 2
(23)
W. Shi, X. Wu, and J. Xia
9
Then proceed with dum+1 ∧ dwm+1 =
d(um + hwm +
=
dum ∧ dwm +
=
dum ∧ dwm +
=
dum ∧ dwm +
=
dum ∧ dwm +
h2 h h ∂ xx um ) ∧ d(wm + ∂ xx um + ∂ xx um+1 ) 2 2 2 h h dum ∧ ∂ xx dum + dum+1 ∧ ∂ xx dum+1 2 2 h h dUm,1 ∧ ∂ xx dUm,1 + dUm,2 ∧ ∂ xx dUm,2 2 2 h h dUm,1 ∧ ∂ x dWm,1 + dUm,2 ∧ ∂ x dWm,2 2 2 2 hP (dUm,i ∧ ∂ x dWm,i ). 2 i=1
That is, dum+1 ∧ dwm+1 − dum ∧ dwm =
2 hP (dUm,i ∧ ∂ x dWm,i ). 2 i=1
(24)
For t = tn,k , we rewrite (24) as n,k n,k n,k dun,k m+1 ∧ dwm+1 − dum ∧ dwm =
2 hP n,k n,k (dUm,i ∧ ∂ x dWm,i ). 2 i=1
(25)
It follows from (17) and (25) that the multi-symplectic conservation law is
+
2 1 P n n (dvn+1 ∧ dun+1 m,i − dvm,i ∧ dum,i )h 2 m,i i=1 µ ¶ r ¡ ¢ P 1 n,k n,k n,k bk (νt )φ0 (νt ) + b¯ k (νt )νt2 φ1 (νt ) (dun,k m+1 ∧ dwm+1 − dum ∧ dwm )τ = 0. k=1 φ0 (ck νt )
(26)
¤ Similarly, we can also obtain the following fully-discrete scheme for (8) h2 n,k n,k un,k ∂ xx un,k = u + hw + m , m m m+1 2 h n,k n,k n,k wm+1 = wn,k m + (∂ xx um + ∂ xx um+1 ), 2 r P k = 1, . . . , r, Umn,k = φ0 (ck νt )unm + ck φ1 (ck νt )τvnm + τ2 a¯ kl (νt )gn,l m, l=1 r n n 2P ¯ un+1 bk (νt )gn,k m = φ0 (νt )um + φ1 (νt )τvm + τ m , k=1 r vn+1 = −νt ωt φ1 (νt )un + φ0 (νt )vn + τ P bk (νt )gn,k . m m m m
(27)
k=1
4. Construction of explicit multi-symplectic schemes Based on the general framework presented in Section 3, a large number of novel multi-symplectic schemes can be constructed. Here, an undoubted fact is that, symplectic Runge-Kutta schemes for the integration of general Hamiltonian systems are usually implicit, and multi-symplectic integrators constructed by concatenating symplectic Runge-Kutta methods and/or symplectic partitioned Runge-Kutta methods are usually implicit. Thus, in practice, one has to solve implicit algebraic equations, often by using some iterative approximation methods. In this case, the resulting integration scheme may no longer be symplectic. The advantage of explicit symplectic integrators, when compared with fully implicit or partly implicit methods, is that they do not require the solution of large and complicated systems of nonlinear algebraic or transcendental equations. Therefore, it is clear that the computational cost of an implicit multi-symplectic method is usually higher than that of an explicit one. Consequently, in this section, we construct two explicit multi-symplectic schemes, which are actually two multi-symplectic extended leap-frog methods. It can be observed that when ω x → 0, ωt → 0, these two new schemes Eleap-frogI and Eleap-frogII reduce to classical leap-frog schemes.
W. Shi, X. Wu, and J. Xia
10
4.1. Extended leap-frog scheme I: An explicit multi-symplectic ERKN scheme In this subsection, we pay attention to constructing an explicit 2-stage MSERKN scheme of order 2. n,k n+1 First of all, since ∂ xx Um,i appear in (19), we observe that we cannot obtain un+1 m,i and vm,i explicitly provided n,k that ∂ xx Um,i could not be expressed explicitly. This means that, even if we use explicit SERKN methods for both spatial and temporal discretizations, the resulting multi-symplectic scheme may not be necessarily explicit. We take the spatial discretization as an example. If we chose c˜ = [0, 1]T , bˆ 1 (ν x ) = aˆ 21 (ν x ), and n,k n,k n,k bˆ 2 (ν x ) = aˆ 22 (ν x ) = 0, then Um,1 = un,k m , U m,2 = um+1 , and the second and third equations in (19) become n,k 2 n,k n,k n,k n,k 2ˆ um+1 =φ0 (ν x )um + φ1 (ν x )hwm + h b1 (ν x )(ω x um + ∂ xx um ), n,k n,k 2 n,k n,k ˜ wm+1 = − ν x ω x φ1 (ν x )un,k m + φ0 (ν x )wm + hb1 (ν x )(ω x um + ∂ xx um ) n,k + hb˜ 2 (ν x )(ω2x un,k m+1 + ∂ xx um+1 ). From this we derive that ∂ xx un,k m = and then ∂ xx un,k m+1 =
n,k n,k un,k m+1 − φ0 (ν x )um − φ1 (ν x )hwm − ω2x un,k m , h2 bˆ 1 (ν x )
n,k n,k un,k m+2 − φ0 (ν x )um+1 − φ1 (ν x )hwm+1 − ω2x un,k m+1 . h2 bˆ 1 (ν x )
Then it follows that n,k n,k 2 n,k n,k ˜ wn,k m+1 = − ν x ω x φ1 (ν x )um + φ0 (ν x )wm + hb1 (ν x )(ω x um + ∂ xx um ) + hb˜ 2 (ν x )(ω2x un,k + ∂ xx un,k ) m+1
m+1
b˜ 1 (ν x ) n,k b˜ 1 (ν x ) b˜ 1 (ν x ) um+1 − φ0 (ν x )un,k φ1 (ν x )wn,k m − ˆ m ˆ ˆ hb1 (ν x ) hb1 (ν x ) b1 (ν x ) b˜ 2 (ν x ) n,k b˜ 2 (ν x ) b˜ 2 (ν x ) + um+2 − φ0 (ν x )un,k − φ1 (ν x )wn,k m+1 m+1 . hbˆ 1 (ν x ) hbˆ 1 (ν x ) bˆ 1 (ν x )
n,k = − ν x ω x φ1 (ν x )un,k m + φ0 (ν x )wm +
Assuming
b˜ 2 (ν x )un,k m+2
bˆ 1 φ0 (ν x ) b˜ 1 (ν x ) = , φ1 (ν x ) ¡ ¢ ¡ ¢ n,k ˜ ˜ − bˆ 1 (ν x )ν2x φ1 (ν x ) + b˜ 1 (ν x )φ0 (ν x ) un,k m + b1 (ν x ) − b2 (ν x )φ0 (ν x ) u
(28)
m+1 . hbˆ 1 (ν x ) + hb˜ 2 (ν x )φ1 (ν x ) Therefore, from (28), together with the symplectic conditions (13), (14) and order conditions of ERKN methods proposed in [18], we can construct several 2-stage SERKN methods of order 2. In this paper, we consider the following 2-stage SERKN methods of order 2 for the two directional discretizations, which can be expressed in Butcher tableaus:
then we have wn,k m+1 =
0 1
0 φ2 (νt )
0 0
¯ t) b(ν
φ2 (νt ) φ0 (νt )φ2 (νt ) φ1 (νt )
0 φ0 (νt )φ2 (νt ) φ1 (νt ) − φ1 (νt )
0 1
0 φ2 (ν x )
0 0
ˆ x) b(ν ˜ x) b(ν
φ2 (ν x ) φ0 (ν x )φ2 (ν x ) φ1 (ν x )
0 φ0 (ν x )φ2 (ν x ) φ1 (ν x ) − φ1 (ν x )
b(νt )
(29)
(30)
11
W. Shi, X. Wu, and J. Xia
1 − φ0 (ν) where φ2 (ν) = . ν2 Note that as ωt → 0, ω x → 0, the above two SERKN methods (29) and (30) for time and space discretizations, respectively, reduce to the same classical 2-stage symplectic RKN (SRKN) method of order 2: 0
0
0
1
1 2
0
1 2
0
1 2
1 2
This coincides with the leap-frog method qn+1 − 2qn + qn−1 = f (qn ) h2 for problem p˙ = f (q), q˙ = p. The leap-frog method can be translated from the following 1-stage SRKN method of order 2 by exchanging the roles of p and q: 1 2
0 (31)
1 2
1 With the discretizations of the SERKN methods (29) and (30), (19) becomes the following scheme Umn,1 = unm , n,1 n,1 Um+1 − 2Umn,1 + Um−1 n,1 , (∂ U) = xx m 2h2 φ2 (ν x ) ¡ 2 n,1 ¢ 0 n,2 n n 2 n,1 n,1 U = φ (ν )u + φ (ν )τv + τ φ (ν ) ω U + ∂ U − V (U ) , 0 t 1 t 2 t xx m m m t m m m n,2 U n,2 − 2U n,2 + Um−1 n,2 (32) (∂ xx U)m = m+1 2 m , 2h φ2 (ν x ) ¡ 2 n,1 ¢ n+1 0 n n 2 n,1 n,1 um = φ0 (νt )um + φ1 (νt )τvm + τ φ2 (νt ) ωt Um + ∂ xx Um − V (Um ) , n+1 vm = − νt¡ωt φ1 (νt¡)unm + φ0 (νt )vnm ¢ ¡ ¢¢ 0 0 + τ b1 (νt ) ω2t Umn,1 + ∂ xx Umn,1 − V (Umn,1 ) + b2 (νt ) ω2t Umn,2 + ∂ xx Umn,2 − V (Umn,2 ) . In fact, the scheme (32) gives the discretizations for (∂ xx u)nm and (∂tt u)nm by (∂ xx u)nm = and (∂tt u)nm =
unm+1 − 2unm + unm−1 2h2 φ2 (ν x )
n n−1 un+1 m − 2um + um , respectively. This also deduces an extended leap-frog scheme: 2 2τ φ2 (νt ) n n−1 un+1 un − 2un + unm−1 0 m − 2um + um = m+1 2 m − V (unm ). 2 2τ φ2 (νt ) 2h φ2 (ν x )
(33)
We refer to the scheme (33) as Eleap-frogI and use this scheme to solve a linear wave equation (Problem 1) in the next section. Meanwhile, the discretizations of w and v are given by wnm+1 =
n un+2 unm+2 − unm m − um , vn+1 . m = 2hφ1 (ν x ) 2τφ1 (νt )
Furthermore, we have the following theorem.
W. Shi, X. Wu, and J. Xia
12
Theorem 4.1. By simply exchanging the roles of u and w in space and u and v in time, the MSERKN scheme given by the following 1-stage SERKN methods of order 2 for the two directional discretizations can be transformed into the MSERKN scheme (33): 1 2
¯ t) b(ν
φ21 ( ν2t )
0 ± 2φ0 ( ν2t )
b(νt )
φ1 ( ν2t )
1 2
0 ± φ21 ( ν2x ) 2φ0 ( ν2x )
ˆ x) b(ν ˜ x) b(ν
φ1 ( ν2x )
(34)
(35)
Proof. Firstly, the symplecticity of the schemes (34) and (35) can be directly obtained from Theorem 3.1. Then, we only show the details for the space discretization, and those for the time discretization are similar. Applying the 1-stage SERKN method (35) with the spatial discretization yields 1 νx νx U1 = φ0 ( 2 )um + 2 φ1 ( 2 )hwm , (36) ˆ x ) f1 , um+1 = φ0 (ν x )um + φ1 (ν x )hwm + h2 b(ν w ˜ m+1 = −ν x ω x φ1 (ν x )um + φ0 (ν x )wm + hb(ν x ) f1 , where f1 = f (U1 ) = ω2x U1 + ∂ xx U1 . From the last two equations in (36), we get ˆ x ) f1 um+1 =φ0 (ν x )um + φ1 (ν x )hwm + h2 b(ν ˆ x ) wm+1 + ν x ω x φ1 (ν x )um − φ0 (ν x )wm =φ0 (ν x )um + φ1 (ν x )hwm + h2 b(ν ˜ x) hb(ν ˆ x) ˆ x) ˆ x) b(ν b(ν b(ν =φ0 (ν x )um + φ1 (ν x )hwm + h wm+1 + ν2x φ1 (ν x )um − h φ (ν )w ˜ x) ˜ x) ˜ x) 0 x m b(ν b(ν b(ν φ1 (ν x ) =um + h (wm + wm+1 ). 1 + φ0 (ν x ) This shows um+1 − um = h
φ1 (ν x ) (wm + wm+1 ). 1 + φ0 (ν x )
(37)
The first and last equations of (36) imply ˜ x ) f (U1 ) wm+1 = − ν x ω x φ1 (ν x )um + φ0 (ν x )wm + hb(ν ˜ x )(ω2x U1 + ∂ xx U1 ) = − ν x ω x φ1 (ν x )um + φ0 (ν x )wm + hb(ν µ ¶ νx νx (38) 2 21 ˜ = − ν x ω x φ1 (ν x )um + φ0 (ν x )wm + hb(ν x ) ω x φ0 ( )um + ω x φ1 ( )hwm + ∂ xx U1 2 2 2 νx =wm + hφ1 ( )∂ xx U1 . 2 Following the approach to obtaining the leap-frog method by exchanging the roles of p and q in the SRKN method (31), we also exchange the roles of u and w. Then it can be deduced from the first equation of (36), (37) and (38) that 1 νx νx W1 = φ0 ( )wm + φ1 ( )h∂ xx um , 2 2 2 φ1 (ν x ) (39) wm+1 = wm + h (∂ xx um + ∂ xx um+1 ), 1 + φ0 (ν x ) νx u )W1 . m+1 = um + hφ1 ( 2
13
W. Shi, X. Wu, and J. Xia
The first equation of (39) implies wm+1/2
µ ¶ νx 1 νx νx φ1 (ν x ) = φ0 ( )wm + φ1 ( )h∂ xx um = φ0 ( ) wm + h ∂ xx um , 2 2 2 2 1 + φ0 (ν x )
(40)
1 φ1 (ν x ) ∂ xx um . ν x wm+1/2 − h φ0 ( 2 ) 1 + φ0 (ν x )
(41)
and then wm =
Inserting (41) into the second equation of (39) gives wm+1 = and furthermore, wm = That is, wm−1/2
1 φ1 (ν x ) ∂ xx um+1 , wm+1/2 + h φ0 ( ν2x ) 1 + φ0 (ν x ) φ1 (ν x ) 1 wm−1/2 + h ∂ xx um . φ0 ( ν2x ) 1 + φ0 (ν x )
µ ¶ φ1 (ν x ) νx 1 νx νx ∂ xx um = φ0 ( )wm − φ1 ( )h∂ xx um . = φ0 ( ) wm − h 2 1 + φ0 (ν x ) 2 2 2
Therefore, it follows from (40) and (42) that νx wm+1/2 + wm−1/2 = 2φ0 ( )wm , 2 ν wm+1/2 − wm−1/2 = φ1 ( x )h∂ xx um . 2 On the other hand, the last equation of (39) means wm+1/2 =
1 2
νx u − um − (um − um−1 ) wm+1/2 − wm−1/2 = m+1 = φ1 ( )h∂ xx um . hφ1 ( ν2x ) 2
(44)
um+1 − um−1 2hφ1 (ν x ) u m+1 − 2um + um−1 . (∂ xx u)m = 2h2 φ2 (ν x ) wm =
For the time discretization, it also can be proved that (∂tt u)n = holds.
(43)
um+1 − um , which leads to hφ1 ( ν2x )
νx um+1 − um + (um − um−1 ) = 2φ0 ( )wm , wm+1/2 + wm−1/2 = hφ ( νx ) 2
It follows from (44) that
(42)
n un+1 − 2un + un−1 un+2 m − um n+1 and v = m 2τ2 φ2 (νt ) 2τφ1 (νt ) ¤
4.2. Extended leap-frog scheme II: An explicit multi-symplectic ERKN-PRK scheme In this subsection, we construct another explicit multi-symplectic scheme based on Theorem 3.3. We apply the SPRK method (22) to the space discretization, which is equivalent to the St¨ormer-Verlet that gives the second order central difference of ∂ xx Umn,k , and apply the explicit 2-stage SERKN method (29) to the time
14
W. Shi, X. Wu, and J. Xia
discretization. Then the multi-symplectic scheme (27) becomes n,1 U = unm , m n,1 n,1 Um+1 − 2Umn,1 + Um−1 n,1 (∂ U) = , xx m h2 ¡ ¢ 0 Umn,2 = φ0 (νt )unm + φ1 (νt )τvnm + τ2 φ2 (νt ) ω2t Umn,1 + ∂ xx Umn,1 − V (Umn,1 ) , n,2 U n,2 − 2Umn,2 + Um−1 n,2 (∂ xx U)m = m+1 , h2 ¡ 2 n,1 ¢ 0 n+1 n n 2 n,1 n,1 um = φ0 (νt )um + φ1 (νt )τvm + τ φ2 (νt ) ωt Um + ∂ xx Um − V (Um ) , vn+1 = − νt¡ωt φ1 (νt¡)unm + φ0 (νt )vnm ¢ ¡ ¢¢ m 0 0 + τ b1 (νt ) ω2t Umn,1 + ∂ xx Umn,1 − V (Umn,1 ) + b2 (νt ) ω2t Umn,2 + ∂ xx Umn,2 − V (Umn,2 ) .
(45)
Just like the situation of the scheme Eleap-frogI, it can be observed that this scheme is equivalent to the following form n n−1 un+1 unm+1 − 2unm + unm−1 0 m − 2um + um = − V (unm ), (46) 2 2 2τ φ2 (νt ) h and the discretizations of w and v are wnm+1 =
n unm+2 − unm un+2 m − um , vn+1 . m = 2h 2τφ1 (νt )
It is also equivalent to the concatenation of the 1-stage SERKN method (34) of order 2 and SPRK method (22) applied to the temporal and spatial discretizations, respectively. We refer to the above scheme (46) as Eleap-frogII, and in the next section, we use this scheme to solve a nonlinear wave equation (Problem 2). 4.3. Analysis of linear stability In this subsection, we pay attention to another significant aspect for those schemes proposed in this paper, namely, the stability. It is known that the attractiveness of explicit schemes relies on the computational efficiency, but such efficiency is always at the cost of stability. Although an application of the linear stability analysis to nonlinear equations cannot be justified, it is found to be effective in numerical experiments. Using the well-known Fourier method [40], we give some linear stability analysis for the explicit multi-symplectic Eleap-frogI scheme, and the stability for scheme Eleap-frogII can be analyzed in the similar way. Set unm = uˆ n exp{iξmh}, vnm = vˆ n exp{iξmh}, Umn,1 = Uˆ n,1 exp{iξmh}, Umn,2 = Uˆ n,2 exp{iξmh}, etc., √ 0 where ξ ∈ R and i = −1. Then, substituting the above expressions into scheme (32) with V (Umn,i ) = Umn,i , we get n,1 n Uˆ = uˆ , ¶ µ iξh −iξh Uˆ n,2 = φ (ν )ˆun + φ (ν )τˆvn + τ2 φ (ν ) ω2 + (e − 2 + e ) − 1 Uˆ n,1 , 0 t 1 t 2 t t 2h2 φ2 (ν x ) n+1 n,2 uˆ = Uˆ , n+1 = −νt ωt φ1 (νt )ˆun + φ0 (νt )ˆvn vˆ µ ¶ ³ ´ ³ ´ (eiξh − 2 + e−iξh ) (eiξh − 2 + e−iξh ) 2 n,1 2 n,2 ˆ ˆ − 1 U + b2 (νt ) ωt + −1 U . +τ b1 (νt ) ωt + 2h2 φ2 (ν x ) 2h2 φ2 (ν x ) Eliminating Uˆ n,1 , Uˆ n,2 , we obtain
µ
uˆ n+1 vˆ n+1
¶
µ =G
uˆ n vˆ n
¶ ,
15
W. Shi, X. Wu, and J. Xia
where
τ2 φ2 (νt )θ 2 1 − τ φ2 (νt ) + h2 φ2 (ν x ) G= G21
φ1 (νt )τ
, 2 τ φ (ν )θ 2 t 1 − τ2 φ2 (νt ) + 2 h φ2 (ν x )
with θ = (cos ξh − 1), and G21 =
τφ1 (νt )θ + τ3 b2 (νt )νt2 φ2 (νt )θ(ω2t − 2) τ3 b2 φ2 (νt )θ2 − φ1 (νt )τ + b2 (νt )φ2 (νt )τ3 (1 − ω2t ). + h2 φ2 (ν x ) h4 φ22 (ν x )
We can easily verify that G11G22 − G12G21 = 1. Hence, the eigenvalues of G = G(ξ, θ, ν x , νt ) satisfy λ2 − (G11 + G22 )λ + (G11G22 − G12G21 ) = 0. That is, λ2 − (G11 + G22 )λ + 1 = 0.
(47)
Obviously, the roots λi (i = 1, 2) of equation (47) satisfy |λi | ≤ 1 if and only if |G11 + G22 | ≤ 2. Then the method satisfies the von Neumann condition, which is necessary for the stability [40]. Assuming |G11 + G22 | ≤ 2 that ¯ ¯ 2 ¯ ¯ ¯1 − τ2 φ2 (νt ) + τ φ2 (νt )θ ¯ ≤ 1, ¯ 2 h φ2 (ν x ) ¯ which is equivalent to
τ2 φ2 (νt )θ ≤ 1. (48) h2 φ2 (ν x ) τ π Denote r = , and suppose 0 < r ≤ α and 0 < νt ≤ ν x ≤ . Since φ2 (ν) is a monotonically decreasing h 2 π 1 function on the interval [0, ] with φ2 (0) = and −2 ≤ θ ≤ 0, we obtain 2 2 −1 ≤ 1 − τ2 φ2 (νt ) +
1 1/2 τ2 φ2 (νt )θ − 2α2 < 1 − τ2 φ2 (νt ) + 2 < 1. 2 φ2 (π/2) h φ2 (ν x ) √ 6 1 1/2 Obviously, if 0 < α ≤ , then − 2α2 ≥ −1, and the inequality (48) is satisfied. π 2 φ2 (π/2) Furthermore, it also gives a sufficient condition for the stability of Eleap-frogI. In fact, conditions 0 < √ 6 π τ2 φ2 (νt )θ α≤ and νt ≤ ν x ≤ make −1 < 1 − τ2 φ2 (νt ) + 2 < 1, which implies that the equation (47) has π 2 h φ2 (ν x ) two distinct conjugate roots λi with |λi | = 1, i = 1, 2. Therefore, the method is linearly stable if √ √ 6 6 π 0 0 is just the initial condition translated through space: u(x, t) = u(x − c0 t, 0). ∂D ± ∂D 0 . An important quantity of dispersive systems is the group velocity defined by, vg (k) = ω (k) = − ∂k ∂ω The group velocity is usually a function of the wave number k and characterizes the speed of energy transport 00 in wave packets. In this case, ω (k) , 0, and waves with different wave numbers travel with different phase velocities v p , which brings the group velocity dispersion (GVD) and ultimately causes the spatial spreading of the wave packet. The dispersive properties of the discretizations can be studied by calculating their numerical dispersion relations. To do so, we start with the discretizations of the linearized sine-Gordon equation (55) associated with the Eleap-frogI and Eleap-frogII schemes in their single equation format D2 un−1 D2 un Eleap − frogI : t m − x m−1 + χunm = 0, 2φ2 (νt ) 2φ2 (ν x ) (57) D2 un−1 Eleap − frogII : t m − D2x unm−1 + χunm = 0, 2φ2 (νt )
22
W. Shi, X. Wu, and J. Xia
n un+1 un − unm m − um and D x unm = m+1 are the finite difference operators. τ h P i(Km−Ωn) Assume a discrete general solution to a linear PDE of the form Ae , K = kh, Ω = ωτ, where
where Dt unm =
K
each mode is a solution provided K and Ω satisfy a numerical dispersion relation DN (Ω, K) = 0. For the equation (57), the corresponding numerical dispersion relations are given by µ√ ¶2 µ √ ¶2 2φ2 (νt ) 2φ2 (ν x ) Ω K sin − sin − χ = 0, Eleap − frogI : φ (ν )τ 2 φ (ν )h 2 µ √2 t ¶2 µ 2 x ¶2 2φ2 (νt ) Ω 2 K sin − sin − χ = 0. Eleap − frogII : φ2 (νt )τ 2 h 2
(58)
Both these two schemes preserve the form of the analytical dispersion relation (56). The exact relationship is given by the following theorem. Theorem 5.1. The multi-symplectic schemes (58) qualitatively preserves the dispersion relation of the linear PDE. Specifically, there exist diffeomorphisms ψ1 and ψ2 satisfying the exact dispersion relationship ¡ ¢ DN (Ω, K) = D ψ1 (Ω), ψ2 (K) = ψ21 − ψ22 − χ = 0, (59) where
¶ µ√ √ ¡ ¢ 2φ2 (νt ) Ω 2φ2 (ν x ) K sin , sin , Eleap − frogI : ψ1 (Ω), ψ2 (K) = φ2 (νt )τ 2 φ2 (ν x )h 2 ¶ µ√ ¡ ¢ 2φ2 (νt ) Ω 2 K sin , sin , Eleap − frogII : ψ1 (Ω), ψ2 (K) = φ2 (νt )τ 2 h 2
(60)
for −π < K, Ω < π. Solving the numerical dispersion relations (58) gives s φ2 (νt )r2 K φ2 (νt )hr Eleap − frogI : Ω(K) = ±2 arcsin sin2 + χ, φ2 (ν x ) 2 2 r K φ (ν )hr Eleap − frogII : Ω(K) = ±2 arcsin 2φ2 (νt )r2 sin2 + 2 t χ. 2 2 The significance of Theorem 5.1 for the non-dispersive equation (55) (with χ = 0) can be seen in Figures 0 00 9–11, which show the dispersion curves Ω(K), the group velocity Ω (K), the group velocity dispersion Ω (K) for the exact (56), and numerical (58). In the experiment, two different values of the ratio r = τ/h, with h = 0.01 are used. The exact relation is given by Ω(K) = ±rK. The significance of Theorem 5.1 for the linearized sine-Gordon equation (55) (with χ = 1) can be seen in √ Figures 12–14. The exact relation is given by Ω(K) = ±r K 2 + h2 . From Figures 9–14, we can see that the existence of diffeomorphisms (60) is not enough to preserve the qualitative features of the analytical solution. All of the schemes introduce numerical dispersions. From Figures 10 and 13, because Ω(K) and K both have direction, we can see that the dispersion relation curves are all monotonically increasing at rates given by their numerical group velocities. Similarly, Figures 11 and 14 display that the group velocity curves (see Figures 10 and 13 ) are monotonically decreasing. Therefore, considering the numerical solutions of the linear wave equations given by Eleap-frogI and Eleap-frogII schemes, we further obtain the following corollaries. ∂DN ± ∂DN 0 at (K, Ω) of the two schemes Corollary 5.1. The numerical group velocities Vg (k) = Ω (K) = − ∂K ∂Ω ∂D ± ∂D 0 Eleap-frogI and Eleap-frogII have the same sign as the group velocity vg (k) = ω (k) = − for the ∂k ∂ω associated pair (k, ω). Corollary 5.2. Both Eleap-frogI and Eleap-frogII schemes introduce negative dispersion, and when the wave number is increasing, the numerical group velocities are monotonically decreasing.
23
W. Shi, X. Wu, and J. Xia
r=0.7
r=0.4 1.5 Eleap−frogI Ω⋅K>0 Eleap−frogI Ω⋅K0 exact Ω⋅K0 Eleap−frogI Ω⋅K0 exact Ω⋅K0 Eleap−frogI Ω⋅K0 exact Ω⋅K0 Eleap−frogI Ω⋅K0 exact Ω⋅K0
’
0.3
Eleap−frogI Ω (K)⋅K>0
0.2
Eleap−frogI Ω (K)⋅K0 Eleap−frogII Ω⋅K0 exact Ω⋅K0 Eleap−frogII Ω⋅K0 exact Ω⋅K0 Eleap−frogII Ω⋅K0 exact Ω⋅K0 Eleap−frogII Ω⋅K0 exact Ω⋅K0
0.4
’
Eleap−frogII Ω (K)⋅K>0
Eleap−frogII Ω’(K)⋅K0
0.2
exact Ω (K)⋅K0
0.4
’
Ω (K)
’
Eleap−frogII Ω (K)⋅K
