Combining finite element space-discretizations with symplectic time-marching schemes for linear Hamiltonian systems

We provide a short introduction to the devising of a special type of methods for numerically approximating the solution of Hamiltonian partial differential equations. These methods use Galerkin space-discretizations which result in a system of ODEs displaying a discrete version of the Hamiltonian structure of the original system. The resulting system of ODEs is then discretized by a symplectic time-marching method. This combination results in high-order accurate, fully discrete methods which can preserve the invariants of the Hamiltonian defining the ODE system. We restrict our attention to linear Hamiltonian systems, as the main results can be obtained easily and directly, and are applicable to many Hamiltonian systems of practical interest including acoustics, elastodynamics, and electromagnetism. After a brief description of the Hamiltonian systems of our interest, we provide a brief introduction to symplectic time-marching methods for linear systems of ODEs which does not require any background on the subject. We describe then the case in which finite-difference space-discretizations are used and focus on the popular Yee scheme (1966) for electromagnetism. Finally, we consider the case of finite-element space discretizations. The emphasis is placed on the conservation properties of the fully discrete schemes. We end by describing ongoing work.

1. Introduction

We present a short introduction to the devising of a particular type of numerical methods for linear Hamiltonian systems. The distinctive feature of these methods is that they are obtained by combining finite element methods for the space-discretization with symplectic methods for the time-discretization. The finite element space discretization of the Hamiltonian system is devised in such a way that it results in a system of ODEs which is also Hamiltonian. This guarantees that, when the system is discretized by a symplectic method, the resulting fully discrete method can conserve the linear and quadratic time-invariants of the original system of ODEs. This is a highly-sought property, especially for long-time simulations.

For a description of the pioneering work and of the development of these methods, the reader can see the Introductions in our papers published in 2017 [1], in 2021 [2], and in 2022 [3] and the references therein. In those papers, other approaches to achieve energy-conserving schemes are also briefly reviewed. Particularly impressive are the schemes proposed in 2019 by Fu and Shu [4] in which the authors construct energy-conserving discontinuous Galerkin (DG) methods for general linear, symmetric hyperbolic systems.

When we began to study this type of numerical methods for linear Hamiltonian systems, we were impressed to see that the methods displayed arbitrary high-order accuracy and preserved (in time) very well the energy and other quantities of physical relevance. Although we were familiar with finite element space discretizations, we were unaware of the relevance and properties of the symplectic time-marching methods. In a way, the paper we present here is the paper we would have liked to read as we plunged into this topic. It is written as an introduction to the subject for numerical analysts of PDEs which, like us, were not familiar with symplectic methods.

We restrict ourselves to the linear case for two reasons. The first is that we have studied those numerical methods for Hamiltonian formulations of acoustic wave propagation, elastodynamics and the Maxwell's equations. So, we do have some experience to share for those systems. The second is that the theory of symplectic time-marching can be covered quite easily and that most of its results also hold in the nonlinear case, although they are not that easy to obtain.

Let us sketch a rough table of contents:

• Hamiltonian and Poisson Systems. (Section 2) We define a general Hamiltonian (and its simple extension, the Poisson) systems, of ODEs or PDEs. We then restrict ourselves to linear systems, give several examples, and discuss the conservation laws associated to them. For the Hamiltonian and Poisson systems, we took material from the books that appeared in 1993 by Olver [5], in 1999 by Marsden and Ratiu [6], in 2004 by Leimkuhler and Reich [7], and in 2010 by Feng and Qin [8]. For the conservation laws, we took the information gathered in our papers on the acoustic wave equation in 2017 [1] elastodynamics in 2021 [2], and electromagnetism in 2022 [3, 9].

• Symplectic and Poisson time-marching methods. (Section 3) We describe the symplectic (and their simple, but useful, extensions to the Poisson) methods for linear Hamiltonian ODEs with special emphasis on how they approximate their linear and quadratic invariants. We took material from the pioneering work of the late 1980's by Feng et al. [10–12], from the popular 1992 review by Sans-Serna [13], from the 2006 book by Hairer et al. [14], and form the 2010 book by Feng and Qin [8].

• Finite Difference space discretizations. (Section 4) We show how to discretize in space the linear Hamiltonian systems by using finite difference methods. In particular, we recast as one of such schemes the popular Yee scheme for electromagnetism proposed in 1966 [15], that is, more than 17 years before the appearance of the first symplectic method. We took material from the work done in 1988 by Ge and Feng [10] and in 1993 by McLachlan [16].

• Finite Element space discretizations. (Section 5) We describe space discretizations based on finite element methods and combine them with symplectic time-marching schemes. We discuss the corresponding conservation properties and show a numerical experiment validating a couple of them. We took material from the work done in 2005 by Groß et al. [17], in 2008 by Xu et al. [18], and by Kirby and Kieu [19], as well as from our own work in 2017 [1], 2021 [2], and 2022 [3, 9].

• Ongoing work. (Section 6) We end by briefly commenting on some ongoing work.

2. Hamiltonian and Poisson systems

2.1. Canonical Hamiltonian systems and their generalization

2.1.1. The canonical nonlinear Hamiltonian systems

A canonical Hamiltonian system is defined on (the phase-space) ℝ²ⁿ in terms of the smooth Hamiltonian functional ${H}:{ℝ}^{2n}\to ℝ$ and a structure matrix J by the system of equations

$\begin{array}{l}\frac{d}{dt}\text{u}=\text{J}{▽}_{\text{u}}{H}(\text{u}),\end{array}$

where¹

$\begin{array}{ll}\text{J}:=\left[\begin{array}{cc}0& -\text{Id}\\ \text{Id}& 0\end{array}\right].& (1)\end{array}$

The Hamiltonian has the remarkable property of remaining constant on the orbits of the system, t ↦ u(t). Indeed,

$\frac{d}{dt}{H}(\text{u})={(▽}_{\text{u}}{H}(\text{u})){}^{\top }\frac{d}{dt}\text{u}={(▽}_{\text{u}}{H}(\text{u})){}^{\top }\text{ J }{(▽}_{\text{u}}{H}(\text{u}))=0,$

since the matrix J is anti-symmetric. This computation suggests an immediate extension of this Property to other functionals. If $\text{u}\mapsto {F}(\text{u})$ is any differentiable functional, we have that

$\frac{d}{dt}{F}(\text{u})={(▽}_{\text{u}}{F}(\text{u})){}^{\top }\frac{d}{dt}\text{u}={(▽}_{\text{u}}{F}(\text{u})){}^{\top }\text{ J }{(▽}_{\text{u}}{H}(\text{u}))=:\left\{{F},{H}\right\},$

and we can see that ${F}(\text{u})$ remains constant on the orbits of the system, t ↦ u(t), if and only if the Poisson bracket, $\left\{{F},{H}\right\}$, is zero. Such quantities are usually called first integrals of the system. They are also called invariants or conserved quantities.²

2.1.2. A generalization

This definition can be generalized by taking a phase space of arbitrary dimension, by allowing J to be an arbitrary anti-symmetric matrix, or by letting it depend on the phase-space variable u, that is, by using J: = J(u). In these cases, the system is called a Poisson system.

A Poisson dynamical system is, see the 1993 books [5, 16], the triplet $({M},{H},\{·,·\})$, where ${M}$ is the phase-space, ${H}$ is the Hamiltonian functional, and {·,·} is the Poisson bracket. The system can then be expressed as

$\frac{d}{dt}\text{u}=\left\{\text{u},{H}\right\}.$

The Poisson bracket {·,·} is a general operator with which we can describe the time-evolution of any given smooth function ${F}:{ℝ}^d\to ℝ$. This operation is defined for a pair of real-smooth functions on the phase-space satisfying: (i) bilinearity, (ii) anti-symmetry, (iii) Jacobi identity:

$\left\{\left\{{F},{G}\right\},{H}\right\}+\left\{\left\{{H},{F}\right\},{G}\right\}+\left\{\left\{{G},{H}\right\},{F}\right\}=0,$

and, (iv) the Leibniz' rule:

$\left\{{F},{G}·{H}\right\}=\left\{{F},{G}\right\}·{H}+{G}·\left\{{F},{H}\right\},$

for ${F},{G},{H}:{ℝ}^d\to ℝ$. For a constant anti-symmetric matrix J, usually called the structure matrix, these four properties are trivially satisfied by the induced Poisson bracket, see the 1993 book by Olver [5, Corollary 7.5].

2.2. Definition of linear Hamiltonian and Poisson systems

As we restrict ourselves to linear Hamiltonian and Poisson systems, let us formally define them. All the definitions are in terms of the structure matrix J (which is anti-symmetric), but, to alleviate the notation, we omit its mention.

Definition 2.1 (Hamiltonian systems). We say that the linear system

$\frac{d}{dt}\text{u}=\text{J}{▽}_{\text{u}}{H}(\text{u})$

is Hamiltonian if the structure matrix J (which is anti-symmetric) is invertible and the Hamiltonian is a quadratic form ${H}(\text{v}):=\frac{1}{2}{\text{v}}^{\top }\text{H}\text{v}$ for some symmetric matrix H.

Definition 2.2 (Poisson systems). We say that the linear system

$\frac{d}{dt}\text{u}=\text{J}{▽}_{\text{u}}{H}(\text{u})$

is a Poisson system if J is a structure matrix (which is anti-symmetric) and the Hamiltonian is a quadratic form ${H}(\text{v}):=\frac{1}{2}{\text{v}}^{\top }\text{H}\text{v}$ for some symmetric matrix H.

We see that a Poisson system is the straightforward generalization of a Hamiltonian system,³ to the case in which the structure matrix is not invertible.⁴

2.3. Examples of linear Hamiltonian systems of ODEs

Let us illustrate and discuss examples of linear Hamiltonian systems and their respective symplectic structure.

Example 1. A textbook example of a linear canonical Hamiltonian system is the system modeling the harmonic oscillator. Its Hamiltonian and structure matrix are

${H}(p,q)=\frac{p^2}{2}+\frac{q^2}{2},\text{J}=\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right].$

This gives rise to the Hamiltonian system of equations

$ṗ=-q,\stackrel{.}{q}=p.$

We know that the restriction of the Hamiltonian to the orbits of this system is constant in time, or, in other words, that the quadratic form ${H}(p,q)=\frac{p^2}{2}+\frac{q^2}{2}$ is a first integral of the system. This implies that the orbits lie on circles. No interesting property can be drawn from a linear first integral of the system since the only one it has is the trivial one.

Example 2. The following is a simple example in which the structure matrix J is not invertible. We work in a three-dimensional space (we take u: = (p, q, r)^⊤) and consider the system associated with the following Hamiltonian and structure matrix:

${H}(\text{u})=\frac{p^2}{2}+\frac{q^2}{2}+\frac{r^2}{2},\text{J}=\left[\begin{array}{ccc}0& -1& 1\\ 1& 0& 0\\ -1& 0& 0\end{array}\right].$

This gives rise to the Poisson system

$\frac{d}{dt}\left[\begin{array}{c}p\\ q\\ r\end{array}\right]=\left[\begin{array}{c}-q+r\\ p\\ -p\end{array}\right].$

Let us examine the conserved quantities of the system. Again, we know that the restriction of the Hamiltonian to the orbits of this system is constant in time, or, in other words, that the quadratic form ${H}(p,q,r)=\frac{p^2}{2}+\frac{q^2}{2}+\frac{r^2}{2}$ is a first integral of the system. This implies that the orbits lie on spheres. Unlike the previous case, this system has one linear first integral,⁵ namely, C(p, q, r): = (0, 1, 1)·(p, q, r), as this reflects the fact that q+r is a constant on the orbits of the system. This has an interesting geometric consequence since it means that the orbits of the system lie on the intersection of a sphere and a plane of normal (0, 1, 1), see Figure 1.

FIGURE 1

Figure 1. Two views of the orbits (orange solid line) of the Poisson system on the phase-space (q, p, r). The orbits are circles (right) lying on planes (left) with normal (0, 1, 1) so that q + r is a constant. The blue sphere represents the constant Hamiltonian value corresponding to ${H}=3^2/2$.

2.4. Hamiltonian PDEs: Definition and examples

To extend the definition of finite-dimensional Hamiltonian and Poisson systems, the triplet $({M},{H},\{·,·\})$, to systems of partial differential equations, we need to introduce a phase-space ${M}$ of smooth functions with suitable boundary conditions (periodic for simplicity), a Hamiltonian functional ${H}$, and a Poisson bracket {·,·}. The definition of the bracket is induced by an anti-symmetric structure operator J, for functionals ${F},{G}$ on the phase-space ${M}$ defined by

$\{ℱ,G\}(\text{u}):={\displaystyle {\int }_{\Omega }\frac{\delta ℱ}{\delta \text{u}}}\text{J}\frac{\delta G}{\delta \text{u}}dx.$

It satisfies bilinearity, anti-symmetry, and the Jacobi identity (Leibniz's rule is dropped). In the definition the operation δ/δu denotes the functional derivative defined for functionals ${F}$ on ${M}$ by

$\begin{array}{l}\frac{\delta {F}}{\delta \text{u}}:{\displaystyle {\int }_{\Omega }\frac{\delta {F}[\text{u}]}{\delta \text{u}}}\delta \text{u}dx=\underset{\epsilon \to 0}{\lim }\frac{{F}[\text{u}+\epsilon \delta \text{u}]-{F}[\text{u}]}{\epsilon }\\ =\frac{d}{d\epsilon }{|}_{\epsilon =0}{F}[\text{u}+\epsilon \delta \text{u}],\end{array}$

for any smooth test function δv in ${M}$. Thus, a Hamiltonian partial differential equation system is given by

$\stackrel{·}{\text{u}}:=\frac{d}{dt}\text{u}=\left\{\text{u},{H}\right\}(\text{u}).$

An important simplification occurs in the case when the structure operator J is a constant, anti-symmetric operator whose definition might include spatial differential operators. In this case, the Jacobi identity is automatically satisfied. See Olver [5, Corollary 7.5].

Let us define and illustrate the Hamiltonian⁶ systems of PDEs we are going to be working with.

Definition 2.3 (Linear Hamiltonian PDE system). We say that the linear system of partial differential equations

$\stackrel{·}{\text{u}}=\left\{\text{u},{H}\right\}(\text{u})=\text{J}\frac{\delta {H}}{\delta \text{u}}.$

is a Hamiltonian system for any constant structure operator J (which is anti-symmetric), inducing the Poisson bracket {·,·}, and the Hamiltonian functional is a quadratic form (meaning that $\partial {H}/\partial \text{v}$ is linear in v).

2.4.1. The scalar wave equation

Let us consider the linear scalar wave equation on a bounded, connected domain Ω ⊂ ℝ^d with smooth boundary ∂Ω

$\begin{array}{ll}\rho ü-▽·(\kappa ▽u)=f,\text{in }\Omega ,t\ge 0,u=0\text{on }\partial \Omega ,& (2)\end{array}$

Where ρ is a scalar-valued function and κ a symmetric, positive definite matrix-valued function. Here the unknown u can be a variety of physical quantities. For example, in acoustics, it represents the deviation from a constant pressure of the medium. In water waves, it represents the distance to the position of rest of the surface of a liquid. Here, we take advantage of its closeness to the vector equations of elastodynamics,⁷ and we formally assume that u has the dimensions of a length. We can then introduce the velocity variable v, and rewrite the equation as a first-order system of equations as follows:

$\begin{array}{ll}\rho \stackrel{·}{u}=\rho v,\rho \stackrel{·}{v}=▽·(\kappa ▽u)+f,\text{in }\Omega ,t\ge 0.& (3)\end{array}$

This system can be written in a Hamiltonian form with Hamiltonian functional, and the canonical structure operator J defined by

${H}\left[u,v\right]={\displaystyle {\int }_{\Omega }}\left(\frac{1}{2}\rho v^2+\frac{1}{2}\kappa ▽u·▽u-fu\right),\text{J}={\rho }^{-1}\left[\begin{array}{cc}0& \text{Id}\\ -\text{Id}& 0\end{array}\right].$

Indeed, observe that the variational derivatives of the Hamiltonian functional are

$\begin{array}{ll}\frac{d}{d\epsilon }{|}_{\epsilon =0}{H}\left[u+\epsilon \delta u,v+\epsilon \delta v\right]& ={\displaystyle {\int }_{\Omega }}(\rho v\delta v-(f+▽·(\kappa ▽u))\delta u),\end{array}$

for δu, δv smooth test functions with compact support on Ω.

A second and equivalent Hamiltonian form can be found as follows. We rewrite the scalar wave equations as a first-order system introducing the vector-valued function q = κ▽u and removing u from the equations, this is known as the mixed formulation

$\begin{array}{ll}\rho \stackrel{·}{v}=▽·\text{q}+f,\stackrel{·}{\text{q}}=\kappa ▽v,\text{in }\Omega ,\forall t\ge 0,v=0\text{on }\partial \Omega .& (4)\end{array}$

This system can also be written in a Hamiltonian form with an equivalent Hamiltonian functional, now in terms of v and q, and the structure operator J defined by

$\begin{array}{l}{H}[v,\text{q}]={\displaystyle {\int }_{\Omega }\left(\frac{1}{2}\rho v^2+\frac{1}{2}{\kappa }^{-1}|\text{q}{|}^2+F·\text{q}\right)},\\ \text{ J}=\left[\begin{array}{cc}0& ({\rho }^{-1}\nabla ·)◦\kappa \\ (\kappa \nabla )◦{\rho }^{-1}& 0\end{array}\right]\end{array}$

The Poisson bracket is then induced by the operator, for functionals F, G

$\left\{F,G\right\}={\displaystyle {\int }_{\Omega }}\left(\frac{\delta F}{\delta \text{q}}\kappa ▽({\rho }^{-1}\frac{\delta G}{\delta v})+\frac{\delta F}{\delta v}{\rho }^{-1}▽·(\kappa \frac{\delta G}{\delta \text{q}})\right)dx.$

Thus,

$\frac{d}{d\epsilon }{|}_{\epsilon =0}\mathscr{H}[v+\epsilon \delta v,\text{q}+\epsilon \delta \text{q}]={\displaystyle {\int }_{\Omega }\left(\rho v\delta v+({\kappa }^{-1}\text{q}-F)·\delta \text{q}\right)},$

where f = ▽·(κF).

In Tables 1, 2, we describe the conservation laws associated with the linear scalar wave equation. See Example 4.36 in Olver [5], and in Walter [32].

TABLE 1

Table 1. Glossary of scalar wave equation quantities.

TABLE 2

Table 2. Conservation laws for the scalar quantity η, $\stackrel{·}{\eta }+▽·f_{\eta }=0$, deduced from the scalar wave equation.

2.4.2. Electromagnetics

Our last example is the Maxwell's equations for the electric and magnetic fields,

$ϵ\stackrel{·}{\text{E}}=▽\times \text{H}-\text{J},\mu \stackrel{·}{\text{H}}=-▽\times \text{E}.$

The Hamiltonian form of these equations corresponds to the Hamiltonian functional and anti-symmetric operator given by

$\mathscr{H}[\text{E},\text{H}]={\displaystyle {\int }_{\Omega }\left(\frac{ϵ}{2}|\text{E}{|}^2+\frac{\mu }{2}|\text{H}{|}^2-{\text{J}}_{\times }·\mu \text{H}\right)},$

$\text{J}=\left[\begin{array}{cc}0& ({ϵ}^{-1}▽\times )◦({\mu }^{-1})\\ -({\mu }^{-1}▽\times )◦({ϵ}^{-1})& 0\end{array}\right],$

Where J_× satisfies ▽ × J_× = J. A second Hamiltonian formulation can be derived by introducing the magnetic vector potential variables A, by μH = ▽ × A, and writing the system as follows

$\stackrel{·}{\text{A}}=-\text{E},ϵ\stackrel{·}{\text{E}}=▽\times ({\mu }^{-1}▽\times \text{A})-\text{J}$

with corresponding Hamiltonian functional and anti-symmetric operator given by

$\mathscr{H}[\text{E},\text{H}]={\displaystyle {\int }_{\Omega }\left(\frac{1}{2}ϵ|\text{E}{|}^2+\frac{\mu }{2}|\text{H}{|}^2-{\text{J}}_{\times }·\mu \text{H}\right)},$

$\text{J}={ϵ}^{-1}\left[\begin{array}{cc}0& \text{Id}\\ -\text{Id}& 0\end{array}\right].$

In Tables 3, 4, we describe the rich set of conservation laws associated to the Maxwell's equations. The first two are associated with the conservation of charge, magnetic and electric. The next three are the classic conservation laws of energy and momenta. Less standard are the last three, which were obtained back in 1964 by Lipkin [21]. See the 2001 paper [33] for a recent update on the classification of conservation laws of Maxwell's equations.

TABLE 3

Table 3. Glossary of electromagnetic quantities, see Sánchez et al. [3].

TABLE 4

Table 4. Conservation laws for the (scalar or vectorial) electromagnetic quantity η, $\stackrel{·}{\eta }+▽·f_{\eta }=S_{\eta }$, deduced from the first two Maxwell's equations.

3. Symplectic and Poisson time-marching methods

As pointed out in the 1992 review by Sanz-Serna [13], the symplectic time-marching methods were devised specifically for the time-discretization of Hamiltonian systems. In this section, we present a simple and direct introduction to the topic, with emphasis on Poisson linear systems, which seems to be new. The approach we take reflects our own path into the subject and is intended for readers that, like us, were utterly and completely oblivious to the relevance of symplectic and Poisson matrices, and were only interested in devising numerical methods for Hamiltonian systems which would capture well their conservation properties. The main objective of this section is to overcome this terrible, terrible shortcoming.

3.1. Symplectic and Poisson matrices

Here, we introduce the symplectic (and their simple extension, the Poisson) matrices which are, as we are soon going to see, deeply related to Hamiltonian (and their simple extension, the Poisson) systems. We use the definitions from the 2006 book by Hairer et al. [14, p. 254].

Definition 3.1 (Symplectic matrices). We say that the invertible matrix E is a symplectic matrix for the structure matrix J (which is anti-symmetric) if J is invertible and if E^⊤J⁻¹E = J⁻¹.

Definition 3.2 (Poisson matrices). We say that the invertible matrix E is a Poisson matrix for the structure matrix J (which is anti-symmetric) if it satisfies EJE^⊤ = J. If, in addition, E^⊤ = Id on the kernel of J, we say that E is a Poisson integrator matrix.

Note that, since

${\text{E}}^{\top }{\text{J}}^{-1}\text{E}={\text{J}}^{-1}\Rightarrow \text{J}={\text{E}}^{-1}\text{J}{\text{E}}^{-\top }\Rightarrow \text{E}\text{J}{\text{E}}^{\top }=\text{J},$

the definition of a Poisson matrix is a straightforward generalization of the definition of a symplectic matrix to the case in which the structure matrix is not necessarily invertible. Note also that the difference between Poisson matrices and Poisson integrator matrices is that the behavior of the transpose of the latter on the kernel of the structure matrix is explicitly specified.

3.1.1. The geometry of Poisson matrices

We next gather what we could call the geometric properties of the Poisson matrices. To state it, we introduce some notation. The kernel of a matrix M is denoted by M⁻¹{0}. Its range is denoted by Mℝ^d. Although we do not assume the structure matrix J to be invertible, we can always define its inverse on part of the space ℝ^d. Indeed, since the restriction of J to Jℝ^d is a one-to-one mapping, we can define J⁻¹ as the inverse of J on Jℝ^d.

Proposition 3.1 (The geometry of Poisson matrices). Let E be a Poisson matrix. Then

$\begin{array}{l}(\text{a})\text{ EJ}{ℝ}^d=\text{J}{ℝ}^d,\\ (\text{b}){\text{E}}^{\top }{\text{J}}^{-1}\text{E}={\text{J}}^{-1}\text{ on J}{ℝ}^d,\\ (\text{c}){\text{E}}^{\top }{\text{J}}^{-1}\left\{0\right\}={\text{J}}^{-1}\left\{0\right\}.\end{array}$

Moreover, if E is also a Poisson integrator matrix, then

$(\text{d})\text{ E}-\text{Id}=\text{JB on }{ℝ}^d,$

for some matrix B: ℝ^d→Jℝ^d.

This result states that a Poisson matrix E is one-to-one on the range of J and behaves as a symplectic matrix⁸ therein, see property (b). Moreover, its transpose is also one-to-one on the kernel of J. If the matrix is also a Poisson integrator, that is, its transpose is the identity on the kernel of J, then E coincides with the identity plus a matrix whose range is included by the range of J.

Proof of Proposition 3.1 Let us prove (a). Since E is a Poisson matrix, we have that EJ = JE^−⊤, and so EJℝ^d⊂Jℝ^d. We also have that E⁻¹J = JE^⊤, and so E⁻¹Jℝ^d⊂Jℝ^d. This proves Property (a).

Let us prove (b). Since the restriction of J to Jℝ^d is a one-to-one mapping, we can define J⁻¹ as the inverse of J on Jℝ^d. Then, since E is a Poisson matrix, on Jℝ^d, we have that J = E⁻¹JE^−⊤ and so, that J⁻¹ = E^⊤J⁻¹E, because EJℝ^d = Jℝ^d by Property (a). This proves Property (b).

Let us now prove (c). Since E is a Poisson matrix, we have that JE^⊤ = E⁻¹J, and so E^⊤J⁻¹{0}⊂J⁻¹{0}. We also have that E^−⊤J⁻¹{0}⊂J⁻¹{0}. This proves Property (c).

It remains to prove (d). Since E is a Poisson integrator, the kernel of E^⊤−Id includes the kernel of J. As a consequence, the range of E−Id is included in the range of J. This proves Property (d).

The proof is now complete.               ⃞

3.1.2. The evolution operator

Our first result states that Hamiltonian (Poisson) systems and symplectic (Poisson) matrix integrators are deeply related in an essential manner.⁹

Proposition 3.2 (The evolution matrix is a Poisson integrator matrix which commutes with A). Consider the Poisson system $\frac{d}{dt}\text{u}=\text{A}\text{u}$, where A = JH. Then,

$\begin{array}{ll}\text{(i)}& \text{The evolution matrix}\text{E}(t):=e^{\text{A}t}\text{is a Poisson integrator matrix,}\\ \text{(ii)}& \left[\text{E}(t),\text{A}\right]:=E(t)\text{A}-\text{A}\text{E}(t)=0.\end{array}$

Proof. Let us prove that E(t) is a Poisson matrix. We have that

$\begin{array}{l}\frac{d}{dt}(\text{E}(t)\text{ J }{\text{E}}^{\top }(t))=\frac{d}{dt}(\text{E}(t))\text{ J }{\text{E}}^{\top }(t)+\text{E}(t)\text{ J }\frac{d}{dt}({\text{E}}^{\top }(t))\\ =\text{E}(t)(\text{AJ}+{\text{JA}}^{\top }){\text{E}}^{\top }(t)\\ =\text{E}(t)(\text{JHJ}+{\text{JHJ}}^{\top }){\text{E}}^{\top }(t)\\ =\text{0},\end{array}$

because J^⊤ = −J. As a consequence, we get that E(t)JE^⊤(t) = E(0)JE^⊤(0) = J. This proves that E(t)JE^⊤(t) = J for all t∈ℝ.

Let us now complete the proof of Property (i) by showing that E^⊤(t)v = v for all v in the kernel of J. Observe that

$\frac{d}{dt}{\text{E}}^{\top }(t)\text{v}={\text{A}}^{\top }{\text{E}}^{\top }\text{v}={\text{E}}^{\top }{\text{A}}^{\top }\text{v}=0.$

Thus, E^⊤(t)v = E^⊤(0)v = v.

Property (ii) clearly holds. This completes the proof.              ⃞

3.1.3. Conserved quantities

Next, we characterize all possible linear and quadratic functionals which remain constant on the orbits of any linear system. It does not have to be a Hamiltonian system.

Proposition 3.3 (Characterization of linear and quadratic first integrals). Let t↦u(t) be any of the orbits of the linear system $\frac{d}{dt}\text{u}=\text{A}\text{u}$. Then,

(i) The mapping ↦ v^Tu(t) is constant if and only if A^Tv = 0,

(ii) If the matrix S is symmetric, the mapping t ↦ u^T(t)S u(t) is constant if and only if S A + A^TS = 0.

For Poisson systems, a couple of simple consequences can be readily obtained. The first is that, in (i), the linear functional t↦v^⊤u(t) is constant in time if and only if A^⊤v = −HJv = 0. In other words, the system has linear first integrals if only and only if either J or H is not invertible. If they are, then so is A^⊤ and the only linear first integral of the system is the trivial one.

The second consequence is that, in (ii), we can take S: = H. Indeed, SA+A^⊤S = HJH+HJ^⊤H = 0, because J is anti-symmetric. This shows that, for any Hamiltonian system, the Hamiltonian ${H}(\text{u}):=\frac{1}{2}{\text{u}}^{\top }\text{H}\text{u}$ is always a quadratic first integral.

Proof. Let us prove Property (i). We have

$\frac{d}{dt}({\text{v}}^{\top }\text{u})={\text{v}}^{\top }\frac{d}{dt}\text{u}={\text{v}}^{\top }\text{A}\text{u}.$

This implies Property(i).

Let us now prove Property (ii). We have

$\begin{array}{l}\frac{d}{dt}{(\text{u}}^{\top }\text{S}\text{u})=(\frac{d}{dt}{\text{u}}^{\top })\text{S}\text{u}+{\text{u}}^{\top }\text{S}(\frac{d}{dt}\text{u})={\text{u}}^{\top }{(\text{A}}^{\top }\text{S}+\text{S}\text{A})\text{u},\end{array}$

and we see that, for the last quantity to be equal to zero, the matrix A^⊤S+SA has to be anti-symmetric. Since this matrix is also symmetric, because S is symmetric, it must be identically zero. This implies Property (ii) and completes the proof.                   ⃞

3.2. Approximate Poisson evolution operators commuting with A

From now on, we focus on finding out when it is possible to guarantee that the linear and quadratic first integrals of the original Hamiltonian system are also maintained constant after discretization. We give a characterization of the Poisson Runge-Kutta methods and show that these methods do conserve the linear and quadratic first integrals of the original Hamiltonian. This exact conservation property is the first main result of this section. It justifies the relevance of Poisson numerical methods for the time-discretizations of Hamiltonian ODEs.

So, consider now the approximation u_n to u(Δtn) given by the one-step numerical scheme

${\text{u}}_{n+1}:={\text{E}}_{\Delta t}{\text{u}}_n.$

If E_Δt is a symplectic matrix, we say that the scheme is symplectic. If E_Δt is a Poisson matrix, we say that the scheme is a Poisson scheme.

Clearly, the discrete evolution matrix E_Δt is an approximation to the exact evolution matrix E(Δt) = e^AΔt. Since, by Proposition 0.0.2, the exact evolution matrix is a Poisson integrator matrix, one wonders what do we gain by requiring that the discrete evolution matrix E_Δt to be also a Poisson integrator matrix.

It is possible to answer that question in a very satisfactory manner if we additionally require that the discrete evolution matrix E_Δt commute with A. As we see next, if this is the case, the original Hamiltonian remains constant on the discrete orbits of the numerical scheme.

Proposition 3.4. Assume that

$\begin{array}{ll}\text{(i)}& \text{The discrete evolution matrix}{\text{E}}_{\Delta t}\text{is a Poisson integrator matrix},\\ \text{(ii)}& \left[{\text{E}}_{\Delta t},\text{A}\right]:={\text{E}}_{\Delta t}\text{A}-\text{A}{\text{E}}_{\Delta t}=0.\end{array}$

Then the mapping $n\mapsto \frac{1}{2}{\text{u}}_n^{\top }\text{H}{\text{u}}_n$ is constant.

Proof. Since

$\begin{array}{ll}\frac{1}{2}{\text{u}}_{n+1}^{\top }\text{H}{\text{u}}_{n+1}=& \frac{1}{2}{\text{u}}_n^{\top }{\text{E}}_{\Delta t}^{\top }\text{H}{\text{E}}_{\Delta t}{\text{u}}_n=\frac{1}{2}{\text{u}}_n^{\top }\text{H}{\text{u}}_n+\frac{1}{2}\Theta ,\end{array}$

Where $\Theta :={\text{u}}_n^{\top }({\text{E}}_{\Delta t}^{\top }\text{H}{\text{E}}_{\Delta t}-\text{H}){\text{u}}_n$, we only have to show that Θ = 0.

We know that Jℝ^d⊕J⁻¹{0} = ℝ^d so we can write u_n = Jw+v, for some w∈ℝ^d and v∈J⁻¹{0}. Then, we insert this expression into the definition of Θ to obtain

$\Theta ={\Theta }_{w,w}+2{\Theta }_{v,w}+{\Theta }_{v,v},$

where

$\begin{array}{l}{\Theta }_{w,w}={\text{w}}^{\top }{\text{J}}^{\top }({\text{E}}_{\Delta t}^{\top }{\text{HE}}_{\Delta t}-\text{H})\text{ J}\text{w},\\ {\Theta }_{v,w}={\text{v}}^{\top }({\text{E}}_{\Delta t}^{\top }{\text{HE}}_{\Delta t}-\text{H})\text{ J}\text{w},\\ {\Theta }_{v,v}={\text{v}}^{\top }({\text{E}}_{\Delta t}^{\top }{\text{HE}}_{\Delta t}-\text{H})\text{v}\\ ={\text{v}}^{\top }({\text{E}}_{\Delta t}^{\top }\text{H}-{\text{HE}}_{\Delta t}^{-1}){\text{E}}_{\Delta t}\text{v}.\\ ={\text{v}}^{\top }({\text{E}}_{\Delta t}^{\top }\text{H}-{\text{HE}}_{\Delta t}^{-1})(\text{v}+\text{JB}\text{v}),\end{array}$

by Property (d) of Proposition 3.1. We can rewrite these quantities in a more convenient manner with the help of the following matrix:

$\Phi :={\text{E}}_{\Delta t}^{\top }\text{H}-\text{H}{\text{E}}_{\Delta t}^{-1}.$

Indeed, we can now write that

$\begin{array}{l}{\Theta }_{w,w}={\text{w}}^{\top }{\text{J}}^{\top }\Phi {\text{E}}_{\Delta t}\text{ J}\text{w},{\Theta }_{v,w}={\text{v}}^{\top }\Phi {\text{E}}_{\Delta t}\text{ J}\text{w}\text{ and}\\ {\Theta }_{v,v}={\text{v}}^{\top }\Phi (\text{v}+\text{JB}\text{v}).\end{array}$

But

$\begin{array}{l}\Phi {\text{E}}_{\Delta t}\text{J}=({\text{E}}_{\Delta t}^{\top }\text{H}-\text{H}{\text{E}}_{\Delta t}^{-1}){\text{E}}_{\Delta t}\text{J}\\ ={\text{E}}_{\Delta t}^{\top }\text{H}{\text{E}}_{\Delta t}\text{J}+{\text{A}}^{\top }\\ =({\text{E}}_{\Delta t}^{\top }\text{H}{\text{E}}_{\Delta t}\text{J}{\text{E}}_{\Delta t}^{\top }+{\text{A}}^{\top }{\text{E}}_{\Delta t}^{\top }){\text{E}}_{\Delta t}^{-\top }\\ =({\text{E}}_{\Delta t}^{\top }\text{H}({\text{E}}_{\Delta t}\text{J}{\text{E}}_{\Delta t}^{\top }-\text{J})+{\text{E}}_{\Delta t}^{\top }\text{H}\text{J}+{\text{A}}^{\top }{\text{E}}_{\Delta t}^{\top }){\text{E}}_{\Delta t}^{-\top }\\ =({\text{E}}_{\Delta t}^{\top }\text{H}({\text{E}}_{\Delta t}\text{J}{\text{E}}_{\Delta t}^{\top }-\text{J})+{\left[{\text{E}}_{\Delta t},\text{A}\right]}^{\top }){\text{E}}_{\Delta t}^{-\top }.\end{array}$

Since E_Δt is a Poisson matrix and commutes with A, we have that Φ = 0 on the range of E_ΔtJ, and, by Property (a) of Proposition 3.1, on the range of J. As a consequence, we get that Θ = v^⊤Φv.

It remains to prove that v^⊤Φv equals to zero. Since Jℝ^d⊕J⁻¹{0} = ℝ^d, we can write that Hv = y+h where y is in the range of J and h in kernel of J. Since we can consider that J⁻¹ is defined on the range of J, we can write that y = J⁻¹x for some element x of the range of J. Then

$\begin{array}{l}{\text{v}}^{\top }\Phi \text{v}={\text{v}}^{\top }({\text{E}}_{\Delta t}^{\top }\text{H}-\text{H}{\text{E}}_{\Delta t}^{-1})\text{v}\hfill \\ ={\text{v}}^{\top }({\text{E}}_{\Delta t}^{\top }\text{H}\text{v})-{({\text{E}}_{\Delta t}^{\top }\text{H}\text{v})}^{\top }\text{v}\hfill \\ ={\text{v}}^{\top }({\text{E}}_{\Delta t}^{\top }(J^{-1}\text{x}+\text{h}))-{({\text{E}}_{\Delta t}^{-\top }(J^{-1}\text{x}+\text{h}))}^{\top }\text{v}\hfill \\ ={\text{v}}^{\top }({\text{E}}_{\Delta t}^{\top }{J}^{-1}\text{x})-{({\text{E}}_{\Delta t}^{-\top }{J}^{-1}\text{x})}^{\top }\text{v}\hfill \\ \text{ \hspace{0.17em}since }{\text{E}}^{\top }=\text{Id on }{J}^{-1}\left\{0\right\},\hfill \\ ={\text{v}}^{\top }(J^{-1}{\text{E}}_{\Delta t}^{-1}\text{x})-{(J^{-1}{\text{E}}_{\Delta t})}^{\top }\text{x}\text{v}\hfill \\ \text{ by (b) of Proposition }3.1,\hfill \\ =0,\hfill \end{array}$

because of the orthogonality of the kernel of J with its range. This completes the proof.               ⃞

3.2.1. Evolution matrices which are rational functions of hA

Next, we consider an important example of evolution matrix E_Δt satisfying the assumptions of the previous result, that is, E_Δt is a Poisson integrator matrix and commutes with A.

Proposition 3.5. Set A: = JH where J is antisymmetric and H is symmetric. Set also E_Δt: = R(ΔtA), where z↦R(z) is a rational function. Then, E_Δt is a Poisson integrator matrix if and only if R(z)R(−z) = 1 and R(0) = 1.

This result has a simple and remarkable geometric interpretation. It states that the Poisson numerical methods whose evolution matrix is a rational function of hA are those for which we can recover the point we were evolving from by just applying the same scheme with time reversed, that is, the schemes for which,

$\text{if}{\text{u}}_{n+1}:={\text{E}}_{\Delta t}{\text{u}}_n\text{then }{\text{u}}_n={\text{E}}_{-\Delta t}{\text{u}}_{n+1}.$

Let us now prove the result.

Proof. Set Z: = ΔtA. By definition, the discrete evolution operator E_Δt: = R(ΔtA) is a Poisson integrator matrix if

(i) R(Z)JR(Z)^⊤ = J,

(ii) R(Z)^⊤ = Id on the kernel of J.

We want to show that these two properties are equivalent to

(1) R(z)R(−z) = 1,

(2) R(0) = 1.

Since R is a rational function, there are polynomials P and Q such that R(z) = (Q(z))⁻¹P(z). In terms of these polynomials, the above conditions read

(a) P(z)P(−z) = Q(z)Q(−z),

(b) P(0) = Q(0).

We only prove that (a) and (b) imply (i) and (ii), as the converse is equally easy.

Let v be an arbitrary element of the kernel of J. So (ΔtHJ)^ℓv = 0, for ℓ>0, and since P is a polynomial matrix we have

$P{(\text{Z})}^{\top }\text{v}=P(-\Delta t\text{H}\text{J})\text{v}=P(0)\text{v}.$

By (b), we get

$P(0)\text{v}=Q(0)\text{v}=Q(-\Delta t\text{H}\text{J})\text{v}=Q{(\text{Z})}^{\top }\text{v},$

and then Property (ii) follows.

Let us prove Property (i). Take u: = Jw, where w is an arbitrary element of ℝ^d. By (a),

$\begin{array}{l}0=Q(\text{Z})Q(-\text{Z})\text{u}-P(\text{Z})P(-\text{Z})\text{u}=\\ Q(\text{Z})Q(-\text{Z})\text{J}\text{w}-P(\text{Z})P(-\text{Z})\text{J}\text{w},\end{array}$

and since

${(-\text{A})}^i\text{J}={(-\text{A})}^{i-1}(-\text{J}\text{H})\text{J}={(-\text{A})}^{i-1}\text{J}({\text{A}}^{\top })=\text{J}{({\text{A}}^{\top })}^i,$

because J is anti-symmetric, we get that

$\begin{array}{l}0=Q(\text{Z})Q(-\text{Z})\text{J}\text{w}-P(\text{Z})P(-\text{Z})\text{J}\text{w}\\ =Q(\text{Z})\text{J}Q({\text{Z}}^{\top })\text{w}-P(\text{Z})\text{J}P({\text{Z}}^{\top })\text{w}.\end{array}$

As this holds for all elements w of ℝ^d, we obtain that

$Q(\text{Z})\text{J}Q({\text{Z}}^{\top })=P(\text{Z})\text{J}P({\text{Z}}^{\top }),$

and Property (i) follows. This completes the proof.              ⃞

3.2.2. Symplectic Runge-Kutta methods

Perhaps the main example of numerical methods whose discrete evolution operator E_Δt is a rational function of ΔtA is the Runge-Kutta methods. In fact, if the Butcher array of the Runge-Kutta method is

then

$R(z):=1+z{b}^{\top }{(Id-zA)}^{-1}11=Q{(z)}^{-1}P(z),$

Where 11 is the s-dimensional column vector whose components are one, and, see [31],

$P(z)=\text{det}(Id-zA\text{ + }z11b^{\top })\text{ and }Q(z)=\text{det}(Id-zA).$

This means that the previous result characterizes all the symplectic Runge-Kutta methods. Two important conclusions can be immediately drawn:

(i) No explicit Runge-Kutta method (that is, Q is not a constant) is symplectic.

(ii) Any Runge-Kutta method for which Q(z) = P(−z) is symplectic. Two notable examples are:

(α) The s-stage Gauss-Legendre Runge-Kutta method, which is of order 2s. The function R(z) is the s-th diagonal Padé approximation to e^z.

(β) Any diagonally implicit Runge-Kutta method for which

$A=\left[\begin{array}{ccccc}\frac{1}{2}{b}_1& 0& 0& ...& 0\\ b_1& \frac{1}{2}{b}_2& 0& ...& 0\\ b_1& b_2& \frac{1}{2}{b}_3& ...& 0\\ ...& ...& ...& ...& 0\\ b_1& b_2& b_3& ...& \frac{1}{2}{b}_s\end{array}\right].$

The function R(z) is given by ${\Pi }_{i=1}^s\frac{1+z{b}_i/2}{1-z{b}_i/2}$.

3.2.3. Conserved quantities

The Poisson integrator Runge-Kutta methods also maintain constant all the linear and quadratic first integrals of the original Poisson system. As pointed out above, this result unveils the importance of these methods for reaching our objective, namely, the conservation of the first integrals of the original Poisson system.

Theorem 3.1 (Exact conservation of linear and quadratic first integrals). All linear and quadratic first integrals of the original Poisson system $\frac{d}{dt}\text{u}=\text{A}\text{u}$ remain constant on the orbits n↦u_n provided by a Poisson integrator Runge-Kutta method. That is,

$\begin{array}{ll}\text{(i)}& \text{The mapping}n\mapsto {\text{v}}^{\top }{\text{u}}_n\text{is constant if }{\text{A}}^{\top }\text{v}=0,\\ \text{(ii)}& \text{The mapping}n\mapsto {\text{u}}_n^{\top }\text{S}{\text{u}}_n\text{is constant if }\text{S}\text{A}+{\text{A}}^{\top }\text{S}=0\text{.}\end{array}$

Proof. Let us prove Property (i). We have

$\begin{array}{l}{\text{v}}^{\top }{\text{u}}_{n+1}={\text{v}}^{\top }R(\Delta t\text{A}){\text{u}}_n\\ ={\text{v}}^{\top }R(0){\text{u}}_n\text{because}{\text{ A}}^{\top }\text{v}=0,\\ ={\text{v}}^{\top }{\text{u}}_n,\end{array}$

because the evolution matrix for Δt = 0 must be the identity, that is, R(0) = Id. This proves Property (i). Note we did not use that the numerical scheme is Poisson, we only used the form of the evolution operator.

Now, let us prove Property (ii). We have

$\begin{array}{l}{\text{u}}_{n+1}^{\top }\text{S}{\text{u}}_{n+1}={\text{u}}_n^{\top }R{(\Delta t\text{A})}^{\top }\text{S}R(\Delta t\text{A}){\text{u}}_n\end{array}$

Since SA = −A^⊤S, we have that SP(ΔtA) = P(−ΔtA^⊤)S, and, as a consequence, that

$Q{(-\Delta t{\text{A}}^{\top })}^{-1}\text{S}=\text{S}Q{(\Delta t\text{A})}^{-1}.$

We can then write that

$\begin{array}{l}{\text{u}}_{n+1}^{\top }\text{S}{\text{u}}_{n+1}={\text{u}}_n^{\top }R{(\Delta t\text{A})}^{\top }P(-\Delta t{\text{A}}^{\top })\text{S}Q{(\Delta t\text{A})}^{-1}{\text{u}}_n\\ ={\text{u}}_n^{\top }R{(\Delta t\text{A})}^{\top }P(-\Delta t{\text{A}}^{\top })Q{(-\Delta t{\text{A}}^{\top })}^{-1}\text{S}{\text{u}}_n\\ ={\text{u}}_n^{\top }R{(\Delta t\text{A})}^{\top }R(-\Delta t{\text{A}}^{\top })\text{S}{\text{u}}_n\\ ={\text{u}}_n^{\top }\text{S}{\text{u}}_n,\end{array}$

since the scheme is a Poisson integrator. This proves Property (ii) and completes the proof.              ⃞

3.3. Poisson integrator matrices which do not commute with A

Now, we turn our attention to discrete Poisson methods whose evolution matrix E_Δt does not commute with the matrix A defining the original system. The loss of this commutativity property has dire consequences, as the exact conservation of the original Hamiltonian cannot be guaranteed anymore. Indeed, as we saw in the proof of Proposition 3.4, when the structure matrix J is invertible and E_Δt is symplectic, we have that

$\begin{array}{ll}\frac{1}{2}{\text{u}}_{n+1}^{\top }\text{H}{\text{u}}_{n+1}=& \frac{1}{2}{\text{u}}_n^{\top }\text{H}{\text{u}}_n-\frac{1}{2}{\text{u}}_{n+1}{\text{J}}^{-1}\left[{\text{E}}_{\Delta t},\text{A}\right]{\text{u}}_n,\end{array}$

Which means that the original Hamiltonian is not maintained constant, as claimed.

Nevertheless, because the matrix evolution E_Δt is a Poisson integrator matrix, we can still get

(i) the existence of a matrix A_Δt such that AJ is symmetric and that ${\text{E}}_{\Delta t}=e^{\Delta t{\text{A}}_{\Delta t}}$,

(ii) the non-drifting property of quadratic first integrals (on bounded orbits) of the original Hamiltonian.

The property (i) states that the discrete dynamics of the numerical scheme is exactly that of a dynamical system (Hamiltonian if E_Δt is symplectic) defined by a matrix A_Δt which we expect to converge to A as the discretization parameter Δt tends to zero. This property will allow us to obtain the non-drifting property (ii).

The non-drifting property (ii) states that, provided that the orbits are bounded, the original quadratic first integrals oscillate around the correct value for all times and are away from it a number proportional to the size of the truncation error of the numerical scheme. This remarkable, non-drifting property is the second main result of this section. This is why we are interested in working with Poisson integrators of this type.

3.3.1. The approximate Poisson system

We begin by showing that there is a Poisson system whose continuous orbits contain the discrete orbits generated by the evolution matrix E_Δt.

Proposition 3.6. Let the mapping Δt → E_Δt, where the discrete evolution matrix E_Δt is a Poisson matrix, be continuous on (0, T) and such that E_{Δt = 0} = Id. Then, for Δt small enough,

(i) The matrix ${\text{A}}_{\Delta t}:=\frac{1}{\Delta t}log{\text{E}}_{\Delta t}$ is such that A_ΔtJ is symmetric.

(ii) The orbits n↦u_n generated by the discrete evolution operator E_Δt lie on the orbits t↦u(t) of the Poisson system $\frac{d}{dt}\text{u}={\text{A}}_{\Delta t}\text{u}$.

Proof. Let us prove Property (i). We have to show that A_ΔtJ is symmetric. But, for Δt small enough, we can write

$\begin{array}{l}{\text{A}}_{\Delta t}\text{J}=\frac{1}{\Delta t}log{\text{E}}_{\Delta t}\text{J}=\frac{1}{\Delta t}{\displaystyle \sum _{\ell =1}^{\infty }}\frac{{(-1)}^{\ell +1}}{\ell }{({\text{E}}_{\Delta t}-\text{Id})}^{\ell }\text{J}.\end{array}$

Using the fact that E_Δt is Poisson, we get that

$({\text{E}}_{\Delta t}-\text{Id})\text{J}=({\text{E}}_{\Delta t}\text{J}-\text{J})=(\text{J}{\text{E}}_{\Delta t}^{-\top }-\text{J})=\text{J}({\text{E}}_{\Delta t}^{-\top }-\text{Id}),$

and so, we can write that

$\begin{array}{l}{\text{A}}_{\Delta t}\text{J}=\frac{1}{\Delta t}{\displaystyle \sum _{\ell =1}^{\infty }}\frac{{(-1)}^{\ell +1}}{\ell }\text{J}{({\text{E}}_{\Delta t}^{-\top }-\text{Id})}^{\ell }\\ =\text{J}\frac{1}{\Delta t}log{\text{E}}_{\Delta t}^{-\top }=\text{J}\frac{1}{\Delta t}log{\text{E}}_{\Delta t}^{-\top }={\text{J}}^{\top }{\text{A}}_{\Delta t}^{\top }.\end{array}$

This proves Property (i).

It remains to prove Property (ii). But the discrete orbits are made of the points

${\text{u}}_n={({\text{E}}_{\Delta t})}^n{\text{u}}_0={(e^{\Delta t{\text{A}}_{\Delta t}})}^n{\text{u}}_0=e^{n\Delta t{\text{A}}_{\Delta t}}{\text{u}}_0,$

Which lie on the orbit $t\mapsto e^{t{\text{A}}_{\Delta t}}{\text{u}}_0$. This proves Property (ii) and completes the proof.                      ⃞

3.3.2. Closeness of the approximate system

Although Proposition 3.3, with A replaced by A_Δt, gives information about the linear and quadratic first integrals of the new system $\frac{d}{dt}\text{u}={\text{A}}_{\Delta t}\text{u}$, we are more interested in knowing how close are those invariants to those of the original system. We have then to relate the matrices A and A_Δt in order to estimate how well the numerical scheme approximates the original quadratic first integrals.

We show that those matrices are as close as the order of accuracy of the numerical scheme. We say that the order of accuracy of the numerical scheme u_n+1 = E_Δtu_n is p if

$||\text{E}(\Delta t)-{\text{E}}_{\Delta t}||\le C{(\Delta t)}^{p+1},$

for some constant C independent of Δt.

Proposition 3.7. Assume that Δt ≤ ln(3/2)/max{||A||, ||A_Δt||}. Then

$||\text{A}-{\text{A}}_{\Delta t}||\le \frac{2}{\Delta t}||\text{E}(\Delta t)-{\text{E}}_{\Delta t}||.$

Proof. We have that

$\begin{array}{l}\text{A}-{\text{A}}_{\Delta t}=\frac{1}{\Delta t}(\log \text{E}(\Delta t)-\log {\text{E}}_{\Delta t})\hfill \\ =\frac{1}{\Delta t}{\displaystyle \sum _{\ell =1}^{\infty }}\frac{{(-1)}^{\ell +1}}{\ell }({(\text{E}(\Delta t)-Id)}^{\ell }-{({\text{E}}_{\Delta t}-Id)}^{\ell })\hfill \\ =\frac{1}{\Delta t}{\displaystyle \sum _{\ell =1}^{\infty }}\frac{{(-1)}^{\ell +1}}{\ell }{\displaystyle \sum _{m=1}^{\ell }}{({\text{E}}_{\Delta t}-Id)}^{m-1}\hfill \\ (\text{E}(\Delta t)-{\text{E}}_{\Delta t}){(\text{E}(\Delta t)-Id)}^{\ell -m},\hfill \end{array}$

because ${\text{R}}^{\ell }-{\text{S}}^{\ell }=\sum _{m=1}^{\ell }{\text{S}}^{m-1}(\text{R}-\text{S}){\text{R}}^{\ell -m}.$ This implies that

$\begin{array}{l}\Vert \text{A}-{\text{A}}_{\Delta t}\Vert \le \frac{1}{\Delta t}{\displaystyle \sum _{\ell =1}^{\infty }\frac{1}{\ell }}[{\displaystyle \sum _{m=1}^{\ell }\Vert {\text{E}}_{\Delta t}-\text{Id}{\Vert }^{m-1}}\Vert \text{E}(\Delta t)-\text{Id}{\Vert }^{\ell -m}]\Vert \\ \text{ E}(\Delta t)-{\text{E}}_{\Delta t}\Vert \\ \le \frac{1}{\Delta t}{\displaystyle \sum _{\ell =1}^{\infty }{\max }_{1\le m\le \ell }}\{\Vert {\text{E}}_{\Delta t}-\text{Id}{\Vert }^{m-1}\Vert \text{E}(\Delta t)-\text{Id}{\Vert }^{\ell -m}\}\\ \Vert \text{E}(\Delta t)-{\text{E}}_{\Delta t}\Vert \end{array}$

Since ||e^tB−Id|| ≤ e^t||B||−1, we get that

$\begin{array}{l}\Vert \text{A}-{\text{A}}_{\Delta t}\Vert \le \frac{1}{\Delta t}{\displaystyle \sum _{\ell =1}^{\infty }}\underset{1\le m\le \ell }{\max }\{{(e^{h\Vert {\text{A}}_{\Delta t}\Vert }-1)}^{m-1}{(e^{h\Vert \text{A}\Vert }-1)}^{\ell -m}\}\hfill \\ \Vert \text{E}(\Delta t)-{\text{E}}_{\Delta t}\Vert \hfill \\ \le \frac{1}{\Delta t}{\displaystyle \sum _{\ell =1}^{\infty }}{(1/2)}^{\ell -1}\Vert \text{E}(\Delta t)-{\text{E}}_{\Delta t}\Vert \hfill \\ =\frac{2}{\Delta t}\Vert \text{E}(\Delta t)-{\text{E}}_{\Delta t}\Vert ,\hfill \end{array}$

because $e^{\Delta tmax\left\{\left|\right|\text{A}\left|\right|,\left|\right|{\text{A}}_{\Delta t}\left|\right|\right\}}\le 3/2$. This completes the proof.              ⃞

3.3.3. Non-drifting of quadratic first integrals

We end by noting that the quadratic first integrals of the original Hamiltonian system do not drift on the discrete orbits generated by the numerical scheme.

Theorem 3.2 (Non-drifting property of the quadratic first integrals). Let the discrete evolution matrix E_Δt be a Poisson matrix, and let n↦u_n denote any of its orbits. Assume that ${\lambda }_{\min }:=\underset{\text{v}\ne 0}{\inf }\left|\frac{{\text{v}}^{\top }{\text{S}}_{\Delta t}\text{v}}{{\text{v}}^{\top }\text{v}}\right|>0$. Then

$|{\text{u}}_n^{\top }\text{S}{\text{u}}_n-{\text{u}}_n^{\top }{\text{S}}_{\Delta t}{\text{u}}_n|\le \frac{1}{{\lambda }_{\min }}\Vert \text{S}-{\text{S}}_{\Delta t}\Vert \left|{\text{u}}_n^{\top }{\text{S}}_{\Delta t}{\text{u}}_n\right|\forall n.$

Proof. Since ${\text{u}}_n^{\top }\text{S}{\text{u}}_n-{\text{u}}_n^{\top }{\text{S}}_{\Delta t}{\text{u}}_n={\text{u}}_n^{\top }(\text{S}-{\text{S}}_{\Delta t}){\text{u}}_n,$ we have that

$\begin{array}{l}|{\text{u}}_n^{\top }\text{S}{\text{u}}_n-{\text{u}}_n^{\top }{\text{S}}_{\Delta t}{\text{u}}_n|\le \underset{\text{v}\ne 0}{\sup }\left|\frac{{\text{v}}^{\top }(\text{S}-{\text{S}}_{\Delta t})\text{v}}{{\text{v}}^{\top }{\text{S}}_{\Delta t}\text{v}}\right|\left|{\text{u}}_n^{\top }{\text{S}}_{\Delta t}{\text{u}}_n\right|\\ \le \underset{\text{v}\ne 0}{\sup }\left|\frac{{\text{v}}^{\top }\text{v}}{{\text{v}}^{\top }{\text{S}}_{\Delta t}\text{v}}\right|\Vert \text{S}-{\text{S}}_{\Delta t}\Vert \left|{\text{u}}_n^{\top }{\text{S}}_{\Delta t}{\text{u}}_n\right|\\ \le \frac{1}{{\lambda }_{\min }}\Vert \text{S}-{\text{S}}_{\Delta t}\Vert \left|{\text{u}}_n^{\top }{\text{S}}_{\Delta t}{\text{u}}_n\right|.\end{array}$

This completes the proof.                ⃞

In the case in which S = H = J⁻¹A, we can take ${\text{S}}_{\Delta t}={\text{J}}^{-1}{\text{A}}_{\Delta t}$ to get that

$\begin{array}{l}|{\text{u}}_n^{\top }\text{S}{\text{u}}_n-{\text{u}}_n^{\top }{\text{S}}_{\Delta t}{\text{u}}_n|\le \frac{\Vert {\text{J}}^{-1}\Vert }{{\lambda }_{\min }}\Vert \text{A}-{\text{A}}_{\Delta t}\Vert \left|{\text{u}}_n^{\top }{\text{S}}_{\Delta t}{\text{u}}_n\right|\\ \le \frac{2\Vert {\text{J}}^{-1}\Vert }{\Delta t{\lambda }_{\min }}\Vert \text{E}-{\text{E}}_{\Delta t}\Vert \left|{\text{u}}_n^{\top }{\text{S}}_{\Delta t}{\text{u}}_n\right|,\end{array}$

and, if the scheme is of order p, we immediately get that

$|{\text{u}}_n^{\top }\text{S}{\text{u}}_n-{\text{u}}_n^{\top }{\text{S}}_{\Delta t}{\text{u}}_n|\le C\frac{2\Vert {\text{J}}^{-1}\Vert }{{\lambda }_{\min }}\left|{\text{u}}_n^{\top }{\text{S}}_{\Delta t}{\text{u}}_n\right|\Delta t^p,$

for Δt small enough. Since $n\mapsto {\text{u}}_n^{\top }{\text{S}}_{\Delta t}{\text{u}}_n$ is constant, we see that ${\text{u}}_n^{\top }\text{S}{\text{u}}_n$ does not drift in time, as claimed. Of course, for this to take place, λ_min has to be strictly positive, which happens when ${\text{v}}^{\top }{\text{S}}_{\Delta t}\text{v}$ is never equal to zero on the unit sphere v^⊤v = 1.

3.4. Separable Hamiltonians and partitioned Runge-Kutta methods

Since there are no explicit symplectic Runge-Kutta methods for general Hamiltonians, one wonders if the situation changes for some subclass of Hamiltonians. It turns out that this is the case for the so-called separable Hamiltonians. Indeed, there are explicit, partitioned Runge-Kutta which are symplectic when applied to separable Hamiltonian systems. This is the third main result of this short introduction.

3.4.1. Separable Hamiltonians

If the Hamiltonian is separable, that is, if

${H}=\frac{1}{2}{(\text{p}}^{\top }{\text{H}}_{pp}\text{p}+{\text{q}}^{\top }{\text{H}}_{qq}\text{q})$

and if the corresponding Hamiltonian system is of the form

$\frac{d}{dt}{\displaystyle \underset{\text{u}}{{\displaystyle \underbrace{\left[\begin{array}{c}\text{p}\\ \text{q}\end{array}\right]}}}}={\displaystyle \underset{\text{J}}{{\displaystyle \underbrace{\left[\begin{array}{cc}0& {\text{J}}_{pq}\\ {\text{J}}_{qp}& 0\end{array}\right]}}}}{\displaystyle \underset{\text{H}}{{\displaystyle \underbrace{\left[\begin{array}{cc}{\text{H}}_{pp}& 0\\ 0& {\text{H}}_{qq}\end{array}\right]}}}}\left[\begin{array}{c}\text{p}\\ \text{q}\end{array}\right]={\displaystyle \underset{\text{A}}{{\displaystyle \underbrace{\left[\begin{array}{cc}0& {\text{A}}_{pq}\\ {\text{A}}_{qp}& 0\end{array}\right]}}}}\left[\begin{array}{c}\text{p}\\ \text{q}\end{array}\right],$

it is interesting to explore schemes which treat the q- and p-components of the unknown u in different ways.

3.4.2. Partitioned Runge-Kutta methods

The so-called partitioned Runge-Kutta methods are of the form

$\begin{array}{l}{\text{E}}_{\Delta t}\text{u}=\text{u}+\Delta t{\displaystyle \sum _{j=1}^s}\left[\begin{array}{cc}0& b_j{\text{A}}_{pq}\\ {\overline{b}}_j{\text{A}}_{qp}& 0\end{array}\right]{\text{U}}_j\text{and}\\ {\text{U}}_i=\text{u}+\Delta t{\displaystyle \sum _{j=1}^s}\left[\begin{array}{cc}0& a_{ij}{\text{A}}_{pq}\\ {\overline{a}}_{ij}{\text{A}}_{qp}& 0\end{array}\right]{\text{U}}_j,\end{array}$

and can be thought as associated to two Butcher arrays,

To emphasize its relation with the standard RK methods, we can rewrite them as

$\begin{array}{l}{\text{E}}_{\Delta t}\text{u}=\text{u}+\Delta t\text{A}{\displaystyle \sum _{j=1}^s}{B}_j{\text{U}}_j\text{where}\\ {\text{U}}_i=\text{u}+\Delta t\text{A}{\displaystyle \sum _{j=1}^s}{A}_{ij}{\text{U}}_{j}i=1,\dots ,s,\end{array}$

where $B_{\ell }:=\left[\begin{array}{cc}{\overline{b}}_{\ell }\text{I}{\text{d}}_{pp}& 0\\ 0& b_{\ell }\text{I}{\text{d}}_{qq}\end{array}\right]$ and $A_{\ell m}:=\left[\begin{array}{cc}{\overline{a}}_{\ell m}\text{I}{\text{d}}_{pp}& 0\\ 0& a_{\ell m}\text{I}{\text{d}}_{qq}\end{array}\right]$. We can see that we recover the standard RK methods whenever ${\overline{b}}_{\ell }=b_{\ell }$ and ${\overline{a}}_{ij}=a_{ij}$.

3.4.3. Conservation of the original linear first integrals

It is clear that these methods have an evolution matrix which does not commute with A. Because of this, all the results in the previous subsection apply to them when they are symplectic. In addition, they maintain constant all the linear first integrals of the original Hamiltonian, not because of their symplecticity but because of the structure of their evolution matrix.

Theorem 3.3 (Exact conservation of linear first integrals). Suppose that the original Hamiltonian is separable. Then, all linear first integrals of the original Hamiltonian are also first integrals of any partitioned Runge-Kutta method.

Proof. We have that ${\text{v}}^{\top }{\text{E}}_{\Delta t}\text{u}={\text{v}}^{\top }(\text{u}+\Delta t\text{A}\sum _{j=1}^s{B}_j{\text{U}}_j)={\text{v}}^{\top }\text{u},$ whenever A^⊤v = 0. This completes the proof.                     ⃞

3.4.4. Symplecticity

We end this section by providing a characterization of the partitioned Runge-Kutta methods which are Poisson (or symplectic), see a detailed proof in Appendix A, and by showing that they can be explicit.

Proposition 3.8. Suppose that the original Hamiltonian is separable. A partitioned Runge-Kutta method is Poisson (or symplectic) if an only if

$b_i{\overline{a}}_{ij}+{\overline{b}}_j{a}_{ji}-b_i{\overline{b}}_j=0i,j=1,\dots ,s.$

When the partitioned Runge-Kutta methods becomes a standard Runge-Kutta method, this result is equivalent to Proposition 3.5, as it can be seen after performing a simple computation.

3.4.5. Explicit methods

An important example symplectic partitioned Runge-Kutta methods have theA-matrices of the Butcher arrays of the following form:

$A=\left[\begin{array}{ccccc}{b}_1& 0& 0& ...& 0\\ b_1& b_2& 0& ...& 0\\ b_1& b_2& b_3& ...& 0\\ ...& ...& ...& ...& 0\\ b_1& b_2& b_3& ...& b_s\end{array}\right]\text{and }\overline{A}=\left[\begin{array}{ccccc}0& 0& 0& ...& 0\\ {\overline{b}}_1& 0& 0& ...& 0\\ {\overline{b}}_1& {\overline{b}}_2& 0& ...& 0\\ ...& ...& ...& ...& 0\\ {\overline{b}}_1& {\overline{b}}_2& {\overline{b}}_3& ...& 0\end{array}\right].$

These methods are charaterized by the two vectors b and $\overline{b}$, with the obvious notation, and can be efficiently implemented as follows:

$\begin{array}{ll}\left[\begin{array}{c}\text{P}\\ \text{Q}\end{array}\right]\leftarrow \text{u},& \text{for}i=1,\dots ,s:\begin{array}{c}\text{P}\leftarrow \text{p}+h{b}_i{\text{A}}_{pq}\text{Q}\\ \text{Q}\leftarrow \text{Q}+h{\overline{b}}_i{\text{A}}_{qp}\text{P}\end{array},\\ {\text{E}}_{\Delta t}\text{u}\leftarrow \left[\begin{array}{c}\text{p}\\ \text{Q}\end{array}\right].\end{array}$

There are two, low accuracy methods which, nevertheless, deserve to be mentioned as they are considered as classic. The first is the so-called the Symplectic Euler method. It is associated with the vectors b = 1, $\overline{b}=1$, and is first-order accurate. It reads:

$\begin{array}{ll}{\text{p}}_{n+1}=& {\text{p}}_n+\Delta t{\text{A}}_{pq}{\text{Q}}_n,\\ {\text{Q}}_{n+1}=& {\text{Q}}_n+\Delta t{\text{A}}_{qp}{\text{P}}_{n+1}.\end{array}$

The second is the Störmer-Verlet method. It is associated to the vectors b = [0, 1]^⊤ and $\overline{b}={\left[1/2,1/2\right]}^{\top }$ and is second-order accurate. It can be obtained by applying the Symplectic Euler method twice by using the Strang splitting, that is,

$\begin{array}{l}\begin{array}{cc}{\text{p}}_{n+1/2}=& {\text{p}}_n+\frac{\Delta t}{2}{\text{A}}_{pq}{\text{Q}}_n,\\ {\text{Q}}_{n+1/2}=& {\text{Q}}_n+\frac{\Delta t}{2}{\text{A}}_{qp}{\text{P}}_{n+1/2},\end{array}\begin{array}{cc}{\text{Q}}_{n+1}=& {\text{Q}}_{n+1/2}+\frac{\Delta t}{2}{\text{A}}_{qp}{\text{P}}_{n+1/2}\\ {\text{p}}_{n+1},=& {\text{p}}_{n+1/2}+\frac{\Delta t}{2}{\text{A}}_{pq}{\text{Q}}_{n+1}.\end{array}\end{array}$

We can also write the method as follows:

$\begin{array}{ll}\begin{array}{cc}{\text{Q}}_{n+1}=& {\text{Q}}_n+\Delta t{\text{A}}_{qp}{\text{P}}_{n+1/2},\\ {\text{p}}_{n+3/2}=& {\text{p}}_{n+1/2}+\Delta t{\text{A}}_{pq}{\text{Q}}_{n+1},\end{array}\text{for}n=0,\dots ,N-1,& (5)\end{array}$

Where ${\text{p}}_{1/2}={\text{p}}_0+\frac{\Delta t}{2}{\text{A}}_{pq}{\text{Q}}_0$ and ${\text{p}}_N={\text{p}}_{N+1/2}-\frac{\Delta t}{2}{\text{A}}_{pq}{\text{Q}}_N$. This is precisely the time-marching scheme used in the well known Yee's scheme [15] for Maxwell's equations.

We end this section by noting that explicit partitioned Runge-Kutta methods of any (even) order of accuracy can be obtained, as was shown in 1990 by Yoshida [23].

3.5. Illustration

We end this section by illustrating the application of the implicit midpoint method (the lowest order, symplectic Runge-Kutta method) and the symplectic Euler method (the lowest order, symplectic partitioned Runge-Kutta method). We consider Example 2, introduced in Section 2.3. Both methods preserve the linear first integral exactly (see Theorem 3.1 for the implicit midpoint method and Theorem 3.3 for the symplectic Euler method), this is observed in Figure 2 where the approximate orbits lie on planes perpendicular to the kernel of J, the vector (0, 1, 1). The implicit midpoint preserves the Hamiltonian exactly (see Theorem 3.1) and the symplectic Euler computes a non-drifting approximation of it (see Theorem 3.2).

FIGURE 2

Figure 2. Two views of the approximate orbits of the Poisson system in Example 2 on the phase-space (p, q, r). The approximate orbits lie on planes with normal vector (0, 1, 1) so that q + r is a constant. The orbits in the dashed-green line were computed using the symplectic Euler method are ellipses and the ones in solid-red were computed using the implicit midpoint method and are circles. The blue sphere represents the constant Hamiltonian value corresponding to ${H}=3^2/2$.

4. Symplectic finite difference methods

The idea of combining symplectic time integrators with finite difference spatial discretizations can be traced back to the late 1980's with the work of Feng and his collaborators [10, 11]. For instance, in 1987, Feng and Qi [11] applied symplectic integrators to the discretization in space by central difference schemes of the linear wave equation. In 1988, Ge and Feng [10] combined symplectic time-marching methods with finite difference schemes for a linear hyperbolic system. In 1993, McLachlan [16] incorporated the Poisson structure from Olver [5] to the design of symplectic finite difference schemes, and applied it to the non-linear waves, the Schrödinger, and the KdV equations. This effort continues. See, for example, the study of multisymplectic finite difference methods [24].

On the other hand, numerical methods (which are essentially symplectic time-marching methods applied to finite difference space discretizations) had been proposed much earlier (than 1980's) when the concept of symplectic time-integrators did not exist or was not systematically studied. A prominent example is Yee's scheme proposed in 1966 [15]. This scheme is also known as the finite-difference time-domain (FDTD) method, for which the acronym was first introduced in Taflove [25]. Later on, Yee's scheme was studied in the multisymplectic framework [26, 27]. However, to the best of our knowledge, no work exists which attempts recasting Yee's scheme as a combination of symplectic time-marching methods with finite difference space discretizations.

In this section, instead of attempting to reach for maximal generality, we focus on the Yee's scheme [15]. We show that the spatial discretization of the scheme leads to a Hamiltonian ODE system, while the time-discretization of the scheme is nothing but the well known, symplectic, second-order accurate Störmer-Verlet method. As a result, all the conservation properties of Section 3 hold for Yee's scheme.

4.1. The time-dependent Maxwell equations

Let us begin by recalling the time-dependent Maxwell's equations in the domain Ω = [0, L₁] × [0, L₂] × [0, L₃]:

$\begin{array}{ll}ϵ\dot{\text{E}}& =\nabla \times \text{H}\text{ in }\Omega ,\\ \mu \dot{\text{H}}& =-\nabla \times \text{E}\text{ in }\Omega ,\end{array}$

With periodic boundary conditions, where E and H represent the electric and the magnetic fields, respectively. If we set u = (E, H)^T, we can rewrite the above equations into a more compact form:

$\begin{array}{ll}\dot{\text{u}}=\text{A}\text{u}=\text{J}\text{H}\text{u},& (\mathrm{6a})\end{array}$

where

$\begin{array}{ll}\text{J}=\left[\begin{array}{cc}0& {ϵ}^{-1}(\nabla \times ){\mu }^{-1}\\ -{\mu }^{-1}(\nabla \times ){ϵ}^{-1}& 0\end{array}\right]\text{ and }\text{H}=\left[\begin{array}{cc}ϵ& 0\\ 0& \mu \end{array}\right].& (\mathrm{6b})\end{array}$

Note Equation (6) is a Hamiltonian dynamical system with the triplet $({M},{H},\left\{·,·\right\})$, where ${M}$ is the phase space, ${H}$ is the Hamiltonian

$\mathscr{H}=\frac{1}{2}{\text{u}}^{\top }\text{H}\text{u}=\frac{1}{2}{\displaystyle {\int }_{\Omega }\left(ϵ|\text{E}{|}^2+\mu |\text{H}{|}^2\right)},$

and the Poisson bracket is defined by

$\{ℱ,G\}:={\displaystyle {\int }_{\Omega }{(\frac{\delta ℱ}{\delta \text{u}})}^{\top }}\text{J}\frac{\delta G}{\delta \text{u}}.$

4.2. The space discretization of Yee's scheme

We next consider the spatial discretization of Yee's scheme [15], for which the electric and the magnetic fields are defined on the following grid points:

$\begin{array}{l}{E}_{1,i,j+\frac{1}{2},k+\frac{1}{2}},E_{2,i+\frac{1}{2},j,k+\frac{1}{2}},E_{3,i+\frac{1}{2},j+\frac{1}{2},k},\\ H_{1,i+\frac{1}{2},j,k},H_{2,i,j+\frac{1}{2},k},H_{3,i,j,k+\frac{1}{2}}.\end{array}$

Namely, the magnetic field is defined on the edges and the electric field is defined on the faces of the cube element [iΔx, (i+1)Δx] × [jΔy, (j+1)Δy] × [kΔz, (k+1)Δz]. Let us denote by

$\begin{array}{l}{V}_E=V_{E_1}\times V_{E_2}\times V_{E_3},V_H=V_{H_1}\times V_{H_2}\times V_{H_3},\end{array}$

the approximation spaces for the electric field E and the magnetic field H, respectively, and let z∈V_E×V_H be the vector representing the discretized electric and magnetic fields. Then, the spatial discretization of Yee's scheme introduces the following Hamiltonian system of ODEs:

$\begin{array}{ll}\dot{\text{z}}=\text{J}\text{H}\text{z},& (\mathrm{7a})\end{array}$

where

$\begin{array}{lll}\text{J}& =\frac{1}{{(\Delta x\Delta y\Delta z)}^2}\left[\begin{array}{cc}0& {({\text{D}}_h^{ϵ})}^{-1}\text{cur}{\text{l}}_h^H{({\text{D}}_h^{\mu })}^{-1}\\ -{({\text{D}}_h^{\mu })}^{-1}\text{cur}{\text{l}}_h^E{({\text{D}}_h^{ϵ})}^{-1}& 0\end{array}\right],& (\mathrm{7b})\end{array}$

$\begin{array}{lll}\text{H}& =(\Delta x\Delta y\Delta z)\left[\begin{array}{cc}{\text{D}}_h^{ϵ}& 0\\ 0& {\text{D}}_h^{\mu }\end{array}\right].& (\mathrm{7c})\end{array}$

In the above equations,

$\begin{array}{l}{\text{D}}_h^{ϵ}:V_E\to V_E,{\text{D}}_h^{\mu }:V_H\to V_H,\end{array}$

are the discretizations of the multiplication operation by ϵ and μ, respectively. Moreover,

$\begin{array}{l}\text{cur}{\text{l}}_h^H:V_H\to V_E,\text{ cur}{\text{l}}_h^E:V_E\to V_H,\end{array}$

are two different finite difference discretizations of the ∇ × operator. These operators will be described in full detail in the last subsection. Here we want to emphasize the abstract structure of the scheme as well as the main properties of these operators which are contained in the following result. A detailed proof is presented in Section 4.5.

Proposition 4.1. The operators ${\text{D}}_h^{ϵ}$ and ${\text{D}}_h^{\mu }$ are symmetric and semi-positive definite. In addition, we have

$\begin{array}{l}{(\text{cur}{\text{l}}_h^H)}^{\top }=\text{cur}{\text{l}}_h^E.\end{array}$

Consequently, the operator J defined in Equation (7b) is anti-symmetric, and the operator H defined in Equation (7c) is semi-positive definite.

As a consequence of the above proposition, the ODE system (Equation 7a) is a Poisson dynamical system. In addition, thanks to the structure of J and H (see Equations 7b, 7c), we know its Hamiltonian is separable.

4.3. The fully discrete Yee's scheme

So we can now time-march this system of ODEs by using a symplectic scheme so that all the results of Section 3.4 hold.

In particular, if in the Equation (5), we replace q by E, P by H, A_qp by ${({\text{D}}_h^{ϵ})}^{-1}\text{cur}{\text{l}}_h^H$ and A_pq by $(-{({\text{D}}_h^{\mu })}^{-1}\text{cur}{\text{l}}_h^E)$, we obtain

$\begin{array}{ll}{\text{E}}_{n+1}={\text{E}}_n+\Delta t{(D_h^{ϵ})}^{-1}{\text{curl}}_h^H{\text{H}}_{n+1/2},& (8)\end{array}$

$\begin{array}{ll}{\text{H}}_{n+3/2}={\text{H}}_{n+1/2}-\Delta t{(D_h^{\mu })}^{-1}\text{cur}{\text{l}}_h^E{\text{E}}_n,& (9)\end{array}$

and we recover Yee's scheme [15].

4.4. Conservation laws

Finally, we identify some conservation quantities associated to Yee's scheme. By Proposition 3.3, this task reduces to identify

• The kernel of A^⊤. Suppose A^⊤v = 0, then v^⊤u is conserved in time.

• Symmetric matrix S such that SA+A^⊤S = 0. Then, the bilinear form u^⊤Su is conserved in time.

Before, we can present the main result, let us introduce some notation. For a given function f_{i, j, k}, we define a shifting operator as follows:

${\tau }_{(s_1,s_2,s_3)}(f_{i,j,k})=f_{i+s_1,j+s_2,k+s_3}.$

In the above formulation, we allow the indexes i, j, and k to be non-integers so that the operator τ_{(s1, s2, s3)} can be applied to f_{i, j, k} on non-integer grid points. With this notation, we define the discretized Hamiltonian energy function

$\begin{array}{l}\mathscr{H}(E_1,E_2,E_3,H_1,H_2,H_3)\\ :=\frac{1}{2}\Delta x\Delta y\Delta z\left({\left({\text{D}}_h^{ϵ}\text{E},\text{E}\right)}_{V_E}+{\left({\text{D}}_h^{\mu }\text{H},\text{H}\right)}_{V_H}\right)\\ :=\frac{1}{2}\Delta x\Delta y\Delta z{\displaystyle \sum _{i,j,k}({\tau }_{(0,\frac{1}{2},\frac{1}{2})}({ϵ}_{i,j,k}{E}_{1,i,j,k}^2)+{\tau }_{(\frac{1}{2},0,\frac{1}{2})}({ϵ}_{i,j,k}{E}_{2,i,j,k}^2)}\\ +{\tau }_{(\frac{1}{2},\frac{1}{2},0)}({ϵ}_{i,j,k}{E}_{3,i,j,k}^2))\\ +\frac{1}{2}\Delta x\Delta y\Delta z{\displaystyle \sum _{i,j,k}({\tau }_{(\frac{1}{2},0,0)}({\mu }_{i,j,k}{H}_{1,i,j,k}^2)+{\tau }_{(0,\frac{1}{2},0)}({\mu }_{i,j,k}{H}_{2,i,j,k}^2)}\\ +{\tau }_{(0,0,\frac{1}{2})}({\mu }_{i,j,k}{H}_{3,i,j,k}^2)).\end{array}$

We also introduce two gradient operators

$\begin{array}{ll}(i)& {\nabla }_{h,0}{\varphi }_{i,j,k}:={\text{T}}_h({\varphi }_{i,j,k}),\\ (ii)& {\nabla }_{h,\frac{1}{2}}{\varphi }_{i+\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2}}:={\text{T}}_h({\varphi }_{i+\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2}}),\end{array}$

Where T_h is the differencing operator given by

${\text{T}}_h:=\left(\frac{{\tau }_{(1,0,0)}-{\tau }_{(0,0,0)}}{\Delta x},\frac{{\tau }_{(0,1,0)}-{\tau }_{(0,0,0)}}{\Delta y},\frac{{\tau }_{(0,0,1)}-{\tau }_{(0,0,0)}}{\Delta z}\right).$

Note that the two gradients are defined on different grid spaces. The gradient ∇_{h, 0} is defined for function living on integer-grid points (i, j, k), while the gradient ${\nabla }_{h,\frac{1}{2}}$ is defined on half-grid points $(i+\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2})$. In addition, we remark that the range of ∇_{h, 0} is V_H, and the range of ${\nabla }_{h,\frac{1}{2}}$ is V_E.

Proposition 4.2 (Conservation of energy and electric/magnetic charges). For the Yee's semi-discrete scheme defined by Equation (7a), the energy ${H}$ is conserved in time. In addition, (the weak form of) the electric and magnetic charges

$\begin{array}{l}{{C}}_{\varphi }(E_1,E_2,E_3,H_1,H_2,H_3):={({\text{D}}_h^{ϵ}\text{E},{\nabla }_{h,2}{\varphi }_1)}_{V_E}+{({\text{D}}_h^{\mu }\text{H},{\nabla }_{h,1}{\varphi }_2)}_{V_H},\end{array}$

is conserved in time. Here, ϕ_i (i = 1, 2) are arbitrary test functions.

4.5. Proof of proposition 4.1 and 4.2

In this subsection, we first define the operators ${\text{D}}_h^{ϵ}$, ${\text{D}}_h^{\mu }$, $\text{cur}{\text{l}}_h^H$, and $\text{cur}{\text{l}}_h^E$. Then we prove Propositions 4.2 and 4.3.

We introduce two types of central difference operators ${\partial }_x^{h,0}$ and ${\partial }_x^{h,\frac{1}{2}}$, which are defined as

$\begin{array}{l}{({\partial }_x^{h,0}f)}_{i+\frac{1}{2},j,k}=\frac{{\tau }_{(1,0,0)}-{\tau }_{(0,0,0)}}{\Delta x}{f}_{i,j,k},\hfill \\ {({\partial }_x^{h,\frac{1}{2}}f)}_{i,j,k}=\frac{{\tau }_{(1/2,0,0)}-{\tau }_{(-1/2,0,0)}}{\Delta x}{f}_{i,j,k}.\hfill \end{array}$

Similarly, we can define ${\partial }_y^{h,0}$, ${\partial }_y^{h,\frac{1}{2}}$, ${\partial }_z^{h,0}$, and ${\partial }_z^{h,\frac{1}{2}}$. It is easy to observe that

$\begin{array}{l}{\partial }_{y,H_1}^{h,0}:V_{H_1}\to V_{E_3},{\partial }_{z,H_1}^{h,0}:V_{H_1}\to V_{E_2},\\ {\partial }_{x,H_2}^{h,0}:V_{H_2}\to V_{E_3},{\partial }_{z,H_2}^{h,0}:V_{H_2}\to V_{E_1},\\ {\partial }_{x,H_3}^{h,0}:V_{H_3}\to V_{E_2},{\partial }_{y,H_3}^{h,0}:V_{H_3}\to V_{E_1}.\end{array}$

Note that we have added the additional sub-indexes H_i to indicate the domain of the operators.

Lemma 4.1 (The anti-symmetry of central differences). We have

$\begin{array}{l}{\partial }_{y,E_3}^{h,\frac{1}{2}}=-{({\partial }_{y,H_1}^{h,0})}^{*}:V_{E_3}\to V_{H_1},{\partial }_{z,E_2}^{h,\frac{1}{2}}=-{({\partial }_{z,H_1}^{h,0})}^{*}:V_{E_2}\to V_{H_1},\\ {\partial }_{x,E_3}^{h,\frac{1}{2}}=-{({\partial }_{x,H_2}^{h,0})}^{*}:V_{E_3}\to V_{H_2},{\partial }_{z,E_1}^{h,\frac{1}{2}}=-{({\partial }_{z,H_2}^{h,0})}^{*}:V_{E_1}\to V_{H_2},\\ {\partial }_{x,E_2}^{h,\frac{1}{2}}=-{({\partial }_{x,H_3}^{h,0})}^{*}:V_{E_2}\to V_{H_3},{\partial }_{y,E_1}^{h,\frac{1}{2}}=-{({\partial }_{y,H_3}^{h,0})}^{*}:V_{E_1}\to V_{H_3}.\end{array}$

Proof. We will only prove that ${\partial }_{y,E_3}^{h,\frac{1}{2}}=-{({\partial }_{y,H_1}^{h,0})}^{*}:V_{E_3}\to V_{H_1}$. The rest is similar. Notice that

$\begin{array}{l}{\langle {\partial }_{y,H_1}^{h,0}{H}_1,T\rangle }_{V_{E_3}}={\displaystyle \sum _{i,j,k}}{({\partial }_{y,H_1}^{h,0}{H}_1)}_{i+\frac{1}{2},j+\frac{1}{2},k}{T}_{i+\frac{1}{2},j+\frac{1}{2},k}\hfill \\ ={\displaystyle \sum _{i,j,k}}\frac{H_{1,i+\frac{1}{2},j+1,k}-H_{1,i+\frac{1}{2},j,k}}{\Delta x}{T}_{i+\frac{1}{2},j+\frac{1}{2},k}\hfill \\ ={\displaystyle \sum _{i,j,k}}{H}_{1,i+\frac{1}{2},j,k}\frac{T_{i+\frac{1}{2},j-\frac{1}{2},k}-T_{i+\frac{1}{2},j+\frac{1}{2},k}}{\Delta x}\hfill \\ ={\langle H_1,-{\partial }_{y,E_3}^{h,\frac{1}{2}}T\rangle }_{V_{H_1}}.\hfill \end{array}$

This proves the claim.              ⃞

Now, we define two discrete curl operators:

$\begin{array}{l}{\text{curl}}_h^E:V_E\to V_H\\ (E_1,E_2,E_3)\mapsto \left[\begin{array}{ccc}0& -{\partial }_{z,E_2}^{h,\frac{1}{2}}& {\partial }_{y,E_3}^{h,\frac{1}{2}}\\ {\partial }_{z,E_1}^{h,\frac{1}{2}}& 0& -{\partial }_{x,E_3}^{h,\frac{1}{2}}\\ -{\partial }_{y,E_1}^{h,\frac{1}{2}}& {\partial }_{x,E_2}^{h,\frac{1}{2}}& 0\end{array}\right]\left[\begin{array}{c}{E}_1\\ E_2\\ E_3\end{array}\right],\end{array}$

and

$\begin{array}{l}{\text{curl}}_h^H:V_H\to V_E\\ (H_1,H_2,H_3)\mapsto \left[\begin{array}{ccc}0& -{\partial }_{z,H_2}^{h,0}& {\partial }_{y,H_3}^{h,0}\\ {\partial }_{z,H_1}^{h,0}& 0& -{\partial }_{x,H_3}^{h,0}\\ -{\partial }_{y,H_1}^{h,0}& {\partial }_{x,H_2}^{h,0}& 0\end{array}\right]\left[\begin{array}{c}{H}_1\\ H_2\\ H_3\end{array}\right].\end{array}$

We also introduce the multiplication operators ${\text{D}}_h^{ϵ}:V_E\to V_E$ such that

$\begin{array}{l}{\text{D}}_h^{ϵ}(E_{1,i,j+\frac{1}{2},k+\frac{1}{2}},E_{2,i+\frac{1}{2},j,k+\frac{1}{2}},E_{3,i+\frac{1}{2},j+\frac{1}{2},k})\\ =({ϵ}_{i,j+\frac{1}{2},k+\frac{1}{2}}{E}_{1,i,j+\frac{1}{2},k+\frac{1}{2}},{ϵ}_{i+\frac{1}{2},j,k+\frac{1}{2}}{E}_{2,i+\frac{1}{2},j,k+\frac{1}{2}},\\ {ϵ}_{i+\frac{1}{2},j+\frac{1}{2},k}{E}_{3,i+\frac{1}{2},j+\frac{1}{2},k}).\end{array}$

The operator ${\text{D}}_h^{\mu }:V_H\to V_H$ can be similarly defined.

Proof of Proposition 4.1 As we have seen, ${\text{D}}_h^{ϵ}$ and ${\text{D}}_h^{\mu }$ are semi-positive definite since ϵ and μ are non-negative. Hence, the matrix H is semi-positive definite by its definition (Equation 7c). Moreover, by Lemma 4.1, we have that ${(\text{cur}{\text{l}}_h^H)}^{*}=\text{cur}{\text{l}}_h^E$, and so, J is anti-symmetric. This completes the proof.                    ⃞

We are next going to prove Proposition 4.2. Before that, we will need a lemma.

Lemma 4.2 (The kernel of the finite difference curl_h operators). We have that

$\begin{array}{l}(\text{i}){\nabla }_{h,1}{\varphi }_{i,j,k}:={\text{T}}_h({\varphi }_{i,j,k})\\ {\text{belongs to the kernel of curl}}_h^H,\\ (\text{ii}){\nabla }_{h,2}{\varphi }_{i+\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2}}:={\text{T}}_h({\varphi }_{i+\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2}})\\ {\text{belongs to the kernel of curl}}_h^E,\end{array}$

Where T_h is the differencing operator given by

${\text{T}}_h:=\left(\frac{{\tau }_{(1,0,0)}-{\tau }_{(0,0,0)}}{\Delta x},\frac{{\tau }_{(0,1,0)}-{\tau }_{(0,0,0)}}{\Delta y},\frac{{\tau }_{(0,0,1)}-{\tau }_{(0,0,0)}}{\Delta z}\right)$

Proof. We prove the lemma by a direct computation. Note that

$\begin{array}{l}cur{l}_h^H◦{\nabla }_{h,1}{\varphi }_{i,j,k}=\left[\begin{array}{ccc}0& -{\partial }_{z,H_2}^{h,0}& {\partial }_{y,H_3}^{h,0}\\ {\partial }_{z,H_1}^{h,0}& 0& -{\partial }_{x,H_3}^{h,0}\\ -{\partial }_{y,H_1}^{h,0}& {\partial }_{x,H_2}^{h,0}& 0\end{array}\right]\left[\begin{array}{c}\frac{{\tau }_{(1,0,0)}-{\tau }_{(0,0,0)}}{\Delta x}\\ \frac{{\tau }_{(0,1,0)}-{\tau }_{(0,0,0)}}{\Delta y}\\ \frac{{\tau }_{(0,0,1)}-{\tau }_{(0,0,0)}}{\Delta z}\end{array}\right]{\varphi }_{i,j,k}\\ =\left[\begin{array}{c}\frac{({\tau }_{(0,1,1)}-{\tau }_{(0,1,0)})-({\tau }_{(0,0,1)}-{\tau }_{(0,0,0)})}{\Delta z\Delta y}-\frac{({\tau }_{(0,1,1)}-{\tau }_{(0,0,1)})-({\tau }_{(0,1,0)}-{\tau }_{(0,0,0)})}{\Delta y\Delta z}\\ \frac{({\tau }_{(1,0,1)}-{\tau }_{(0,0,1)})-({\tau }_{(1,0,0)}-{\tau }_{(0,0,0)})}{\Delta x\Delta z}-\frac{({\tau }_{(1,0,1)}-{\tau }_{(1,0,0)})-({\tau }_{(0,0,1)}-{\tau }_{(0,0,0)})}{\Delta x\Delta z}\\ \frac{({\tau }_{(1,1,0)}-{\tau }_{(1,0,0)})-({\tau }_{(0,1,0)}-{\tau }_{(0,0,0)})}{\Delta y\Delta x}-\frac{({\tau }_{(1,1,0)}-{\tau }_{(0,1,0)})-({\tau }_{(1,0,0)}-{\tau }_{(0,0,0)})}{\Delta x\Delta y}\end{array}\right]x\\ {\varphi }_{i,j,k}=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right].\end{array}$

Hence, Lemma 4.2-(i) holds. Lemma 4.2-(ii) can be proven in a similar manner. This completes the proof.                     ⃞

Proof of Proposition 4.2 Note that we have

$\begin{array}{l}{H}(\text{z})=\frac{1}{2}{\text{z}}^{\top }\text{H}\text{z}.\end{array}$

Hence, to prove that ${H}$ is preserved in time, by Proposition 3.3, we only need to show that HA+A^⊤H = 0, which is true since HA+A^⊤H = HJH+H^⊤J^⊤H, J^⊤ = −J, and H^⊤ = H. Thus, the Hamiltonian ${H}$ is preserved in time on the discrete orbits of the numerical method.

Next, let us prove the conservation of the electric and magnetic charges. By Proposition 3.3, this task reduces to identify the kernel of A^⊤. Since

$\begin{array}{l}{\text{A}}^{\top }=\left[\begin{array}{cc}0& -\text{cur}{\text{l}}_h^H{({\text{D}}_h^{\mu })}^{-1}\\ cur{l}_h^E{({\text{D}}_h^{ϵ})}^{-1}& 0\end{array}\right],\end{array}$

by Lemma 4.2, we have

$\begin{array}{l}\left[\begin{array}{c}{\text{D}}_h^{ϵ}{\nabla }_{h,2}\varphi \\ {\text{D}}_h^{\mu }{\nabla }_{h,1}\psi \end{array}\right]\in ker({\text{A}}^{\top }).\end{array}$

This completes the proof.                    ⃞

5. Symplectic-Hamiltonian finite element methods

This section presents the recent development of the so-called symplectic-Hamiltonian finite element methods for linear Hamiltonian partial differential equations. First, we discuss our results on energy-preserving hybridizable discontinuous Galerkin (HDG) methods for the scalar wave equation [1, 28], and then describe our further contributions to finite element discretizations of the linear elastic equation [2] and the Maxwell equations [3, 9]. We analyze in detail four finite element discretizations of the linear scalar wave (Equation 2) based on standard Galerkin methods, mixed methods, and HDG methods, and prove their structure-preserving properties.

5.1. Notation

Let us first introduce some standard finite element notation. For a given bounded domain Ω⊂ℝ^d, let ${{T}}_h=\left\{K\right\}$ be a family of regular triangulations of $\bar{\Omega }$, with mesh-size parameter h the maximum diameter over the elements. We also assume for simplicity the elements are only simplices. For a domain D∈ℝ^d we denote by (·, ·)_D the L² inner product for scalar, vector and tensor valued functions

${(w,v)}_D:={\displaystyle {\int }_{D}w}v,{(w,v)}_D:={\displaystyle {\int }_D}w·v,{(\underset{\_}{w},\underset{\_}{v})}_D:={\displaystyle {\int }_D}\underset{\_}{w}:\underset{\_}{v},$

for (w, v)∈L²(D), (w, v)∈L²(D)^d, and (w, v)∈L²(D)^d×d. We extend these definitions for the inner product in (d−1)-dimensional domains. For discontinuous Galerkin methods, we define by $\partial {{T}}_h$ the set of all element boundaries ∂K, for $K\in {{T}}_h$, and by ${{F}}_h$ the set of all the faces F of the triangulation ${{T}}_h$. Inner product definitions are extended to these sets by

${(w,v)}_{{{T}}_h}:={\displaystyle \sum _{K\in {{T}}_h}}{(w,v)}_K,{\langle w,v\rangle }_{\partial {{T}}_h}:={\displaystyle \sum _{K\in {{T}}_h}}{\langle w,v\rangle }_{\partial K}.$

for properly defined scalar-valued functions w, v. Similar definitions are given for vector- and tensor-valued functions.

5.2. The continuous Galerkin method

We now present a standard finite element discretization of the linear scalar wave equation using H¹-conforming approximations for the displacement and velocity variables. Let V_h be a continuous piece-wise polynomial subspace of $H_0^1(\Omega )$, then the semi-discrete Galerkin method, corresponding to the primal formulation of the wave (Equation 3), reads as follows: Find (u_h, v_h)∈V_h×V_h such that

$\begin{array}{ll}{(\rho {\dot{u}}_h,w\text{)}}_{\Omega }=\text{ (}\rho v_h,w{)}_{\Omega },\forall w\in V_h,& (\mathrm{10a})\end{array}$

$\begin{array}{ll}{\text{(}\rho {\dot{v}}_h,w)}_{\Omega }=-{(\kappa \nabla u_h,\nabla w)}_{\Omega }+{(f,w)}_{\Omega },\forall w\in V_h.& (\mathrm{10b})\end{array}$

The system is initialized by projections of the initial conditions u₀ and v₀ onto the finite-dimensional space V_h. It is clear that the semi-discrete method inherits the Hamiltonian structure of the continuous (Equation 3) with triplet $({{M}}_h,{{H}}_h,\left\{·,·\right\})$, where the discrete phase-space is ${{M}}_h=V_h\times V_h$, and the Hamiltonian and Poisson brackets are the respective restrictions to this space, i.e.

$\begin{array}{l}{\mathscr{H}}_h[u_h,v_h]=\frac{1}{2}{(\rho v_h,v_h)}_{\Omega }+\frac{1}{2}{(\kappa \nabla u_h,\nabla u_h)}_{\Omega }-{(f,u_h)}_{\Omega },\hfill \\ \{F,G\}={\left({\rho }^{-1}\frac{\delta F}{\delta (u_h,v_h)}J\frac{\delta G}{\delta (u_h,v_h)}\right)}_{\Omega },\hfill \end{array}$

for linear functionals F = F[u_h, v_h], G = G[u_h, v_h], for u_h, v_h∈V_h, and where denotes the canonical anti-symmetric structure matrix $\text{J}\in {ℝ}^{2\text{dim}(V_h)}$.

Equivalently, we can formulate the method in its matrix form based on a given ρ-orthonormal basis {ϕ_i} of the space V_h and define the matrix and vector

${({\text{S}}_{\kappa })}_{ij}:={(\kappa \nabla {\varphi }_i,\nabla {\varphi }_j)}_{\Omega }{\text{f}}_i={(f,{\varphi }_i)}_{\Omega }.$

Then, the evolution variables are represented in this basis by means of the coefficients u,v, i.e., $u_h(x,t)=\sum _i{u}_i(t){\varphi }_i(x)$ and $v_h(x,t)=\sum _i{v}_i(t){\varphi }_i(x)$. We can write the Hamiltonian functional as a function of the coefficients y = (u,v)^⊤ as

${H}\left[u_h,v_h\right]=\frac{1}{2}{v}^{\top }v+\frac{1}{2}{u}^{\top }{\text{S}}_{\kappa }u-u^{\top }\text{f}=\frac{1}{2}{\text{y}}^{\top }\left[\begin{array}{cc}{S}_{\kappa }& 0\\ 0& Id\end{array}\right]\text{y}-{\text{y}}^{\top }\left[\begin{array}{c}\text{f}\\ 0\end{array}\right].$

Therefore, computing the gradient of the Hamiltonian respect to the coefficients y = (u,v)^⊤, we conclude that the standard Galerkin method (Equation 10) is equivalent to

$\left[\begin{array}{c}\dot{u}\\ \dot{v}\end{array}\right]=\left[\begin{array}{cc}0& Id\\ -\text{Id}& 0\end{array}\right]\left[\begin{array}{c}{\text{S}}_{\kappa }u-\text{f}\\ v\end{array}\right]=\left[\begin{array}{cc}0& Id\\ -\text{Id}& 0\end{array}\right]\nabla {H},$

from where its Hamiltonian structure is evident. Note that the resulting system of differential equations is a canonical Hamiltonian system.

The semi-method is then discretized in time using a symplectic time-marching scheme. Here we write down the discretization by an s-stage explicit partitioned Runge-Kutta scheme with coefficients b and $\bar{b}$. The fully discrete scheme is written in terms of the variables yⁿ = (uⁿ,vⁿ)≈(u(tⁿ),v(tⁿ)), for tⁿ = nΔt, n∈ℕ, and time-step Δt, by the iterations

$\begin{array}{l}{\left[\text{P},\text{Q}\right]}^{\top }\leftarrow {\text{y}}^n=(u^n,v^n),\\ \text{P}\leftarrow \text{P}+b_i\Delta t\text{Q},\text{Q}\leftarrow \text{Q}-{\bar{b}}_i\Delta t{\text{S}}_{\kappa }\text{P},\text{for}1\le i\le s,\\ (u^{n+1},v^{n+1})={\text{y}}^{n+1}\leftarrow {\left[\text{P},\text{Q}\right]}^{\top }.\end{array}$

5.3. Mixed methods

In Kirby and Kieu [19] the authors introduce a mixed method for the scalar wave (Equation 4) and prove its Hamiltonian structure, here we review their results and prove that the resulting system of differential equations is a Poisson system.

Let $W_h\in L^2(\Omega )$ and v_h⊂H(div;Ω) be mixed finite element spaces and define the semi-discrete mixed method as follows: Find (v_h, q_h)∈W_h×v_h solution of

$\begin{array}{ll}{(\rho {\dot{v}}_h,w)}_{\Omega }=\text{ (}\nabla ·{\text{q}}_h,w{)}_{\Omega }+{(f,w)}_{\Omega },\forall w\in W_h,& (\mathrm{11a})\end{array}$

$\begin{array}{ll}({\kappa }^{-1}{\stackrel{.}{\text{q}}}_h,\text{r})=-{(v_h\nabla ·\text{r})}_{\Omega },\forall \text{r}\in {\text{V}}_h.& (\mathrm{11b})\end{array}$

The system is initialized by (v_h(0), q_h(0)), where v_h(0) is taken as the L²-projection of the initial data v₀ onto W_h, and q_h(0)∈v_h is solution of the problem

$\begin{array}{l}{(\nabla ·{\text{q}}_h(0),w)}_{\Omega }={(\nabla ·\kappa \nabla u_0,w)}_{\Omega },\forall w\in W_h,\hfill \\ {({\kappa }^{-1}{\text{q}}_h(0),\text{r})}_{\Omega }+{(u_h(0),\nabla ·\text{r})}_{\Omega }=0,\forall \text{r}\in {\text{V}}_h.\hfill \end{array}$

Moreover, the displacement approximation can be computed by $u_h(t)=u_h(0)+\underset{0}{\overset{t}{\int }}{v}_h(s)ds$. The semi-discrete method is Hamiltonian with triplet $({{M}}_h,{{H}}_h,\left\{·,·\right\})$ with ${{M}}_h=W_h\times {\text{V}}_h$,

$\begin{array}{l}\mathscr{H}[v_h,{\text{q}}_h]=\frac{1}{2}(\rho v_h,v_h)+\frac{1}{2}{({\kappa }^{-1}{\text{q}}_h,{\text{q}}_h)}_{\Omega }-{({\text{f}}_{div},{\text{q}}_h)}_{\Omega }\hfill \\ \{F,G\}={\left(\frac{\delta F}{\delta (v_h,{\text{q}}_h)}\text{J}\frac{\delta G}{\delta (v_h,{\text{q}}_h)}\right)}_{\Omega },\hfill \\ \text{ J}=\left[\begin{array}{cc}0& -({\rho }^{-1}\nabla ·)°\kappa \\ (\kappa \nabla )°{\rho }^{-1}& 0\end{array}\right]\hfill \end{array}$

Let {ϕ_i} be a ρ-orthonormal basis of W_h and {ψ_i} be a κ⁻¹-orthonormal basis of v_h, and denote by v and q the coefficients associated to the solution of system (Equation 11), namely, $v_h(x,t)=\sum _i{v}_i(t){\varphi }_i(x)$ and ${\text{q}}_h(x,t)=\sum _i{q}_i(t){\text{ψ}}_i(x)$, and define the matrix

${\text{B}}_{ik}:={(\nabla ·{\text{ψ}}_k,{\varphi }_i)}_{\Omega }.$

We write the Hamiltonian functional in term of the coefficients (v,q)^⊤ = y by

${H}\left[v_h,{\text{q}}_h\right]=\frac{1}{2}{v}^{\top }v+\frac{1}{2}{q}^{\top }q+q^{\top }{f}_{\text{div}}=\frac{1}{2}{\text{y}}^{\top }\left[\begin{array}{cc}Id& 0\\ 0& Id\end{array}\right]\text{y}+{\text{y}}^{\top }\left[\begin{array}{c}0\\ f_{\text{div}}\end{array}\right]$

and thus the matrix system equivalent to the method (Equation 11) is as follows

$\begin{array}{ll}\left[\begin{array}{c}\dot{v}\\ \dot{q}\end{array}\right]& =\left[\begin{array}{cc}0& B\\ -{\text{B}}^{\top }& 0\end{array}\right]\left[\begin{array}{c}v\\ q\end{array}\right]=\left[\begin{array}{cc}0& B\\ -{\text{B}}^{\top }& 0\end{array}\right]\left[\begin{array}{c}\partial H/\partial v\\ \partial H/\partial q\end{array}\right].\end{array}$

Finally, the time-discretization is carried out by a symplectic time integrator. For instance, consider the Butcher array with coefficients c,b,A of a symplectic Runge-Kutta method, then the evolution system for the vector variable yⁿ = (^vn,qn)⊤≈(v(^{tn),q(tn))⊤}, for tⁿ = nΔt, n∈ℕ and time-step Δt,

$Q(\Delta t\text{A}){\text{y}}^{n+1}=P(\Delta t\text{A}){\text{y}}^n,\text{where }\text{A}=\left[\begin{array}{cc}0& \text{B}\\ -{\text{B}}^{\top }& 0\end{array}\right]$

and where P(z) = det(I−zA+zeb^⊤) and Q(z) = det(I−zA).

5.4. Hybridizable discontinuous Galerkin methods

Note that both finite element discretizations introduced above, the continuous Galerkin and the mixed method, inherit the Hamiltonian structure property of the continuous equations due to the conformity of their finite element sub-spaces. Non-conforming finite element discretizations, such as discontinuous Galerkin methods, present more challenges. Here we discuss two hybridizable discontinuous Galerkin schemes for the approximation of solutions of the linear scalar wave equation, the first one a dissipative scheme and the second one a method that inherits the Hamiltonian property.

5.4.1. A dissipative HDG scheme

Let us consider hybridizable discontinuous Galerkin methods for the formulation of velocity-gradient variables [29, 30]. We define the discontinuous element spaces

$\begin{array}{l}{W}_h=\{w\in L^2(\Omega ):w{|}_K\in W(K),\forall K\in {{T}}_h\},\hfill \\ {\text{V}}_h=\{\text{r}\in L^2{(\Omega )}^d:\text{r}{|}_K\in \text{V}(K),\hfill \\ \forall K\in {{T}}_h\},\hfill \\ M_h=\{\mu \in L^2({{F}}_h):{\mu }_F\in M(F),\forall F\in {{F}}_h\},\hfill \end{array}$

Where W(K), v(K), M(F) are local (polynomial) spaces. The semi-discrete HDG method then is as follows: Find $(v_h,{\text{q}}_h,{\hat{v}}_h)\in W_h\times {\text{V}}_h\times M_h$ such that

$\begin{array}{lll}{(\rho {\dot{v}}_h,w)}_{{{T}}_h}& =-{({\text{q}}_h,\nabla w)}_{{{T}}_h}+{\langle {\hat{\text{q}}}_h·\text{n},w\rangle }_{\partial {{T}}_h}+{(f,w)}_{{{T}}_h},& (\mathrm{12a})\end{array}$

$\begin{array}{ll}\forall w\in W_h,& (\mathrm{12b})\end{array}$

$\begin{array}{llll}{({\kappa }^{-1}{\dot{\text{q}}}_h,\text{r})}_{{{T}}_h}& =-{(v_h,\nabla ·\text{r})}_{{{T}}_h}+{\langle {\hat{v}}_h,\text{r}·\text{n}\rangle }_{\partial {{T}}_h},& \forall \text{r}\in {\text{V}}_h,& (\mathrm{12c})\end{array}$

$\begin{array}{llll}{\langle {\hat{\text{q}}}_h·\text{n},\mu \rangle }_{\partial {{T}}_h}& =0,& \mu \in M_h,& (\mathrm{12d})\end{array}$

$\begin{array}{llll}{\hat{\text{q}}}_h·\text{n}& :={\text{q}}_h·\text{n}-\tau (v_h-{\hat{v}}_h),& \text{on}\partial {{T}}_h.& (12e)\end{array}$

The formulation is no longer Hamiltonian, due to the definition of the numerical traces. In fact, we can prove that the discrete energy defined by

${ℰ}_h:=\frac{1}{2}{(\rho v_h,v_h)}_{{{T}}_h}+\frac{1}{2}{({\kappa }^{-1}{\text{q}}_h,{\text{q}}_h)}_{{{T}}_h}.$

is dissipative,

${\stackrel{.}{ℰ}}_h=-{\langle \tau (v_h-{\hat{v}}_h),v_h-{\hat{v}}_h\rangle }_{\partial {{T}}_h}+{(f,v_h)}_{{{T}}_h}.$

5.4.2. A symplectic Hamiltonian HDG scheme

In Sánchez et al. [1] we rewrote the method using the primal formulation (Equation 3), that is, we reintroduce the displacement variables and compute a steady-state approximation q_h. The semi-discrete method is find $(u_h,v_h,{\text{q}}_h,{\hat{u}}_h)\in W_h\times W_h\times {\text{V}}_h\times M_h$ such that

$\begin{array}{ll}(\rho {\dot{u}}_{h}w){{T}}_h=(\rho v_{h}w,){{T}}_h,\forall w\in W_h,& (\mathrm{13a})\end{array}$

$\begin{array}{ll}(\rho {\dot{v}}_h,w){{T}}_h=-(\text{q},\nabla w){{T}}_h+{\langle {\widehat{\text{q}}}_h·\text{n},w\rangle }_{\partial {{T}}_h}+\text{ (}f,w,{)}_{{{T}}_h},\forall w\in W_h,& (\mathrm{13b})\end{array}$

$\begin{array}{ll}{({\kappa }^{-1}\text{q},\text{r})}_{{{T}}_h}=-{(u_h,\nabla ·\text{r})}_{{{T}}_h}+{\langle {\widehat{u}}_h\text{r}·\text{n}\rangle }_{\partial {{T}}_h},\forall \text{r}\in {\text{V}}_h,& (\mathrm{13c})\end{array}$

$\begin{array}{ll}{\langle {\widehat{\text{q}}}_h·\text{n},\mu \rangle }_{{{T}}_h\backslash \Gamma }=0,\mu \in M_h,& (\mathrm{13d})\end{array}$

$\begin{array}{llll}{\hat{\text{q}}}_h·\text{n}& :={\text{q}}_h·\text{n}-\tau (u_h-{\hat{u}}_h),& \text{on }\partial {{T}}_h.& (13e)\end{array}$

This system is Hamiltonian with triplet $({{M}}_h,{{H}}_h,\left\{·,·\right\})$, where the phase space is ${{M}}_h=W_h\times W_h$, {·, ·} is the canonical Poisson bracket and ${{H}}_h$ is the discrete Hamiltonian given by

${{H}}_h\left[u_h,v_h\right]=\frac{1}{2}{(v_h,v_h)}_{{{T}}_h}+\frac{1}{2}{({\kappa }^{-1}{\text{q}}_h,{\text{q}}_h)}_{{{T}}_h}$

$+\frac{1}{2}{\langle \tau (u_h-{\hat{u}}_h),u_h-{\hat{u}}_h\rangle }_{\partial {{T}}_h}-{(f,u_h)}_{{{T}}_h}.$

Let {ϕ_i} the ρ-orthonormal basis of W_h, {ψ_k} the κ⁻¹-orthonormal basis of v_h, and {η_m} a basis of M_h. Define the matrices

$\begin{array}{l}{\text{C}}_{km}:={\langle {\text{ψ}}_k·\text{n},{\eta }_m\rangle }_{\partial {{T}}_h},{\text{B}}_{ik}:={(\nabla ·{\text{ψ}}_k,{\varphi }_i)}_{{{T}}_h},\\ {({\text{S}}_{\tau })}_{ij}:={\langle \tau {\varphi }_i,{\varphi }_j\rangle }_{\partial {{T}}_h},{({\text{E}}_{\tau })}_{im}:={\langle \tau {\varphi }_i,{\eta }_m\rangle }_{\partial {{T}}_h},\\ {({\text{G}}_{\tau })}_{mn}:={\langle \tau {\eta }_m,{\eta }_n\rangle }_{\partial {{T}}_h}.\end{array}$

Then, the variables are written in terms of the coefficients vectors $u,v,q,\hat{u}$ as follows

$\begin{array}{l}{u}_h(x,t)={\displaystyle \sum _i}{u}_i(t){\varphi }_i(x),v_h(x,t)={\displaystyle \sum _i}{v}_i(t){\varphi }_i(x),\\ {\text{q}}_h(x,t)={\displaystyle \sum _k}{q}_k(t){\text{ψ}}_k(x),{\hat{u}}_h(x,t)={\displaystyle \sum _m}{\hat{u}}_m(t){\eta }_m(x).\end{array}$

The Hamiltonian is then rewritten in terms of the coefficients as

${{H}}_h\left[u_h,v_h\right]=\frac{1}{2}{v}^{\top }v+\frac{1}{2}{q}^{\top }q+\left[\begin{array}{c}{u}^{\top },{\hat{u}}^{\top }\end{array}\right]\left[\begin{array}{cc}{\text{S}}_{\tau }& -{\text{E}}_{\tau }\\ -{\text{E}}_{\tau }^{\top }& {\text{G}}_{\tau }\end{array}\right]\left[\begin{array}{c}u\\ \hat{u}\end{array}\right]-u^{\top }f,$

and where the coefficients solve

$\begin{array}{l}\left[\begin{array}{c}\dot{u}\\ \dot{v}\end{array}\right]=\left[\begin{array}{cc}0& Id\\ -\text{Id}& 0\end{array}\right]\left[\begin{array}{c}{\text{S}}_{\tau }u-\text{B}q-{\text{E}}_{\tau }\hat{u}\\ v\end{array}\right]+\left[\begin{array}{c}0\\ f\end{array}\right]=\left[\begin{array}{cc}0& Id\\ -\text{Id}& 0\end{array}\right]\left[\begin{array}{c}\partial {H}/\partial u\\ \partial {H}/\partial v\end{array}\right].\end{array}$

In effect, it is clear that $\partial {H}/\partial v=v$, and

$\begin{array}{l}\frac{\partial \mathscr{H}}{\partial u}={\text{S}}_{\tau }u+\frac{\partial }{\partial u}\left(\frac{1}{2}{q}^{\top }q-{\widehat{u}}^{\top }{\text{E}}_{\tau }^{\top }u+\frac{1}{2}{\widehat{u}}^{\top }{\text{G}}_{\tau }\widehat{u}\right)-f\\ ={\text{S}}_{\tau }u+\frac{\partial }{\partial u}\left(\frac{1}{2}{q}^{\top }q-{\widehat{u}}^{\top }{\text{C}}^{\top }q-\frac{1}{2}{\widehat{u}}^{\top }{\text{G}}_{\tau }\widehat{u}\right)-f.\end{array}$

We rewrite (Equations 13c, 13d) in matrix form, take derivative with respect to u, and multiply on the left by $\left[q^{\top },-{\hat{u}}^{\top }\right]$ we obtain

$\left[q^{\top },-{\hat{u}}^{\top }\right]\left[\begin{array}{cc}Id& -\text{C}\\ {\text{C}}^{\top }& {\text{G}}_{\tau }\end{array}\right]\left[\begin{array}{c}\partial q/\partial u\\ \partial \hat{u}/\partial u\end{array}\right]=\left[q^{\top },-{\hat{u}}^{\top }\right]\left[\begin{array}{c}-{\text{B}}^{\top }\\ {\text{E}}_{\tau }^{\top }\end{array}\right],$

which can be expressed as a column vector

$\frac{\partial }{\partial u}\left(\frac{1}{2}{q}^{\top }q-{\hat{u}}^{\top }{\text{C}}^{\top }q-\frac{1}{2}{\hat{u}}^{\top }{\text{G}}_{\tau }\hat{u}\right)=-\text{B}q-{\text{E}}_{\tau }\hat{u}.$

Therefore, this proves that

$\frac{\partial {H}}{\partial u}={\text{S}}_{\tau }u-\text{B}q-{\text{E}}_{\tau }\hat{u}-f$

and thus the Hamiltonian form of Equation (13) follows.

5.5. A numerical example

We manufacture a numerical example to test the energy-conserving properties of the schemes presented in this section. Consider the rectangular domain with an obstacle Ω depicted in Figure 3. We solve the scalar wave equation with Homogeneous Neumann boundary conditions at all domain boundaries and initial conditions set as $u_0(x,y)=exp(-0.5{(x-4)}^2/{(1/5)}^2)$ and $v_0(x,y)=-25(x-4)exp(-0.5{(x-4)}^2/{(1/5)}^2)$.

FIGURE 3

Figure 3. Computational mesh of the rectangular domain Ω with an obstacle (second row, right) used for all the computation. The first row and second (left and middle) show snapshots at t = 0, 2, 4, 6, 8 of the velocity variables v_h approximated by the Hamiltonian HDG method (Equation 13) with polynomial approximation spaces of order 3 and the implicit midpoint.

In Figure 3, we present the computational mesh of the domain used in our computations. We compare six numerical schemes:

• CG-symp: Continuous Galerkin with polynomial order 3, and implicit midpoint scheme (symplectic).

• CG-diss: Continuous Galerkin with polynomial order 3 and singly DIRK method of order 3 (nonsymplectic), see Table 5.

• Mixed-symp: Mixed method with Raviart-Thomas spaces of order 3, and Störmer-Verlet scheme (symplectic PRK method of order 2).

• Mixed-diss: Mixed method with Raviart-Thomas spaces of order 3 and singly DIRK method of order 3 (nonsymplectic), see Table 5.

• HDG-symp: Hamiltonian HDG method (Equation 13) with polynomial order 3, and implicit midpoint scheme (symplectic).

• HDG-diss: HDG method (Equation 12) (non-Hamiltonian) with polynomial order 3, and implicit midpoint scheme (symplectic).

TABLE 5

Table 5. Butcher array for a singly diagonally implicit Runge-Kutta method, with $\gamma =(1+\sqrt{3})/2$.

We illustrate the evolution of the numerical approximation given by the Hamiltonian HDG method with polynomial approximations of order 3 and implicit midpoint in Figure 3. Snapshots are presented at times t = 0, 2, 4, 6, 8.

In Figure 4 we present a comparison of the numerical schemes computing the physical quantities energy (see Table 1), linear momentum, and the norm of the pseudomomentum for the six methods. In the first column, we plot the energy loss, that is, $|{ℰ}_h(0)-{ℰ}_h(t^n)|$, for the energy ${ℰ}_h(t^n)$ computed with each numerical scheme at the time tⁿ. In the second column, we plot the linear momentum loss and, in the third column the evolution of the pseudomomentum. We observe the improved approximation of the energy and pseudomomentum of the symplectic Hamiltonian finite element methods. The linear momentum, which is a first integral of the system, is preserved for each of the numerical schemes.

FIGURE 4

Figure 4. Plots of the approximated physical quantities: Energy loss (left), linear momentum loss (middle), and pseudomomentum. The first row shows a comparison of the continuous Galerkin methods using the implicit midpoint (dot-dashed blue line) and an sDIRK method of order 3 (dashed red line). The second row shows a comparison of the Mixed method using the Raviart-Thomas element with a polynomial approximation of order 3, using the Störmer-Verlet scheme and the sDIRK method of order 3. The third row shows a comparison of the Hamiltonian HDG scheme (Equation 13) and the dissipative HDG scheme (Equation 12), for both we use the implicit midpoint.

6. Ongoing work

As we have seen, the application of a symplectic time-marching method to a Hamiltonian system of ODEs can guarantee the preservation of all its linear and quadratic invariants. On the other hand, to obtain such a system of ODEs, one uses space discretizations which do not guarantee the preservation of all the linear and quadratic invariants of the original Hamiltonian PDEs! Indeed, as far as we can tell, it is not well understood how to obtain all the discrete versions of those invariants for any given space discretization, by finite difference or finite element methods. In fact, although it is almost automatic how to find discrete versions of the energy, the discrete equivalents of other conserved quantities, like the Lipkin's invariants for electromagnetism, for example, remain elusive. This topic constitutes the subject of ongoing work.

Recently, a new class of symplectic Discontinuous Galerkin methods were found whose distinctive feature is the use of time operators for defining their numerical traces, see Cockburn et al. [9]. Work to find how useful this type of methods are is under way.

The combination of Galerkin space discretizations with symplectic time-marching methods to a variety of systems of nonlinear Hamiltonian problems, including the Schrödinger and KdV equations, finite deformation elasticity, water waves and nonlinear wave equations, is also the subject of ongoing work.

Author contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

BC was supported in part by the NSF through grant no. DMS-1912646. MS was partially supported by FONDECYT Regular grant no. 1221189 and by Centro Nacional de Inteligencia Artificial CENIA, FB210017, Basal ANID Chile.

Acknowledgments

The authors would like to express their gratitude to EC for the invitation to write this article. The authors also want to thank Peter Olver for many useful clarifications on the various concepts used for the systems under consideration and for bringing to their attention the Olver [5, Example 4.36] and Walter [32] about conservation laws for the scalar wave equation.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

1. ^What we are denoting by J was originally denoted by J⁻¹. We decided not to keep such notation because it suggests the invertibility of J and this property does not hold in many interesting cases. For example, when extending these definitions from ODEs to PDEs, J⁻¹ naturally becomes a non-invertible differential operator, as already noted in the pioneering 1988 work by Ge and Feng [10, Section 3]. This prompted them to simply drop the notation J⁻¹. No wonder most authors dealing with Hamiltonian PDEs, see, for example, the 1993 book by Olver [5], automatically replace the original J⁻¹ by J, just as we are doing it here.

2. ^This holds when the functional ${F}$ only depends on the phase variable u. If the functional ${F}$ depends also on time, that is, if ${F}={F}(\text{u},t)$, the above calculation must be modified to read $\frac{d}{dt}{F}=\left\{{F},{H}\right\}+{\partial }_t{F}$.

3. ^This generalization is so simple that, in our first reading, the difference between a “Hamiltonian” and a “Poisson” dynamical system completely escaped us.

4. ^A small remark on the matrix A: = JH is in order. For the canonical structure matrix J, it is standard to say that the matrix A is Hamiltonian if J^⊤A = H is symmetric. If J is only required to be invertible, one could say that A is a Hamiltonian matrix if J⁻¹A = H is symmetric. Finally, if J is not invertible, one is tempted to say that A is a “Poissonian” matrix if AJ = JHJ is symmetric. We did not find this terminology in the current literature and so dissuaded ourselves to introduce it here, especially because we do not need it in a significant way.

5. ^This is an example of a Casimir: a function C whose gradient ∇C lies in the kernel of J. The functional C is then conserved for any Poisson system with structure matrix J. A Casimir is also called a distinguished function in the 1993 book by Olver [5]. The history of the term Casimir can be found in the notes to Chapter 6 of that book.

6. ^In the framework of PDEs, the fine distinction made between “Hamiltonian” and “Poisson” dynamical systems of ODEs does not seem too popular. In the 2006 book by Harier et al. [14, Section VII.2], a “Poisson” dynamical system is defined as the dynamical system obtained when (in our notation) the structure matrix J is a non-constant matrix which depends on the unknown u, J = J(u). We were perplexed by the explicit absence of the case in which J is constant but not necessarily invertible. However, later we understood that the “constant” case was trivially contained in the more interesting “non-constant” one.

7. ^Everything done for the scalar wave equation can be extended to the equations of elastodynamics very easily. See, for example, the Hamiltonian formulations in our recent work [2] and the six independent conservation laws in the 1976 work by Fletcher [20].

8. ^For a geometric interpretation of the symplecticity of an operator, see the 1992 review by Sanz-Serna [13], the 2006 book by Hairer et al. [14] and the 2010 book by Feng and Qin [8]. In the many discussions with our dynamical system colleagues, we always got the impression that symplecticity enhances the numerical schemes in much better ways than the simple conservation of energy can. We admit that we are more ready to believe their statements than to really understand their arguments. Perhaps because it is difficult to visualize orbits of infinite-dimensional dynamical systems. To be pragmatic, we restrict ourselves to the operational definition of symplecticity given by equality (b).

9. ^To us, the intuition of what is a Poisson integrator matrix can be extracted directly from Proposition 3.2: If the matrix A is of the form JH, with J anti-symmetric and H symmetric, then the matrix e^tA is a Poisson integrator matrix. We have the impression that the properties defining Poisson integrator matrices can be formulated precisely to match this remarkable relation.

References

1. SánchezMA, Ciuca C, Nguyen NC, Peraire J, Cockburn B. Symplectic Hamiltonian HDG methods for wave propagation phenomena. J Comput Phys. (2017) 350:951–73. doi: 10.1016/j.jcp.2017.09.010

CrossRef Full Text | Google Scholar

2. SánchezMA, Cockburn B, Nguyen NC, Peraire J. Symplectic Hamiltonian finite element methods for linear elastodynamics. Comput Methods Appl Mech Engrg. (2021) 381:113843. doi: 10.1016/j.cma.2021.113843

CrossRef Full Text | Google Scholar

3. SánchezMA, Du S, Cockburn B, Nguyen NC, Peraire J. Symplectic Hamiltonian finite element methods for electromagnetics. Comput Methods Appl Mech Engrg. (2022) 396:114969. doi: 10.1016/j.cma.2022.114969

CrossRef Full Text | Google Scholar

4. Fu G, Shu CW. Optimal energy-conserving discontinuous Galerkin methods for linear symmetric hyperbolic systems. J Comput Phys. (2019) 394:329–63. doi: 10.1016/j.jcp.2019.05.050

CrossRef Full Text | Google Scholar

5. Olver PJ. Applications of Lie Groups to Differential Equations. Vol. 107 of Graduate Texts in Mathematics. 2nd ed. New York, NY: Springer-Verlag (1993).

Google Scholar

6. Marsden JE, Ratiu TS. Introduction to Mechanics and Symmetry. Vol. 17 of Texts in Applied Mathematics. 2nd ed. New York, NY: Springer-Verlag (1999).

Google Scholar

7. Leimkuhler B, Reich S. Simulating Hamiltonian Dynamics. Vol. 14 of Cambridge Monographs on Applied and Computational Mathematics. Cambridge: Cambridge University Press (2004).

Google Scholar

8. Feng K, Qin MZ. Symplectic geometric algorithms for Hamiltonian systems. In: Zhejiang Science and Technology Publishing House, Hangzhou: Springer, Heidelberg; (2010).

Google Scholar

9. Cockburn B, Du S, Sánchez MA. Discontinuous Galerkin methods with time-operators in their numerical traces for time-dependent electromagnetics. Comput Methods Appl Math. (2022) 22:775–96. doi: 10.1515/cmam-2021-0215

CrossRef Full Text | Google Scholar

10. Ge Z, Feng K. On the Approximations of linear Hamiltonian systems. J Comput Math. (1988) 6:88–97.

Google Scholar

11. Feng K, Qin MZ. The symplectic methods for the computation of Hamiltonian equations. In: Numerical methods for partial differential equations (Shanghai, 1987). vol. 1297 of Lecture Notes in Math. Berlin: Springer (1987). p. 1–37.

Google Scholar

12. Feng K, Wu HM, Z M. Symplectic difference schemes for linear Hamiltonian systems. J Comput Math. (1990) 8:371–80.

Google Scholar

13. Sanz-Serna JM. Symplectic integrators for Hamiltonian problems: an overview. In: Acta Numerica 1992 Acta Numer. Cambridge: Cambridge Univ. Press (1992). p. 243–86.

Google Scholar

14. Hairer E, Lubich C, Wanner G. Geometric numerical integration : structure-preserving algorithms for ordinary differential equations. In: Springer Series in Computational Mathematics. Berlin; Heidelberg; New York, NY: Springer. (2006).

Google Scholar

15. Yee K. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. IEEE Trans Antennas Propagat. (1966) 14:302–7. doi: 10.1109/TAP.1966.1138693

CrossRef Full Text | Google Scholar

16. McLachlan R. Symplectic integration of Hamiltonian wave equations. Numer Math. (1993) 66:465–92. doi: 10.1007/BF01385708

CrossRef Full Text | Google Scholar

17. Groß M, Betsch P, Steinmann P. Conservation properties of a time FE method. Part IV: Higher order energy and momentum conserving schemes. Internat J Numer Methods Engrg. (2005) 63:1849–97. doi: 10.1002/nme.1339

CrossRef Full Text | Google Scholar

18. Xu Y, van der Vegt JJW, Bokhove O. Discontinuous Hamiltonian finite element method for linear hyperbolic systems. J Sci Comput. (2008) 35:241–65. doi: 10.1007/s10915-008-9191-y

CrossRef Full Text | Google Scholar

19. Kirby RC, Kieu TT. Symplectic-mixed finite element approximation of linear acoustic wave equations. Numer Math. (2015) 130:257–91. doi: 10.1007/s00211-014-0667-4

CrossRef Full Text | Google Scholar

20. Fletcher DC. Conservation laws in linear elastodynamics. Arch Rational Mech Anal. (1976) 60:329–53. doi: 10.1007/BF00248884

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Lipkin DM. Existence of a new conservation law in electromagnetic theory. J Math Phys. (1964) 5:695–700. doi: 10.1063/1.1704165

CrossRef Full Text | Google Scholar

22. Tang Y, Cohen AE. Optical chirality and its interaction with matter. Phys Rev Lett. (2010) 104:163901. doi: 10.1103/PhysRevLett.104.163901

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Yoshida H. Construction of higher order symplectic integrators. Phys Lett A. (1990) 150:262–8. doi: 10.1016/0375-9601(90)90092-3

CrossRef Full Text | Google Scholar

24. Bridges JT, Reich S. Numerical methods for Hamiltonian PDEs. J Phys A. (2006) 39:5287–320. doi: 10.1088/0305-4470/39/19/S02

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Taflove A. Application of the finite-difference time-domain method to sinusoidal steady-state electromagnetic-penetration problems. IEEE Trans Electromagn Compatibil. (1980) 3:191–202. doi: 10.1109/TEMC.1980.303879

CrossRef Full Text | Google Scholar

26. Sun Y, Tse PSP. Symplectic and multisymplectic numerical methods for Maxwell's equations. J Comput Phys. (2011) 230:2076–94. doi: 10.1016/j.jcp.2010.12.006

CrossRef Full Text | Google Scholar

27. Stern A, Tong Y, Desbrun M, Marsden JE. Geometric computational electrodynamics with variational integrators and discrete differential forms. In: Geometry, Mechanics, and Dynamics. vol. 73 of Fields Inst. Commun. New York, NY: Springer (2015). p. 437–75.

Google Scholar

28. Cockburn B, Fu Z, Hungria A, Ji L, Sánchez MA, Sayas FJ. Stormer-Numerov HDG methods for acoustic waves. J Sci Comput. (2018) 75:597–624. doi: 10.1007/s10915-017-0547-z

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Nguyen NC, Peraire J, Cockburn B. High-order implicit hybridizable discontinuous Galerkin methods for acoustics and elastodynamics. J Comput Phys. (2011) 230:3695–718. doi: 10.1016/j.jcp.2011.01.035

CrossRef Full Text | Google Scholar

30. Cockburn B, Quenneville-Bélair V. Uniform-in-time superconvergence of HDG methods for the acoustic wave equation. Math Comp. (2014) 83:65–85. doi: 10.1090/S0025-5718-2013-02743-3

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Hairer E, Wanner G. Solving ordinary differential equations. II, volume 14 of Springer Series in Computational Mathematics. 2nd ed. Berlin: Springer-Verlag (1996).

Google Scholar

32. Walter AS. Nonlinear invariant wave equations. In: Invariant Wave Equations (Proc. Ettore Majorana International School of Mathematical Physics Erice, 1977), Lecture Notes in Physics, Vol. 73. Berlin: Springer (1978). p. 197–249.

Google Scholar

33. Anco SC, Pohjanpelto J. Classification of local conservation laws of Maxwell's equations. Acta Appl. Math. (2001) 69:285–327. doi: 10.1023/A:1014263903283

CrossRef Full Text | Google Scholar

Appendix

1. Symplectic partitioned Runge-Kutta methods

Here, we provide a proof of Proposition 3.8 which characterizes the Symplectic Partitioned Runge-Kutta (PRK) methods for ODEs. We use the notation of Section 3.4.

We must show that $\Theta :={\text{u}}^{1,\top }{(\text{E}}_{\Delta t}\text{J}{\text{E}}_{\Delta t}^{\top }-\text{J}){\text{u}}^2$ is zero for any u¹ and u² in ℝ^d. By Theorem 3.3, ${\text{E}}_{\Delta t}^{\top }\text{v}=\text{v}$ for any element v of the kernel of J. Thus, we can take u¹ and u² in the range of J. On that subspace, J is invertible and it is enough to prove that $\Theta :={\text{u}}^{1,\top }{(\text{E}}_{\Delta t}^{\top }{\text{J}}^{-1}{\text{E}}_{\Delta t}^{\top }-{\text{J}}^{-1}{)\text{u}}^2$ is zero.

Inserting the definition of the PRK numerical method E_Δt, we obtain

$\begin{array}{l}\Theta =\Delta t{\displaystyle \sum _{j=1}^s{\text{u}}^{1,\top }}{\text{J}}^{-1}\text{A}{B}_j{\text{U}}_j^2+\Delta t{\displaystyle \sum _{i=1}^s{\text{U}}_i^{1,\top }}{B}_i{\text{A}}^{\top }{\text{J}}^{-1}{\text{u}}^2\\ +{(\Delta t)}^2{\displaystyle \sum _{i,j=1}^s{\text{U}}_i^{1,\top }}{B}_i{\text{A}}^{\top }{\text{J}}^{-1}\text{A}{B}_j{\text{U}}_j^2\\ =\Delta t{\displaystyle \sum _{j=1}^s{\text{u}}^{1,\top }}\text{H}{B}_j{\text{U}}_j^2+\Delta t{\displaystyle \sum _{i=1}^s{\text{U}}_i^{1,\top }}{B}_i\text{H}{\text{u}}^2\\ -{(\Delta t)}^2{\displaystyle \sum _{i,j=1}^s{\text{U}}_i^{1,\top }}{B}_i\text{HJH}{B}_j{\text{U}}_j^2,\end{array}$

because A = JH. Since $\text{u}={\text{U}}_{\ell }-\Delta t\text{A}\sum _{m=1}^s{A}_{\ell m}{\text{U}}_m,$ we get that

$\begin{array}{l}\Theta =\Delta t{\displaystyle \sum _{j=1}^s{({\text{U}}_j^1-\Delta t\text{JH}{\displaystyle \sum _{i=1}^s{A}_{ji}}{\text{U}}_i^1)}^{\top }}\text{H}{B}_j{\text{U}}_j^2\\ +\Delta t{\displaystyle \sum _{i=1}^s{\text{U}}_i^{1,\top }}{B}_i\text{H}({\text{U}}_i^2-\Delta t\text{JH}{\displaystyle \sum _{j=1}^s{A}_{ij}}{\text{U}}_j^2)\\ -{(\Delta t)}^2{\displaystyle \sum _{i,j=1}^s{\text{U}}_i^{1,\top }}{B}_i{\text{HJ}}^{\top }\text{H}{B}_j{\text{U}}_j^2,\end{array}$

and so

$\begin{array}{ll}\Theta =& \Delta t{\displaystyle \sum _{\ell =1}^s}{\text{U}}_{\ell }^{1,\top }{\text{H}}_{\ell }{\text{u}}_{\ell }^2+{(\Delta t)}^2{\displaystyle \sum _{i,j=1}^s}{\text{u}}_i^{1,\top }{\text{S}}_{ij}{\text{u}}_j^2,\end{array}$

where

$\begin{array}{lll}{\text{H}}_{\ell }& =\text{H}{B}_{\ell }-B_{\ell }\text{H}\text{and }{\text{S}}_{ij}& =A_{ji}\text{H}\text{J}\text{H}{B}_j+B_i\text{H}\text{J}\text{H}{A}_{ij}-B_i\text{H}\text{J}\text{H}{B}_j.\end{array}$

Next, we incorporate the information of the Hamiltonian being separable:

$\begin{array}{l}\text{H}=\left[\begin{array}{cc}{\text{H}}_{pp}& 0\\ 0& {\text{H}}_{qq}\end{array}\right]\text{ and HJH}=\left[\begin{array}{cc}0& {\Lambda }_{pq}\\ -{\Lambda }_{pq}^{\top }& 0\end{array}\right]\text{ where }\\ {\Lambda }_{pq}={\text{H}}_{pp}{\text{J}}_{pq}{\text{H}}_{qq},\end{array}$

and, with ${}^{\top }{B}_{\ell }:{=}^{\top }\left[\begin{array}{cc}{\overline{b}}_{\ell }\text{I}{\text{d}}_{pp}& 0\\ 0& b_{\ell }\text{I}{\text{d}}_{qq}\end{array}\right]=\left[\begin{array}{cc}{b}_{\ell }\text{I}{\text{d}}_{pp}& 0\\ 0& {\overline{b}}_{\ell }\text{I}{\text{d}}_{qq}\end{array}\right]$, conclude that

$\begin{array}{l}{\text{H}}_{\ell }=\text{H}(B_{\ell }-B_{\ell })=\text{0 and}\\ {\text{S}}_{ij}n=\text{ HJH}{(}^{\top }{A}_{ji}{B}_j{+}^{\top }{B}_i{A}_{ij}{-}^{\top }{B}_i{B}_j)\\ =\text{0},\end{array}$

by hypothesis. This completes the proof of Proposition 3.8.