- Email: [email protected]
Multi-symplectic integration methods for Hamiltonian PDEs
Future Generation Computer Systems 19 (2003) 395–402
Multi-symplectic integration methods for Hamiltonian PDEs Brian E. Moore a,∗,1 , Sebastian Reich b a
Department of Mathematics and Statistics, University of Surrey, Guildford, Surrey GU2 7XH, UK b Department of Mathematics, Imperial College, London SW2 7BZ, UK
Abstract Recent results on numerical integration methods that exactly preserve the symplectic structure in both space and time for Hamiltonian PDEs are discussed. The Preissman box scheme is considered as an example, and it is shown that the method exactly preserves a multi-symplectic conservation law and any conservation law related to linear symmetries of the PDE. Local energy and momentum are not, in general, conserved exactly, but semi-discrete versions of these conservation laws are. Then, using Taylor series expansions, one obtains a modified multi-symplectic PDE and modified conservation laws that are preserved to higher order. These results are applied to the nonlinear Schrödinger (NLS) equation and the sine-Gordon equation in relation to the numerical approximation of solitary wave solutions. © 2002 Elsevier Science B.V. All rights reserved. Keywords: Multi-symplectic PDEs; Preissman box scheme; Conservation laws; Modified equations; Backward error analysis
1. Introduction When approximating the solution of a differential equation using a numerical integrator, it is often desirable to preserve as many properties of the exact evolution equation as possible. In particular, it is advantageous to preserve the symplectic structure of Hamiltonian systems 0 d Id Jyt = ∇y H(y) with J = , (1) −Id 0d where y(t) ∈ R2d and H : R2d → R is the Hamiltonian (energy) of the system. Numerical schemes that have this property, known as symplectic integrators, have proven to be both accurate and efficient in the ∗ Corresponding author. E-mail addresses: [email protected] (B.E. Moore), [email protected] (S. Reich). 1 This work was accomplished while the author was visiting Imperial College, London.
long-time approximation of solutions to Hamiltonian ODEs [3,10,13,18,20]. In addition, symplectic integrators lead to a useful interpretation of backward error analysis, in which a modified equation, the equation solved by the numerical scheme to higher order than the original equation, is used to understand the discretization error induced by the numerical scheme. Recently, the idea of symplectic integration has been extended in a novel way to Hamiltonian PDEs [8,14], by exploiting the multi-symplectic structure of the PDE Kzt + Lzx = ∇z S(z),
(2)
where K and L are skew-symmetric, z(x, t) is the vector of state variables, and S(z) is a smooth function. Specifically, Bridges and Reich [8] apply a symplectic integrator to each independent variable to get a multi-symplectic integrator which preserves the multi-symplectic structure of (2). This approach has paved the way to a useful backward error analysis for multi-symplectic integrators, such that a modified
0167-739X/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 7 - 7 3 9 X ( 0 2 ) 0 0 1 6 6 - 8
396
B.E. Moore, S. Reich / Future Generation Computer Systems 19 (2003) 395–402
multi-symplectic PDE is satisfied by the numerical solution to higher order in both space and time [15,16]. The modified PDE can also be used to derive modified conservation laws of energy and momentum, leading to a better understanding of the numerical scheme. The focus of the present paper is on the multisymplectic Preissman box scheme, which is the space–time version of the implicit midpoint scheme for ODEs. This scheme has been considered in the context of multi-symplectic PDEs in [8,19], and it is commonly used in the area of hydraulics [1,2]. It is shown that this scheme preserves a conservation law, which is directly related to linear symmetries associated with (2) through Noether theory, as well as a multi-symplectic conservation law, exactly. Backward error analysis is also used to show that the modified PDE satisfies these conservation laws and the conservation laws of energy and momentum to higher order. The outline is as follows. We first introduce the conservation laws satisfied by (2) in Section 2, and discuss the box scheme as well as its conservation properties in Section 3. Then, Section 4 considers the modified equations and conservation laws obtained through backward error analysis. The results are then applied in Section 5 to a model problem known as the nonlinear Schrödinger (NLS) equation and soliton solutions of the sine-Gordon equation.
2. Conservation laws There are several properties of (2) that take the form of conservation laws and should be considered when using a numerical integrator to approximate solutions. First, according to Bridges [4,5], a multi-symplectic conservation law ∂t ω + ∂x κ = 0,
where ω = dz ∧ K dz and
κ = dz ∧ L dz is derived directly from (2). Here, we use wedge product notation such that dz denotes the vector of differentials, and let ∂t and ∂x denote differentiation with respect to time and space. One can restrict this conservation law to the space of solutions by making the ansatz dz = zt dt + zx dx.
This yields the conservation law ∂t (zTt Kzx ) + ∂x (zTt Lzx ) = 0.
(3)
Following the analysis in [4,5], we can also derive local conservation laws of energy and momentum. Using the time invariance of (2), an energy conservation law can easily be derived by taking the product of zTt with (2). Then we obtain the energy conservation law ∂t E + ∂x F = 0,
where E = S(z) + zTx Lz/2 and
F = −zTt Lz/2 are known, respectively, as the energy density and the energy flux. Similarly, the spatial invariance of (2) can be used to take the product of zTx with (2), which gives the momentum conservation law ∂x G + ∂t I = 0,
where G = S(z) + zTt Kz/2 and
I = −zTx Kz/2. It is interesting to note that the conservation law (3) is equivalent to the one obtained from differentiation of the energy and momentum conservation laws, i.e. ∂x (Et + Fx ) − ∂t (Gx + It ) = 0.
(4)
According to Bridges [4], there may also be additional conservation laws for (2) related to linear symmetries. Essentially, this follows from Noether theory and derivation of the multi-symplectic PDE (2) from a Lagrangian functional 1 L = Ldt dx for L = zT (Kzt + Lzx ) − S(z). 2 (5) More specifically, take a linear one-parameter family of linear coordinate transformations [17] given by the group action Gε (z) = eεA z,
(6)
which is chosen in such a way that it is symplectic with respect to both ω and κ. Symplecticity is equivalent to AT K + KA = 0
and
AT L + LA = 0.
(7)
Using these two identities, the invariance of a Lagrangian under such a change of variables leads to ∂L 0= = S (z)Az dt dx ∂ε ε=0
B.E. Moore, S. Reich / Future Generation Computer Systems 19 (2003) 395–402
and we obtain the invariance condition
S (z)Az = 0.
397
where (8)
Direct application of (8) to the multi-symplectic PDE formulation (2) yields (Az)T Kzt + (Az)T Lzx = (Az)T ∇z S(z) = 0
zn+1/2,i+1/2 = 41 (zn+1,i+1 + zn+1,i + zn,i+1 + zn,i ). We will call this scheme a multi-symplectic integrator because it preserves exactly a multi-symplectic conservation law [8]
and using (7), this can be written as a conservation law
n+1/2,i n,i+1/2 δ+ + δ+ =0 t ω xκ
∂t (zT KAz) + ∂x (zT LAz) = 0.
for
(9)
Notice here that each of these conservation laws is preserved locally, and these formulations are made possible by the multi-symplectic structure of the PDE. This is in contrast to a Hamiltonian formulation which yields global conservation of densities and neglects the fluxes. We make the distinction here because the preservation of a local conservation law numerically is a stronger result than the global conservation that has been achieved in the past, and these are the issues discussed as we consider a numerical method.
ωn+1/2,i = dzn+1/2,i ∧ K dzn+1/2,i and κn,i+1/2 = dzn,i+1/2 ∧ K dzn,i+1/2 . Eq. (10) can also be derived using a discrete variational principle to get ∂ L = 0, ∂zn,i
L=
n
Ln+1/2,i+1/2
i
for the approximation 3. A multi-symplectic box scheme
n+1/2,i Ln+1/2,i+1/2 = 21 (zn+1/2,i+1/2 )T Kδ+ t z
We now present a centered scheme and discuss the preservation of these conservation laws. First, introduce the notation zn,i , which denotes the numerical approximation of z(xn , ti ) for n = 0, 1, 2, . . . , N and i = 0, 1, 2, . . . , T , where N is the number of grid points and T is the number of time steps. We also set x = xn+1 −xn and t = ti+1 −ti . In addition, define the midpoints zn+1/2,i = 21 (zn+1,i + zn,i ) n,i+1/2
z
=
1 n,i+1 2 (z
and
n,i
+z )
and the forward difference approximations zn+1,i − zn,i and x zn,i+1 − zn,i n,i zt ≈ δ+ . z = t t Using the implicit midpoint scheme [20] to discretize (2) in space and time separately yields the multi-symplectic Preissman box scheme [8]
n,i zx ≈ δ+ = xz
n+1/2,i n,i+1/2 Kδ+ + Lδ+ = ∇z S(zn+1/2,i+1/2 ), t z xz
(10)
n,i+1/2 + 21 (zn+1/2,i+1/2 )T Lδ+ xz
−S(zn+1/2,i+1/2 ) of (5). This is a generalization of the discrete multi-symplectic formulation of Marsden et al. [14] for first-order field theories. Provided the PDE (2) is invariant under an appropriate linear symmetry, the Preissman box scheme also satisfies a discrete version of the conservation law (9) exactly. To make this clear, multiply (10) by (Azn+1/2,i+1/2 )T to get n+1/2,i (zn+1/2,i+1/2 )T KAδ+ t z
n,i+1/2 +(zn+1/2,i+1/2 )T LAδ+ = 0, xz
where we make use of the symplectic and invariance conditions (7) and (8). Then it is easy to show that this is equivalent to the discrete conservation law n+1/2,i T δ+ ) KAzn+1/2,i ) t ((z
n,i+1/2 T +δ+ ) LAzn,i+1/2 ) = 0. x ((z
For linear PDEs, where S(z) = (zT Bz)/2 for B symmetric, this scheme also preserves fully discrete
398
B.E. Moore, S. Reich / Future Generation Computer Systems 19 (2003) 395–402
energy and momentum conservation laws exactly [8]. However, for nonlinear systems in general, exact conservation of discrete energy and momentum conservation laws cannot be guaranteed, but certain semi-discretized energy and momentum conservation laws are still preserved exactly. For example, similar to the analysis of [16], we can take the product of n+1/2 T (zt ) with the spatially discrete equation n+1/2
Kzt
n n+1/2 + Lδ+ ) x z = ∇S(z
to get n+1/2 T
n n+1/2 ) Lδ+ ), x z = ∂t S(z
(zt
n+1/2 T n+1/2 ) Kzt
because K skew-symmetric implies (zt = 0. Now, since n T n −δ+ x ((z ) Lzt ) n+1/2 T
n n+1/2 T n ) Lδ+ ) Lδ+ x z − (z x zt ,
= (zt
concept of modified equations through backward error analysis.
4. The modified equations Based on the ideas of [16], we can perform a useful backward error analysis for the discretization (10), such that the modified equations are also multi-symplectic and satisfy modified versions of the conservation laws in Section 2. Using Taylor series expansions about ti+1/2 we find that δ+ t z(ti ) = zt (ti+1/2 ) +
and an analogous expansion holds for δ+ x z(xn ), which implies that the scheme is second-order in both space and time. Truncating these expansions and substituting into (10) yields the modified system of equations Kzt +
we have
2 4 1 24 t zttt (ti+1/2 ) + O(t )
2 1 24 t Kzttt
+ Lzx +
2 1 24 x Lzxxx
= ∇z S(z), (11)
n+1/2 T
n n+1/2 T n 1 ) Lδ+ ) Lδ+ x z = 2 ∂t ((z xz )
(zt
n T n − 21 δ+ x ((z ) Lzt ).
This leads to the spatially discrete energy conservation law n ∂t En+1/2 + δ+ xF =0
˜ zt + L˜ ˜ zx = ∇z˜ S(˜ ˜ z) K˜
with E
n+1/2 n
F =
= S(z
n+1/2
)−
which can be viewed as a generalized (higher order) multi-symplectic PDE and is satisfied by the numerical solution up to O(t 4 + x4 ). Eq. (11) can also be written in the form of a standard (first-order) multi-symplectic PDE
n 1 n+1/2 T ) Lδ+ xz 2 (z
for z˜ = [z, p, q, r, s]T and and
n 1 n T 2 (z ) Lzt .
A semi-discrete momentum conservation law is derived in the same way by reversing the roles of space and time, i.e. discretize the PDE in time only and take i+1/2 T the product of the resulting equation with (zx ) . Thus, a semi-discrete momentum conservation law i ∂x Gi+1/2 + δ+ t I =0
with
˜ K=
is satisfied for i Gi+1/2 = S(zi+1/2 ) − 21 (zi+1/2 )T Kδ+ t z
S˜ = S +
and
I i = 21 (zi )T Kzix . To understand the energy and momentum conservation of fully discretized equations we must resort to the
L˜ =
2 T 2 T 1 1 24 t q Kp + 24 x s Lr
K
0
1 2 24 t K
0
0
1 − 24 t 2 K
0
0
1 2 24 t K
0
0
0
0
0
0
0
0
L
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1 − 24 x2 L
1 2 24 x L
0
0
0
0
0 0 0, 0 0 0 0 1 2 24 x L 0 0 , 0 0
B.E. Moore, S. Reich / Future Generation Computer Systems 19 (2003) 395–402
399
where p = zt , q = pt , r = zx and s = rx . The modified conservation laws of energy and momentum are easily found using this modified multi-symplectic PDE or (11) directly. The modified form of (3) can then be computed from these modified energy and momentum conservation laws using (4). The modified PDE (11) can also be derived from a modified Lagrangian density
used to investigate computational solitary wave solutions for the sine-Gordon equation.
˜ 1 = 1 (zT Kzt − L 2
where ψ is complex-valued. Letting ψ = v + iw, we can rewrite this as the Hamiltonian system of equations
2 T 1 24 t zt Kztt
+ zT Lzx
1 − 24 x2 zTx Lzxx ) − S(z).
In general, the modified Lagrangian density for any number of modifications ρ, can be written as ˜ρ = L
ρ j=0
+
(−1)j t 2j (2)2j+1 (2j + 1)!
ρ j=0
j
j+1
(∂t z)T K∂t
z
(−1)j x2j (∂j z)T L∂xj+1 z − S(z), (2)2j+1 (2j + 1)! x
j
where ∂η denotes the jth derivative with respect to η. Given a linear symmetry, this modified Lagrangian density is also invariant under the transformation (6). For example, the associated conservation law for ρ = 1 is easily obtained by multiplying the modified equation (11) by (Az)T , which yields zT KAzt +
2 T T 1 24 t z KAzttt + z LAzx 1 + 24 x2 zT LAzxxx = 0.
Then the identity zT KAzttt = ∂t (zT KAztt − 21 zTt KAzt )
5.1. The nonlinear Schrödinger equation Consider the NLS equation iψt + ψxx + V (|ψ|2 )ψ = 0,
vt = −wxx − wV (v2 + w2 ), wt = vxx + vV (v2 + w2 ). Defining ψx = σ + iφ, this system of equations can also be written as the multi-symplectic PDE (2) where z = [v, w, σ, φ]T J 02 02 −I2 K= , L= 02 02 I2 0 2 and S(z) = 21 [σ 2 + φ2 + V(v2 + w2 )], where J is given in (1) with d = 1. This system is invariant under the action of the one-parameter group of rotations SO(2) [7], given by Gθ (z) = Rθ z, for cos θ − sin θ 0 0 cos θ 0 0 sin θ . Rθ = 0 cos θ − sin θ 0 0
sin θ
0
cos θ
and a similar identity for zT LAzxxx , imply the modified conservation law
Now, define the matrix A such that
0 = ∂t [zT KAzt +
Az =
2 T 2 T 1 1 24 t z KAztt − 48 t zt KAzt ] 1 1 +∂x [zT LAzx + 24 x2 zT LAzxx − 48 x2 zTx LAzx ].
5. Applications We can now apply these results to some model problems. We first consider the NLS equation as a Hamiltonian PDE with a rotational symmetry. The second example outlines how backward error analysis can be
d Gθ (z)|θ=0 = [−w, v, −φ, σ]T . dθ
Then zT KAz = w2 + v2
and
zT LAz = 2(vφ − wσ)
satisfy the conservation law (9). Now it is straight forward to derive the results of the previous sections for a Preissman box scheme discretization of this equation, including discrete conservation laws and modified equations with their conservation laws.
400
B.E. Moore, S. Reich / Future Generation Computer Systems 19 (2003) 395–402
5.2. Soliton solutions The NLS equation has a family of solitary wave solutions as do several other equations that can be written in the from (2), such as the KdV equation, Boussinesq equations, and nonlinear wave equations (cf. [4–7,21]). In addition, a thorough analysis of the error growth due to time discretization for these solutions has been performed for the KdV equation [9] and the NLS equation [12]. Specifically, Durán and Sanz-Serna [12] have shown that the implicit midpoint time discretization has better error propagation mechanisms than other non-conservative schemes. In particular, they have shown that the time-discrete solutions of the NLS equation, which have been initialized with an exact soliton profile, are made up of a modified solitary wave, that have a different wave speed and amplitude than the exact soliton, as well as a complementary term that grows linearly in time and higher order terms. Using the results of the previous sections, it may be possible to show that the numerical solutions of the fully discretized NLS equation initialized with an exact soliton profile are made up of a modified soliton, which has corrections in both space and time, along with the complementary error and higher order terms. From a slightly different perspective, one may be able to show that the modified equations have solitary wave solutions. Then, initializing the numerical scheme with the new modified soliton profile would yield numerical solutions that are preserved to higher order in both space and time. Results in that spirit have been achieved concerning soliton solutions of the sine-Gordon equation utt = uxx − sin u,
(12)
which can be written in the from (2) with z = [u, v, w]T
0 −1 0 0 K = 1 0 0, L = 0 0 0 0 1 1 1 S(z) = v2 − w2 − cos u. 2 2
0 0 0
−1 0 0
and
It is well known that there is a family of solitary wave
solutions for (12) given by
±(x − ct − ξ0 ) , u(x, t) = 4 tan −1 ±exp (1 − c2 )1/2
(13)
where |c| < 1 is the wave speed, and ξ0 is given by the initial conditions [11]. For our purposes we take ξ0 = 0. Using the explicit first-order in space and time multi-symplectic Euler box scheme [16] to discretize gives n,i n,i δ+ + δ+ = − sin un,i , t v xw n,i δ− = −wn,i , xu
n,i δ− = vn,i , t u
(14)
where we have used the backward differences n,i = δ− xu
un,i − un−1,i un,i − un,i−1 n,i and δ− . = t u x t
Notice that eliminating v and w gives the familiar second-order scheme δ2t un,i = δ2x un,i − sin un,i ,
(15)
− where δ2η = δ+ η δη for η = x, t. Hence, the numerical n,i solution u is accurate to second-order in both t and x, but using the scheme (14) yields first-order error because v and w are only approximated to first-order. To make this more clear, backward error analysis [16] can be used for (14) to get the modified equations
vt + 21 tvtt + wx + 21 xwxx = − sin u,
(16)
ut − 21 tutt = v,
(17)
ux − 21 xuxx = −w.
(18)
Since the numerical solution is already second-order accurate in u, only (17) and (18) need to be verified along the numerical solution to check for second-order accuracy with respect to the modified equations. Computing higher order corrections gets increasingly more complicated. However, the modified equation for (15), given by utt +
2 1 12 t utttt
= uxx + 21 x2 uxxxx − sin u
(19)
is easily obtained using Taylor expansions and is satisfied by the numerical solution to fourth-order. Now it is reasonable to ask if this modified system of equations has soliton solutions, which can be answered by
B.E. Moore, S. Reich / Future Generation Computer Systems 19 (2003) 395–402
401
Fig. 1. Error plotted against t for the original initial condition (dashed), and the modified initial condition for one modification (chain) and three modifications (solid).
setting u(x, t) = u(x − ct) = u(ξ) in (19) and finding solutions of the boundary value problem (c2 − 1)u −
2 1 12 (x
− c4 t 2 )u = − sin u.
(20)
Though it is not obvious how to solve this equation for a heteroclinic connection, we were able to find such a connection to high accuracy numerically. This approximation can be used to initialize the scheme (15) and this in turn yields a numerical solution for the PDE that is accurate to fourth-order in u with respect to both t and x. Using log–log scale in Fig. 1, we plot the solution error against t. For these results we have kept the ratio t/x fixed, so change in t directly corresponds to change in x. The dashed line shows that the numerical solution, obtained by initializing (14) with un,0 = u(xn , 0), vn,0 = ut (xn , 0) and
wn,0 = −ux (xn , 0) for u given in (13), is first-order accurate. The chain line shows that the same initial condition for (14) with v and w modified according to (17) and (18) yields second-order accuracy in both t and x. Finally, using a numerically computed solution of (20) as the initial condition for the scheme (15), the solid line shows approximately fourth-order convergence. The accuracy here is limited by the fact that the reference soliton solution had to be computed numerically, and this is a delicate procedure due to the heteroclinic nature of the solution.
Acknowledgements Sebastian Reich gratefully acknowledges partial financial support by European Commission funding for
402
B.E. Moore, S. Reich / Future Generation Computer Systems 19 (2003) 395–402
the Research Training Network “Mechanics and Symmetry in Europe”. References [1] M.B. Abbott, Computational Hydraulics, Pitman, London, 1979. [2] M.B. Abbott, D.R. Basco, Computational Fluid Dynamics, Longman Scientific, Harlow, 1989. [3] G. Benettin, A. Giorgilli, On the Hamiltonian interpolation of near to the identity symplectic mappings with application to symplectic integration algorithms, J. Stat. Phys. 74 (1994) 1117–1143. [4] T.J. Bridges, Multi-symplectic structures and wave propagation, Math. Proc. Camb. Phil. Soc. 121 (1997) 147–190. [5] T.J. Bridges, A geometric formulation of the conservation of wave action and its implications for signature and the classification of instabilities, Proc. Roy. Soc. London Ser. A 453 (1997) 1365–1395. [6] T.J. Bridges, G. Derks, Unstable eigenvalues and the linearization about solitary waves and fronts with symmetry, Proc. Roy. Soc. London Ser. A 455 (1999) 2427–2469. [7] T.J. Bridges, G. Derks, The symplectic Evans matrix, and the instability of solitary waves and fronts with symmetry, Arch. Rat. Mech. Anal. 156 (2001) 1–87. [8] T.J. Bridges, S. Reich, Multi-symplectic integrators: numerical schemes for Hamiltonian PDEs that conserve symplecticity, Phys. Lett. A 284 (2001) 184–193. [9] J. De Frutos, J.M. Sanz-Serna, Accuracy and conservation properties in numerical integration: the case of the Korteweg– de Vries equation, Numer. Math. 75 (1997) 421–445. [10] P. Deuflhard, J. Hermans, B. Leimkuhler, A.E. Mark, S. Reich, R.D. Skeel (Eds.), Computational Molecular Dynamics: Challenges, Methods, Ideas, Lecture Notes in Computational Science and Engineering, vol. 4, Springer, Heidelberg, 1999. [11] P.G. Drazin, Solitons, Lecture Note Series 85, London Mathematical Society, Cambridge University Press, Cambridge, 1983. [12] A. Durán, J.M. Sanz-Serna, The numerical integration of relative equilibrium solutions. The nonlinear Schrödinger equation, IMA J. Numer. Anal. 20 (2000) 235–261. [13] E. Hairer, C. Lubich, On the life-span of backward error analysis, Numer. Math. 76 (1997) 441–462. [14] J.E. Marsden, G.P. Patrick, S. Shkoller, Multi-symplectic geometry, variational integrators, and nonlinear PDEs, Comm. Math. Phys. 199 (1998) 351–395. [15] B.E. Moore, A modified equations approach for multisymplectic integration methods, Ph.D. Thesis, University of Surrey, in preparation.
[16] B.E. Moore, S. Reich, Backward error analysis for multisymplectic integration methods, submitted for publication. [17] P.J. Olver, Applications of Lie Groups to Differential Equations, Springer, New York, 1986. [18] S. Reich, Backward error analysis for numerical integrators, SIAM J. Numer. Anal. 36 (1999) 1549–1570. [19] S. Reich, Multi-symplectic Runge–Kutta collocation methods for Hamiltonian wave equations, J. Comput. Phys. 157 (2000) 473–499. [20] J.M. Sanz-Serna, M.P. Calvo, Numerical Hamiltonian Problems, Chapman & Hall, London, 1994. [21] G.B. Whitham, Linear and Nonlinear Waves, Wiley/ Interscience, New York, 1974.
Brian E. Moore studied at Colorado Christian University from 1993 until 1997 when he was awarded a Bachelor of Science degree in Mathematics. In 1999, he received his Master of Science degree in mathematical and computer sciences at Colorado School of Mines. He is currently in the third year of his PhD work in applied mathematics, which is sponsored by the University of Surrey, and he is working under the supervision of Prof. Sebastian Reich at Imperial College, London. His research interests include dynamical systems, numerical analysis, and nonlinear PDEs.
Sebastian Reich studied at the Technical University Dresden from 1982 until 1986, where he received his master degrees in electrical engineering and mathematics. He earned his PhD in 1990. His research interests include numerical analysis, dynamical systems, molecular dynamics, and geophysical fluid dynamics. He has held research positions at the Karl-Weierstrass Institute and the Konrad-Zuse-Zentrum in Berlin, the University of British Columbia and Simon Fraser University in Vancouver, and the University of Surrey, Guildford. He is currently Professor in computational and mathematical modeling at Imperial College, London.
Multi-symplectic integration methods for Hamiltonian PDEs Brian E. Moore a,∗,1 , Sebastian Reich b a
Department of Mathematics and Statistics, University of Surrey, Guildford, Surrey GU2 7XH, UK b Department of Mathematics, Imperial College, London SW2 7BZ, UK
Abstract Recent results on numerical integration methods that exactly preserve the symplectic structure in both space and time for Hamiltonian PDEs are discussed. The Preissman box scheme is considered as an example, and it is shown that the method exactly preserves a multi-symplectic conservation law and any conservation law related to linear symmetries of the PDE. Local energy and momentum are not, in general, conserved exactly, but semi-discrete versions of these conservation laws are. Then, using Taylor series expansions, one obtains a modified multi-symplectic PDE and modified conservation laws that are preserved to higher order. These results are applied to the nonlinear Schrödinger (NLS) equation and the sine-Gordon equation in relation to the numerical approximation of solitary wave solutions. © 2002 Elsevier Science B.V. All rights reserved. Keywords: Multi-symplectic PDEs; Preissman box scheme; Conservation laws; Modified equations; Backward error analysis
1. Introduction When approximating the solution of a differential equation using a numerical integrator, it is often desirable to preserve as many properties of the exact evolution equation as possible. In particular, it is advantageous to preserve the symplectic structure of Hamiltonian systems 0 d Id Jyt = ∇y H(y) with J = , (1) −Id 0d where y(t) ∈ R2d and H : R2d → R is the Hamiltonian (energy) of the system. Numerical schemes that have this property, known as symplectic integrators, have proven to be both accurate and efficient in the ∗ Corresponding author. E-mail addresses: [email protected] (B.E. Moore), [email protected] (S. Reich). 1 This work was accomplished while the author was visiting Imperial College, London.
long-time approximation of solutions to Hamiltonian ODEs [3,10,13,18,20]. In addition, symplectic integrators lead to a useful interpretation of backward error analysis, in which a modified equation, the equation solved by the numerical scheme to higher order than the original equation, is used to understand the discretization error induced by the numerical scheme. Recently, the idea of symplectic integration has been extended in a novel way to Hamiltonian PDEs [8,14], by exploiting the multi-symplectic structure of the PDE Kzt + Lzx = ∇z S(z),
(2)
where K and L are skew-symmetric, z(x, t) is the vector of state variables, and S(z) is a smooth function. Specifically, Bridges and Reich [8] apply a symplectic integrator to each independent variable to get a multi-symplectic integrator which preserves the multi-symplectic structure of (2). This approach has paved the way to a useful backward error analysis for multi-symplectic integrators, such that a modified
0167-739X/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 7 - 7 3 9 X ( 0 2 ) 0 0 1 6 6 - 8
396
B.E. Moore, S. Reich / Future Generation Computer Systems 19 (2003) 395–402
multi-symplectic PDE is satisfied by the numerical solution to higher order in both space and time [15,16]. The modified PDE can also be used to derive modified conservation laws of energy and momentum, leading to a better understanding of the numerical scheme. The focus of the present paper is on the multisymplectic Preissman box scheme, which is the space–time version of the implicit midpoint scheme for ODEs. This scheme has been considered in the context of multi-symplectic PDEs in [8,19], and it is commonly used in the area of hydraulics [1,2]. It is shown that this scheme preserves a conservation law, which is directly related to linear symmetries associated with (2) through Noether theory, as well as a multi-symplectic conservation law, exactly. Backward error analysis is also used to show that the modified PDE satisfies these conservation laws and the conservation laws of energy and momentum to higher order. The outline is as follows. We first introduce the conservation laws satisfied by (2) in Section 2, and discuss the box scheme as well as its conservation properties in Section 3. Then, Section 4 considers the modified equations and conservation laws obtained through backward error analysis. The results are then applied in Section 5 to a model problem known as the nonlinear Schrödinger (NLS) equation and soliton solutions of the sine-Gordon equation.
2. Conservation laws There are several properties of (2) that take the form of conservation laws and should be considered when using a numerical integrator to approximate solutions. First, according to Bridges [4,5], a multi-symplectic conservation law ∂t ω + ∂x κ = 0,
where ω = dz ∧ K dz and
κ = dz ∧ L dz is derived directly from (2). Here, we use wedge product notation such that dz denotes the vector of differentials, and let ∂t and ∂x denote differentiation with respect to time and space. One can restrict this conservation law to the space of solutions by making the ansatz dz = zt dt + zx dx.
This yields the conservation law ∂t (zTt Kzx ) + ∂x (zTt Lzx ) = 0.
(3)
Following the analysis in [4,5], we can also derive local conservation laws of energy and momentum. Using the time invariance of (2), an energy conservation law can easily be derived by taking the product of zTt with (2). Then we obtain the energy conservation law ∂t E + ∂x F = 0,
where E = S(z) + zTx Lz/2 and
F = −zTt Lz/2 are known, respectively, as the energy density and the energy flux. Similarly, the spatial invariance of (2) can be used to take the product of zTx with (2), which gives the momentum conservation law ∂x G + ∂t I = 0,
where G = S(z) + zTt Kz/2 and
I = −zTx Kz/2. It is interesting to note that the conservation law (3) is equivalent to the one obtained from differentiation of the energy and momentum conservation laws, i.e. ∂x (Et + Fx ) − ∂t (Gx + It ) = 0.
(4)
According to Bridges [4], there may also be additional conservation laws for (2) related to linear symmetries. Essentially, this follows from Noether theory and derivation of the multi-symplectic PDE (2) from a Lagrangian functional 1 L = Ldt dx for L = zT (Kzt + Lzx ) − S(z). 2 (5) More specifically, take a linear one-parameter family of linear coordinate transformations [17] given by the group action Gε (z) = eεA z,
(6)
which is chosen in such a way that it is symplectic with respect to both ω and κ. Symplecticity is equivalent to AT K + KA = 0
and
AT L + LA = 0.
(7)
Using these two identities, the invariance of a Lagrangian under such a change of variables leads to ∂L 0= = S (z)Az dt dx ∂ε ε=0
B.E. Moore, S. Reich / Future Generation Computer Systems 19 (2003) 395–402
and we obtain the invariance condition
S (z)Az = 0.
397
where (8)
Direct application of (8) to the multi-symplectic PDE formulation (2) yields (Az)T Kzt + (Az)T Lzx = (Az)T ∇z S(z) = 0
zn+1/2,i+1/2 = 41 (zn+1,i+1 + zn+1,i + zn,i+1 + zn,i ). We will call this scheme a multi-symplectic integrator because it preserves exactly a multi-symplectic conservation law [8]
and using (7), this can be written as a conservation law
n+1/2,i n,i+1/2 δ+ + δ+ =0 t ω xκ
∂t (zT KAz) + ∂x (zT LAz) = 0.
for
(9)
Notice here that each of these conservation laws is preserved locally, and these formulations are made possible by the multi-symplectic structure of the PDE. This is in contrast to a Hamiltonian formulation which yields global conservation of densities and neglects the fluxes. We make the distinction here because the preservation of a local conservation law numerically is a stronger result than the global conservation that has been achieved in the past, and these are the issues discussed as we consider a numerical method.
ωn+1/2,i = dzn+1/2,i ∧ K dzn+1/2,i and κn,i+1/2 = dzn,i+1/2 ∧ K dzn,i+1/2 . Eq. (10) can also be derived using a discrete variational principle to get ∂ L = 0, ∂zn,i
L=
n
Ln+1/2,i+1/2
i
for the approximation 3. A multi-symplectic box scheme
n+1/2,i Ln+1/2,i+1/2 = 21 (zn+1/2,i+1/2 )T Kδ+ t z
We now present a centered scheme and discuss the preservation of these conservation laws. First, introduce the notation zn,i , which denotes the numerical approximation of z(xn , ti ) for n = 0, 1, 2, . . . , N and i = 0, 1, 2, . . . , T , where N is the number of grid points and T is the number of time steps. We also set x = xn+1 −xn and t = ti+1 −ti . In addition, define the midpoints zn+1/2,i = 21 (zn+1,i + zn,i ) n,i+1/2
z
=
1 n,i+1 2 (z
and
n,i
+z )
and the forward difference approximations zn+1,i − zn,i and x zn,i+1 − zn,i n,i zt ≈ δ+ . z = t t Using the implicit midpoint scheme [20] to discretize (2) in space and time separately yields the multi-symplectic Preissman box scheme [8]
n,i zx ≈ δ+ = xz
n+1/2,i n,i+1/2 Kδ+ + Lδ+ = ∇z S(zn+1/2,i+1/2 ), t z xz
(10)
n,i+1/2 + 21 (zn+1/2,i+1/2 )T Lδ+ xz
−S(zn+1/2,i+1/2 ) of (5). This is a generalization of the discrete multi-symplectic formulation of Marsden et al. [14] for first-order field theories. Provided the PDE (2) is invariant under an appropriate linear symmetry, the Preissman box scheme also satisfies a discrete version of the conservation law (9) exactly. To make this clear, multiply (10) by (Azn+1/2,i+1/2 )T to get n+1/2,i (zn+1/2,i+1/2 )T KAδ+ t z
n,i+1/2 +(zn+1/2,i+1/2 )T LAδ+ = 0, xz
where we make use of the symplectic and invariance conditions (7) and (8). Then it is easy to show that this is equivalent to the discrete conservation law n+1/2,i T δ+ ) KAzn+1/2,i ) t ((z
n,i+1/2 T +δ+ ) LAzn,i+1/2 ) = 0. x ((z
For linear PDEs, where S(z) = (zT Bz)/2 for B symmetric, this scheme also preserves fully discrete
398
B.E. Moore, S. Reich / Future Generation Computer Systems 19 (2003) 395–402
energy and momentum conservation laws exactly [8]. However, for nonlinear systems in general, exact conservation of discrete energy and momentum conservation laws cannot be guaranteed, but certain semi-discretized energy and momentum conservation laws are still preserved exactly. For example, similar to the analysis of [16], we can take the product of n+1/2 T (zt ) with the spatially discrete equation n+1/2
Kzt
n n+1/2 + Lδ+ ) x z = ∇S(z
to get n+1/2 T
n n+1/2 ) Lδ+ ), x z = ∂t S(z
(zt
n+1/2 T n+1/2 ) Kzt
because K skew-symmetric implies (zt = 0. Now, since n T n −δ+ x ((z ) Lzt ) n+1/2 T
n n+1/2 T n ) Lδ+ ) Lδ+ x z − (z x zt ,
= (zt
concept of modified equations through backward error analysis.
4. The modified equations Based on the ideas of [16], we can perform a useful backward error analysis for the discretization (10), such that the modified equations are also multi-symplectic and satisfy modified versions of the conservation laws in Section 2. Using Taylor series expansions about ti+1/2 we find that δ+ t z(ti ) = zt (ti+1/2 ) +
and an analogous expansion holds for δ+ x z(xn ), which implies that the scheme is second-order in both space and time. Truncating these expansions and substituting into (10) yields the modified system of equations Kzt +
we have
2 4 1 24 t zttt (ti+1/2 ) + O(t )
2 1 24 t Kzttt
+ Lzx +
2 1 24 x Lzxxx
= ∇z S(z), (11)
n+1/2 T
n n+1/2 T n 1 ) Lδ+ ) Lδ+ x z = 2 ∂t ((z xz )
(zt
n T n − 21 δ+ x ((z ) Lzt ).
This leads to the spatially discrete energy conservation law n ∂t En+1/2 + δ+ xF =0
˜ zt + L˜ ˜ zx = ∇z˜ S(˜ ˜ z) K˜
with E
n+1/2 n
F =
= S(z
n+1/2
)−
which can be viewed as a generalized (higher order) multi-symplectic PDE and is satisfied by the numerical solution up to O(t 4 + x4 ). Eq. (11) can also be written in the form of a standard (first-order) multi-symplectic PDE
n 1 n+1/2 T ) Lδ+ xz 2 (z
for z˜ = [z, p, q, r, s]T and and
n 1 n T 2 (z ) Lzt .
A semi-discrete momentum conservation law is derived in the same way by reversing the roles of space and time, i.e. discretize the PDE in time only and take i+1/2 T the product of the resulting equation with (zx ) . Thus, a semi-discrete momentum conservation law i ∂x Gi+1/2 + δ+ t I =0
with
˜ K=
is satisfied for i Gi+1/2 = S(zi+1/2 ) − 21 (zi+1/2 )T Kδ+ t z
S˜ = S +
and
I i = 21 (zi )T Kzix . To understand the energy and momentum conservation of fully discretized equations we must resort to the
L˜ =
2 T 2 T 1 1 24 t q Kp + 24 x s Lr
K
0
1 2 24 t K
0
0
1 − 24 t 2 K
0
0
1 2 24 t K
0
0
0
0
0
0
0
0
L
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1 − 24 x2 L
1 2 24 x L
0
0
0
0
0 0 0, 0 0 0 0 1 2 24 x L 0 0 , 0 0
B.E. Moore, S. Reich / Future Generation Computer Systems 19 (2003) 395–402
399
where p = zt , q = pt , r = zx and s = rx . The modified conservation laws of energy and momentum are easily found using this modified multi-symplectic PDE or (11) directly. The modified form of (3) can then be computed from these modified energy and momentum conservation laws using (4). The modified PDE (11) can also be derived from a modified Lagrangian density
used to investigate computational solitary wave solutions for the sine-Gordon equation.
˜ 1 = 1 (zT Kzt − L 2
where ψ is complex-valued. Letting ψ = v + iw, we can rewrite this as the Hamiltonian system of equations
2 T 1 24 t zt Kztt
+ zT Lzx
1 − 24 x2 zTx Lzxx ) − S(z).
In general, the modified Lagrangian density for any number of modifications ρ, can be written as ˜ρ = L
ρ j=0
+
(−1)j t 2j (2)2j+1 (2j + 1)!
ρ j=0
j
j+1
(∂t z)T K∂t
z
(−1)j x2j (∂j z)T L∂xj+1 z − S(z), (2)2j+1 (2j + 1)! x
j
where ∂η denotes the jth derivative with respect to η. Given a linear symmetry, this modified Lagrangian density is also invariant under the transformation (6). For example, the associated conservation law for ρ = 1 is easily obtained by multiplying the modified equation (11) by (Az)T , which yields zT KAzt +
2 T T 1 24 t z KAzttt + z LAzx 1 + 24 x2 zT LAzxxx = 0.
Then the identity zT KAzttt = ∂t (zT KAztt − 21 zTt KAzt )
5.1. The nonlinear Schrödinger equation Consider the NLS equation iψt + ψxx + V (|ψ|2 )ψ = 0,
vt = −wxx − wV (v2 + w2 ), wt = vxx + vV (v2 + w2 ). Defining ψx = σ + iφ, this system of equations can also be written as the multi-symplectic PDE (2) where z = [v, w, σ, φ]T J 02 02 −I2 K= , L= 02 02 I2 0 2 and S(z) = 21 [σ 2 + φ2 + V(v2 + w2 )], where J is given in (1) with d = 1. This system is invariant under the action of the one-parameter group of rotations SO(2) [7], given by Gθ (z) = Rθ z, for cos θ − sin θ 0 0 cos θ 0 0 sin θ . Rθ = 0 cos θ − sin θ 0 0
sin θ
0
cos θ
and a similar identity for zT LAzxxx , imply the modified conservation law
Now, define the matrix A such that
0 = ∂t [zT KAzt +
Az =
2 T 2 T 1 1 24 t z KAztt − 48 t zt KAzt ] 1 1 +∂x [zT LAzx + 24 x2 zT LAzxx − 48 x2 zTx LAzx ].
5. Applications We can now apply these results to some model problems. We first consider the NLS equation as a Hamiltonian PDE with a rotational symmetry. The second example outlines how backward error analysis can be
d Gθ (z)|θ=0 = [−w, v, −φ, σ]T . dθ
Then zT KAz = w2 + v2
and
zT LAz = 2(vφ − wσ)
satisfy the conservation law (9). Now it is straight forward to derive the results of the previous sections for a Preissman box scheme discretization of this equation, including discrete conservation laws and modified equations with their conservation laws.
400
B.E. Moore, S. Reich / Future Generation Computer Systems 19 (2003) 395–402
5.2. Soliton solutions The NLS equation has a family of solitary wave solutions as do several other equations that can be written in the from (2), such as the KdV equation, Boussinesq equations, and nonlinear wave equations (cf. [4–7,21]). In addition, a thorough analysis of the error growth due to time discretization for these solutions has been performed for the KdV equation [9] and the NLS equation [12]. Specifically, Durán and Sanz-Serna [12] have shown that the implicit midpoint time discretization has better error propagation mechanisms than other non-conservative schemes. In particular, they have shown that the time-discrete solutions of the NLS equation, which have been initialized with an exact soliton profile, are made up of a modified solitary wave, that have a different wave speed and amplitude than the exact soliton, as well as a complementary term that grows linearly in time and higher order terms. Using the results of the previous sections, it may be possible to show that the numerical solutions of the fully discretized NLS equation initialized with an exact soliton profile are made up of a modified soliton, which has corrections in both space and time, along with the complementary error and higher order terms. From a slightly different perspective, one may be able to show that the modified equations have solitary wave solutions. Then, initializing the numerical scheme with the new modified soliton profile would yield numerical solutions that are preserved to higher order in both space and time. Results in that spirit have been achieved concerning soliton solutions of the sine-Gordon equation utt = uxx − sin u,
(12)
which can be written in the from (2) with z = [u, v, w]T
0 −1 0 0 K = 1 0 0, L = 0 0 0 0 1 1 1 S(z) = v2 − w2 − cos u. 2 2
0 0 0
−1 0 0
and
It is well known that there is a family of solitary wave
solutions for (12) given by
±(x − ct − ξ0 ) , u(x, t) = 4 tan −1 ±exp (1 − c2 )1/2
(13)
where |c| < 1 is the wave speed, and ξ0 is given by the initial conditions [11]. For our purposes we take ξ0 = 0. Using the explicit first-order in space and time multi-symplectic Euler box scheme [16] to discretize gives n,i n,i δ+ + δ+ = − sin un,i , t v xw n,i δ− = −wn,i , xu
n,i δ− = vn,i , t u
(14)
where we have used the backward differences n,i = δ− xu
un,i − un−1,i un,i − un,i−1 n,i and δ− . = t u x t
Notice that eliminating v and w gives the familiar second-order scheme δ2t un,i = δ2x un,i − sin un,i ,
(15)
− where δ2η = δ+ η δη for η = x, t. Hence, the numerical n,i solution u is accurate to second-order in both t and x, but using the scheme (14) yields first-order error because v and w are only approximated to first-order. To make this more clear, backward error analysis [16] can be used for (14) to get the modified equations
vt + 21 tvtt + wx + 21 xwxx = − sin u,
(16)
ut − 21 tutt = v,
(17)
ux − 21 xuxx = −w.
(18)
Since the numerical solution is already second-order accurate in u, only (17) and (18) need to be verified along the numerical solution to check for second-order accuracy with respect to the modified equations. Computing higher order corrections gets increasingly more complicated. However, the modified equation for (15), given by utt +
2 1 12 t utttt
= uxx + 21 x2 uxxxx − sin u
(19)
is easily obtained using Taylor expansions and is satisfied by the numerical solution to fourth-order. Now it is reasonable to ask if this modified system of equations has soliton solutions, which can be answered by
B.E. Moore, S. Reich / Future Generation Computer Systems 19 (2003) 395–402
401
Fig. 1. Error plotted against t for the original initial condition (dashed), and the modified initial condition for one modification (chain) and three modifications (solid).
setting u(x, t) = u(x − ct) = u(ξ) in (19) and finding solutions of the boundary value problem (c2 − 1)u −
2 1 12 (x
− c4 t 2 )u = − sin u.
(20)
Though it is not obvious how to solve this equation for a heteroclinic connection, we were able to find such a connection to high accuracy numerically. This approximation can be used to initialize the scheme (15) and this in turn yields a numerical solution for the PDE that is accurate to fourth-order in u with respect to both t and x. Using log–log scale in Fig. 1, we plot the solution error against t. For these results we have kept the ratio t/x fixed, so change in t directly corresponds to change in x. The dashed line shows that the numerical solution, obtained by initializing (14) with un,0 = u(xn , 0), vn,0 = ut (xn , 0) and
wn,0 = −ux (xn , 0) for u given in (13), is first-order accurate. The chain line shows that the same initial condition for (14) with v and w modified according to (17) and (18) yields second-order accuracy in both t and x. Finally, using a numerically computed solution of (20) as the initial condition for the scheme (15), the solid line shows approximately fourth-order convergence. The accuracy here is limited by the fact that the reference soliton solution had to be computed numerically, and this is a delicate procedure due to the heteroclinic nature of the solution.
Acknowledgements Sebastian Reich gratefully acknowledges partial financial support by European Commission funding for
402
B.E. Moore, S. Reich / Future Generation Computer Systems 19 (2003) 395–402
the Research Training Network “Mechanics and Symmetry in Europe”. References [1] M.B. Abbott, Computational Hydraulics, Pitman, London, 1979. [2] M.B. Abbott, D.R. Basco, Computational Fluid Dynamics, Longman Scientific, Harlow, 1989. [3] G. Benettin, A. Giorgilli, On the Hamiltonian interpolation of near to the identity symplectic mappings with application to symplectic integration algorithms, J. Stat. Phys. 74 (1994) 1117–1143. [4] T.J. Bridges, Multi-symplectic structures and wave propagation, Math. Proc. Camb. Phil. Soc. 121 (1997) 147–190. [5] T.J. Bridges, A geometric formulation of the conservation of wave action and its implications for signature and the classification of instabilities, Proc. Roy. Soc. London Ser. A 453 (1997) 1365–1395. [6] T.J. Bridges, G. Derks, Unstable eigenvalues and the linearization about solitary waves and fronts with symmetry, Proc. Roy. Soc. London Ser. A 455 (1999) 2427–2469. [7] T.J. Bridges, G. Derks, The symplectic Evans matrix, and the instability of solitary waves and fronts with symmetry, Arch. Rat. Mech. Anal. 156 (2001) 1–87. [8] T.J. Bridges, S. Reich, Multi-symplectic integrators: numerical schemes for Hamiltonian PDEs that conserve symplecticity, Phys. Lett. A 284 (2001) 184–193. [9] J. De Frutos, J.M. Sanz-Serna, Accuracy and conservation properties in numerical integration: the case of the Korteweg– de Vries equation, Numer. Math. 75 (1997) 421–445. [10] P. Deuflhard, J. Hermans, B. Leimkuhler, A.E. Mark, S. Reich, R.D. Skeel (Eds.), Computational Molecular Dynamics: Challenges, Methods, Ideas, Lecture Notes in Computational Science and Engineering, vol. 4, Springer, Heidelberg, 1999. [11] P.G. Drazin, Solitons, Lecture Note Series 85, London Mathematical Society, Cambridge University Press, Cambridge, 1983. [12] A. Durán, J.M. Sanz-Serna, The numerical integration of relative equilibrium solutions. The nonlinear Schrödinger equation, IMA J. Numer. Anal. 20 (2000) 235–261. [13] E. Hairer, C. Lubich, On the life-span of backward error analysis, Numer. Math. 76 (1997) 441–462. [14] J.E. Marsden, G.P. Patrick, S. Shkoller, Multi-symplectic geometry, variational integrators, and nonlinear PDEs, Comm. Math. Phys. 199 (1998) 351–395. [15] B.E. Moore, A modified equations approach for multisymplectic integration methods, Ph.D. Thesis, University of Surrey, in preparation.
[16] B.E. Moore, S. Reich, Backward error analysis for multisymplectic integration methods, submitted for publication. [17] P.J. Olver, Applications of Lie Groups to Differential Equations, Springer, New York, 1986. [18] S. Reich, Backward error analysis for numerical integrators, SIAM J. Numer. Anal. 36 (1999) 1549–1570. [19] S. Reich, Multi-symplectic Runge–Kutta collocation methods for Hamiltonian wave equations, J. Comput. Phys. 157 (2000) 473–499. [20] J.M. Sanz-Serna, M.P. Calvo, Numerical Hamiltonian Problems, Chapman & Hall, London, 1994. [21] G.B. Whitham, Linear and Nonlinear Waves, Wiley/ Interscience, New York, 1974.
Brian E. Moore studied at Colorado Christian University from 1993 until 1997 when he was awarded a Bachelor of Science degree in Mathematics. In 1999, he received his Master of Science degree in mathematical and computer sciences at Colorado School of Mines. He is currently in the third year of his PhD work in applied mathematics, which is sponsored by the University of Surrey, and he is working under the supervision of Prof. Sebastian Reich at Imperial College, London. His research interests include dynamical systems, numerical analysis, and nonlinear PDEs.
Sebastian Reich studied at the Technical University Dresden from 1982 until 1986, where he received his master degrees in electrical engineering and mathematics. He earned his PhD in 1990. His research interests include numerical analysis, dynamical systems, molecular dynamics, and geophysical fluid dynamics. He has held research positions at the Karl-Weierstrass Institute and the Konrad-Zuse-Zentrum in Berlin, the University of British Columbia and Simon Fraser University in Vancouver, and the University of Surrey, Guildford. He is currently Professor in computational and mathematical modeling at Imperial College, London.