Multisymplecticity of hybridizable discontinuous Galerkin methods

TL;DR

This paper establishes necessary and sufficient conditions for hybridizable discontinuous Galerkin methods to satisfy multisymplectic conservation laws, enabling structure-preserving discretizations of Hamiltonian PDEs with high-order accuracy.

Contribution

It provides a comprehensive characterization of multisymplecticity conditions for HDG methods, including common finite element variants, extending structure-preserving discretizations to high-order methods.

Findings

HDG methods satisfy multisymplectic conservation laws under specific conditions

Common finite element methods like mixed and nonconforming are multisymplectic when hybridized

Multisymplecticity holds for high-order methods on unstructured meshes

Abstract

In this paper, we prove necessary and sufficient conditions for a hybridizable discontinuous Galerkin (HDG) method to satisfy a multisymplectic conservation law, when applied to a canonical Hamiltonian system of partial differential equations. We show that these conditions are satisfied by the "hybridized" versions of several of the most commonly-used finite element methods, including mixed, nonconforming, and discontinuous Galerkin methods. (Interestingly, for the continuous Galerkin method in dimension greater than one, we show that multisymplecticity only holds in a weaker sense.) Consequently, these general-purpose finite element methods may be used for structure-preserving discretization (or semidiscretization) of canonical Hamiltonian systems of ODEs or PDEs. This establishes multisymplecticity for a large class of arbitrarily-high-order methods on unstructured meshes.

Equations

$${\partial }_{\mu }{u}^i={\varphi }_{\mu }^i(\cdot ,u,\sigma ),-{\partial }_{\mu }{\sigma }_i^{\mu }=f_i(\cdot ,u,\sigma ),$$

$${\partial }_{\mu }{u}^i=\frac{\partial H}{\partial {\sigma }_i^{\mu }},-{\partial }_{\mu }{\sigma }_i^{\mu }=\frac{\partial H}{\partial u^i}.$$

$${\dot{q}}^i=\frac{\partial H}{\partial p_i},-{\dot{p}}_i=\frac{\partial H}{\partial q^i}.$$

$$H(x,u,\sigma )=\frac{1}{2}{a}_{\mu \nu }(x){\sigma }^{\mu }{\sigma }^{\nu }+F(x,u),$$

$${\partial }_{\mu }u=a_{\mu \nu }{\sigma }^{\nu },-{\partial }_{\mu }{\sigma }^{\mu }=\frac{\partial F}{\partial u}.$$

$$-{\partial }_{\mu }{a}^{\mu \nu }{\partial }_{\nu }u=\frac{\partial F}{\partial u},$$

$$\mathrm{grad}u=a^{-1}\sigma ,-\mathrm{div}\sigma =\frac{\partial F}{\partial u},$$

$$-\mathrm{div}(a\mathrm{grad}u)=\frac{\partial F}{\partial u},$$

$$-\mathrm{div}(a\mathrm{grad}u)+cu=-{\partial }_{\mu }{a}^{\mu \nu }{\partial }_{\nu }u+cu=f.$$

$${\partial }_{\mu }{v}^i$$

$$-{\partial }_{\mu }{\tau }_i^{\mu }$$

$${\partial }_{\mu }{\omega }^{\mu }=0,$$

$${\partial }_{\mu }{\omega }^{\mu }$$

$$=\mathrm{d}{\varphi }_{\mu }^i\wedge \mathrm{d}{\sigma }_i^{\mu }-\mathrm{d}{u}^i\wedge \mathrm{d}{f}_i$$

$$=\mathrm{d}({\varphi }_{\mu }^i\mathrm{d}{\sigma }_i^{\mu }+f_i\mathrm{d}{u}^i),$$

$$\partial_{\mu}\bigl{(}\mathrm{d}u^{i}(v,\tau)\bigr{)}=\partial_{\mu}v^{i}=\mathrm{d}\phi^{i}_{\mu}(v,\tau)=\mathrm{d}(\partial_{\mu}u)(v,\tau).$$

$${\varphi }_{\mu }^i\mathrm{d}{\sigma }_i^{\mu }+f_i\mathrm{d}{u}^i=\frac{\partial H}{\partial {\sigma }_i^{\mu }}\mathrm{d}{\sigma }_i^{\mu }+\frac{\partial H}{\partial u^i}\mathrm{d}{u}^i=\mathrm{d}H,$$

$${\partial }_{\mu }u=a_{\mu \nu }{\sigma }^{\nu },-{\partial }_{\mu }{\sigma }^{\mu }=\frac{\partial F}{\partial u},$$

$${\partial }_{\mu }{\omega }^{\mu }$$

$$=\mathrm{d}({\partial }_{\mu }u)\wedge \mathrm{d}{\sigma }^{\mu }+\mathrm{d}u\wedge \mathrm{d}({\partial }_{\mu }{\sigma }^{\mu })$$

$$\displaystyle=a_{\mu\nu}\,\mathrm{d}\sigma^{\nu}\wedge\mathrm{d}\sigma^{\mu}+\mathrm{d}u\wedge\biggl{(}-\frac{\partial^{2}F}{\partial u^{2}}\,\mathrm{d}u-\frac{\partial^{2}F}{\partial\sigma\partial u}\,\mathrm{d}\sigma\biggr{)}.$$

$${\int }_K{\partial }_{\mu }{\omega }^{\mu }{\mathrm{d}}^{m}x={\int }_{\partial K}{\omega }^{\mu }{\mathrm{d}}^{m-1}{x}_{\mu },$$

$${\int }_{\partial K}{\omega }^{\mu }{\mathrm{d}}^{m-1}{x}_{\mu }=0,\forall K⋐U.$$

$${\omega }^{\mu }=\mathrm{d}{u}^i\wedge \mathrm{d}{\sigma }_i^{\mu }=\mathrm{d}{u}^i\otimes \mathrm{d}{\sigma }_i^{\mu }-\mathrm{d}{\sigma }_i^{\mu }\otimes \mathrm{d}{u}^i,$$

$${\int }_{\partial K}(\mathrm{d}{u}^i\otimes \mathrm{d}{\sigma }_i^{\mu }){\mathrm{d}}^{m-1}{x}_{\mu }={\int }_{\partial K}(\mathrm{d}{\sigma }_i^{\mu }\otimes \mathrm{d}{u}^i){\mathrm{d}}^{m-1}{x}_{\mu }.$$

$${\partial }_{\mu }v=a_{\mu \nu }{\tau }^{\nu },-{\partial }_{\mu }{\tau }^{\mu }=\frac{{\partial }^{2}F}{\partial u^2}v.$$

$${\int }_{\partial K}v{\tau }^{\prime \mu }{\mathrm{d}}^{m-1}{x}_{\mu }$$

$$\displaystyle=\int_{K}\biggl{(}a_{\mu\nu}\tau^{\nu}\tau^{\prime\mu}-v\frac{\partial^{2}F}{\partial u^{2}}v^{\prime}\biggr{)}\,\mathrm{d}^{m}x.$$

$$\int_{\partial K}v^{\prime}\tau^{\mu}\,\mathrm{d}^{m-1}x_{\mu}=\int_{K}\biggl{(}a_{\mu\nu}\tau^{\prime\nu}\tau^{\mu}-v^{\prime}\frac{\partial^{2}F}{\partial u^{2}}v\biggr{)}\,\mathrm{d}^{m}x.$$

$${\int }_{\partial K}v{\tau }^{\prime \mu }{\mathrm{d}}^{m-1}{x}_{\mu }={\int }_{\partial K}{v}^{\prime }{\tau }^{\mu }{\mathrm{d}}^{m-1}{x}_{\mu }.$$

$${\int }_{\partial K}v\mathrm{grad}{v}^{\prime }\cdot \mathbf{n}={\int }_{\partial K}{v}^{\prime }\mathrm{grad}v\cdot \mathbf{n}$$

$${\partial }_{\mu }(v^i{\tau }_i^{\prime \mu })={\partial }_{\mu }(v^{\prime i}{\tau }_i^{\mu }),$$

$${\int }_{\partial K}{v}^i{\tau }_i^{\prime \mu }{\mathrm{d}}^{m-1}{x}_{\mu }={\int }_{\partial K}{v}^{\prime i}{\tau }_i^{\mu }{\mathrm{d}}^{m-1}{x}_{\mu },$$

Full text

1. Introduction

1.1. Motivation and background

Hamiltonian systems of ordinary differential equations (ODEs) and partial differential equations (PDEs) are ubiquitous in applications, especially in the modeling of physical systems.

One essential property of a Hamiltonian ODE is that its time flow is symplectic: that is, it conserves a closed, nondegenerate $2$-form on phase space. This has motivated the development of symplectic integrators: one-step numerical integrators that, when applied to Hamiltonian ODEs, are also symplectic maps. It turns out that these methods have several numerical advantages that result from preserving the symplectic structure. Furthermore, most of these methods (such as symplectic partitioned Runge–Kutta methods) may be applied to general systems of ODEs, whether or not the user is aware of any Hamiltonian/symplectic structure—but if such a structure is present, then the methods will automatically preserve them. For a comprehensive survey of structure-preserving numerical integrators, including symplectic integrators, see Hairer et al. [22].

Similarly, Hamiltonian PDEs satisfy a multisymplectic conservation law. Since symplecticity is a desirable property for numerical integration of canonical Hamiltonian ODEs, it is natural to seek numerical methods for canonical Hamiltonian PDEs whose solutions satisfy the multisymplectic conservation law, in an appropriate sense. There has been important work on multisymplectic methods over the past two decades, particularly by Marsden and collaborators [27, 29, 28, 25] and Reich and collaborators [37, 38, 10, 35, 16]. However, most of these methods have consisted either of tensor products of symplectic Runge–Kutta-type methods on rectangular grids or of relatively low-order, finite-difference-type methods on unstructured meshes. For the variational integrators of Marsden et al., one must also know the (Lagrangian) geometric structure of the PDE, in advance, in order to devise the method.

The impact of multisymplecticity on solutions to PDEs, and on their discretizations, is not fully understood. However, multisymplecticity is known to be necessary to preserve traveling waves of hyperbolic equations [30], and compact multisymplectic methods can preserve dispersion relations much better than non-multisymplectic methods or noncompact finite difference methods [5, 34, 33]. For boundary value problems, multisymplecticity restricts the types of bifurcations that can occur [31, 32]. Because it is a local property, multisymplecticity is a strictly stronger property than the symplecticity obtained by integrating over space.

There has been some previous work on the application of finite element methods to certain problems on structured (especially rectangular) meshes. Guo et al. [21] considered the 2-D nonlinear Poisson equation on a regular rectangular grid, meshed with biased triangles, using the continuous Galerkin method with linear shape functions, and they showed that the degrees of freedom satisfied a multisymplectic finite-difference scheme. Zhen et al. [44] did the same for first-order rectangular Lagrange elements on a regular grid. Chen [11] used first- and second-order rectangular elements to derive Lagrangian variational methods (in the sense of Marsden et al. [27]) and applied these to the sine-Gordon equation on a regular grid, while pointing out that higher-order elements could also be used in principle. In all of these examples, however, finite elements were really only used as a tool to construct a finite-difference stencil on a regular, 2-D rectangular grid.

As finite element methods are traditionally formulated, there is a serious conceptual obstacle to discussing multisymplecticity. Namely, many classical finite element methods are posed on spaces of global functions, making it difficult to make sense of local properties like the multisymplectic conservation law. Indeed, the interpretation of multisymplecticity is much more straightforward for finite difference or finite volume methods with local stencils, or for tensor products of 1-D integrators on a rectangular grid, which may explain why the previous work has been focused on such methods.

Hybrid finite element methods provide a way around this obstacle, since they consist of local problems coupled through their boundary traces, where the boundary traces are allowed to be independent variables. (Oftentimes, these boundary traces are interpreted as “Lagrange multipliers” enforcing weak continuity between local regions.) While hybrid methods have a long history (see the comprehensive work by Brezzi and Fortin [7]), the recent work of Cockburn et al. [13] has shown that a wide variety of finite element methods—including those not previously thought of as hybrid methods—may be “hybridized” within a unified framework. Such methods are called hybridizable discontinuous Galerkin (HDG) methods, and they include not only several classical mixed and hybrid methods, but also hybridized versions of continuous and discontinuous Galerkin methods, nonconforming methods, and others. This framework provides precisely the local structure needed to examine multisymplecticity of the finite element methods in this class.

1.2. Organization of the paper

The paper is organized as follows:

• •

In Section 2, we review systems of PDEs in a particular canonical form. This form includes the de Donder–Weyl equations for a Hamiltonian and many elliptic and hyperbolic variational PDEs. We recall the multisymplectic conservation law for classical (i.e., smooth) solutions and illustrate how this manifests concretely for a class of semilinear elliptic PDEs in mixed form. We also discuss the relationship between multisymplecticity of solutions and reciprocity principles in physical systems.

• •

In Section 3, we develop a hybrid “flux formulation” for these systems of PDEs. As in Cockburn et al. [13], this yields a collection of weak problems on non-overlapping subdomains, coupled only through approximate traces on their shared boundaries. Our framework includes not only the linear second-order elliptic PDEs considered by Cockburn et al. [13], but also a more general class of nonlinear systems of PDEs, including canonical Hamiltonian PDEs.

Within this formulation, we establish criteria for solutions to satisfy weak and strong versions of the multisymplectic conservation law; the distinction is shown to be related to weak and strong conservativity of numerical fluxes (cf. Arnold et al. [4], Cockburn et al. [13]). In addition to the subsequent numerical applications, we also use a domain-decomposition argument to write the weak problem (in the sense of distributions) in this flux formulation, thereby establishing multisymplecticity for weak solutions to Hamiltonian PDEs.

• •

In Section 4, we examine several particular classes of HDG methods, including the hybridized Raviart–Thomas (RT-H), Brezzi–Douglas–Marini (BDM-H), local discontinuous Galerkin (LDG-H), continuous Galerkin (CG-H), nonconforming (NC-H), and interior penalty (IP-H) methods. These methods are posed in the framework of Section 3, and their multisymplecticity is then examined. Each of these methods, except for CG-H, is proved to be strongly multisymplectic. A counterexample shows that CG-H is only weakly multisymplectic, resulting from the fact that it is only weakly conservative.

2. Canonical and multisymplectic systems of PDEs

2.1. Canonical systems of PDEs

Given a domain $U\subset\mathbb{R}^{m}$, consider a system of first-order PDEs having the form

[TABLE]

where $\mu=1,\ldots,m$ and $i=1,\ldots,n$. Here, $u=u^{i}(x)$ and $\sigma=\sigma_{i}^{\mu}(x)$ are unknown functions on $U$, while $\phi=\phi^{i}_{\mu}(x,u,\sigma)$ and $f=f_{i}(x,u,\sigma)$ are given functions on $U\times\mathbb{R}^{n}\times\mathbb{R}^{mn}$. We abbreviate $\partial_{\mu}\coloneqq\partial/\partial x^{\mu}$ and adopt the Einstein index convention of summing over repeated indices—so, for instance, the expression $\partial_{\mu}\sigma_{i}^{\mu}$ in (1) has an implied sum over $\mu$ and may therefore be interpreted as the divergence of $\sigma_{i}$.

Among these is the important class of (canonical) Hamiltonian systems,

[TABLE]

where $H=H(x,u,\sigma)$ is a function called the Hamiltonian. In the special case $m=1$, the resulting system of ODEs yields Hamilton’s equations of classical mechanics, which are usually written as

[TABLE]

The equations (2) are called the de Donder–Weyl equations (de Donder [15], Weyl [43]). These canonical systems are an important special case of a more general class of Hamiltonian systems of PDEs, cf. Bridges [8, 9].

Throughout this section, we assume that all of the functions above are smooth. Later, in Section 3, we will relax this assumption in order to introduce a weak formulation of (1).

Example \theexample (semilinear elliptic PDE).

Let $n=1$, so that $u=u(x)$ is a scalar field and $\sigma=\sigma^{\mu}(x)$ a vector field on $U\subset\mathbb{R}^{m}$. Consider

[TABLE]

where $a=a^{\mu\nu}(x)$ is symmetric and positive-definite with matrix inverse a_{\mu\nu}(x)\coloneqq\bigl{(}a^{\mu\nu}(x)\bigr{)}^{-1} at each $x\in U$, and where $F$ is arbitrary. Then the de Donder–Weyl equations for this Hamiltonian are

[TABLE]

From the first of these equations, we have $\sigma^{\mu}=a^{\mu\nu}\partial_{\nu}u$, so substituting this into the second yields

[TABLE]

which is a second-order semilinear elliptic PDE in divergence form. Note that we could also have written the de Donder–Weyl equations in the equivalent, coordinate-free form,

[TABLE]

so the substitution $\sigma=a\operatorname{grad}u$ yields

[TABLE]

which is an equivalent expression for the second-order PDE above.

An important special case is when $F(x,u)=f(x)u-\frac{1}{2}c(x)u^{2}$ for given $f$ and $c$ on $U$. In this case, we obtain a linear second-order elliptic PDE,

[TABLE]

In particular, if $a^{\mu\nu}\equiv\delta^{\mu\nu}$ (where $\delta$ is the Kronecker delta, i.e., $a$ is the identity matrix) and $c\equiv 0$, then this simply becomes Poisson’s equation $-\Delta u=f$ on $U$.

We will regularly return to this example throughout the paper.

2.2. The multisymplectic conservation law

Define the collection of canonical 2-forms $\omega^{\mu}\coloneqq\mathrm{d}u^{i}\wedge\mathrm{d}\sigma_{i}^{\mu}$ on $\mathbb{R}^{n}\times\mathbb{R}^{mn}$, for $\mu=1,\ldots,m$.

Notation \thenotation.

Unless otherwise stated, differential forms and exterior differential operators (such as $\mathrm{d}$, $\wedge$, etc.) are assumed to be on $\mathbb{R}^{n}\times\mathbb{R}^{mn}$, where $x\in U$ (if it appears) is fixed. Differentiation with respect to $x$ will always be denoted using the previously-defined $\partial_{\mu}$ notation.

Definition \thedefinition.

Let $(u,\sigma)$ be a solution to (1). A (first) variation of $(u,\sigma)$ is a solution $(v,\tau)$ of the linearized problem

[TABLE]

The system (1) is multisymplectic if \partial_{\mu}\Bigl{(}\omega^{\mu}\bigl{(}(v,\tau),(v^{\prime},\tau^{\prime})\bigr{)}\Bigr{)}=0 for any pair of variations $(v,\tau)$ and $(v^{\prime},\tau^{\prime})$. This is abbreviated by

[TABLE]

where it is understood that $\omega^{\mu}$ is evaluated on variations of solutions to (1). The equation (4) is called the multisymplectic conservation law.

Lemma \thelemma.

The system (1) is multisymplectic if and only if the $1$-form $\phi^{i}_{\mu}\,\mathrm{d}\sigma_{i}^{\mu}+f_{i}\,\mathrm{d}u^{i}$ on $\mathbb{R}^{n}\times\mathbb{R}^{mn}$ is closed for each $x\in U$.

Proof.

Using (1), we calculate

[TABLE]

so the first expression vanishes if and only if the last expression vanishes. ∎

Remark \theremark.

Certain steps in this calculation, such as $\partial_{\mu}(\mathrm{d}u^{i})=\mathrm{d}(\partial_{\mu}u^{i})$, are seen to be valid by evaluating both sides on arbitrary variations of $(u,\sigma)$:

[TABLE]

Similar calculations with differential forms will appear throughout this paper, where they are interpreted as holding for arbitrary variations of solutions.

Corollary \thecorollary.

Every Hamiltonian system is multisymplectic.

Proof.

From (2), we have

[TABLE]

which is exact and therefore closed. ∎

Corollary \thecorollary.

Every multisymplectic system is Hamiltonian.

Proof.

Since $\mathbb{R}^{n}\times\mathbb{R}^{mn}$ is simply connected, its first de Rham cohomology is trivial. Hence, the closed $1$-form $\phi^{i}_{\mu}\,\mathrm{d}\sigma_{i}^{\mu}+f_{i}\,\mathrm{d}u^{i}$ is exact, i.e., it equals $\mathrm{d}H$ for some Hamiltonian $H$. More precisely, for each fixed $x\in U$, the $1$-form equals $\mathrm{d}H_{x}$ for some function $H_{x}$ on $\mathbb{R}^{n}\times\mathbb{R}^{mn}$, and these may be combined into a single Hamiltonian $H(x,u,\sigma)\coloneqq H_{x}(u,\sigma)$. ∎

Remark \theremark.

Section 2.2 depends entirely on the fact that $\mathbb{R}^{n}\times\mathbb{R}^{mn}$ has trivial first de Rham cohomology. However, it is possible to define canonical Hamiltonian systems on more general spaces—in particular, on the dual jet bundle of some fiber bundle over $U$ (cf. Gotay [17] and references therein). In this setting, the argument of Section 2.2 holds only if the fibers of this bundle have trivial first de Rham cohomology. However, a weaker statement—that every multisymplectic system is locally Hamiltonian—still holds, by Poincaré’s lemma. By contrast, Section 2.2 remains true even in this more general setting.

Example \theexample (semilinear elliptic PDE, continued).

Let us see how the multisymplectic conservation law manifests in the class of semilinear elliptic PDEs we encountered in Section 2.1. For the system

[TABLE]

we calculate

[TABLE]

The first term vanishes by the symmetry of $a$ and the antisymmetry of $\wedge$, while the remaining terms vanish since $\mathrm{d}u\wedge\mathrm{d}u=0$ (again, the antisymmetry of $\wedge$) and since $F=F(x,u)$ does not depend on $\sigma$.

2.3. Integral form of the multisymplectic conservation law

Given an arbitrary subdomain $K\Subset U$, the divergence theorem implies that

[TABLE]

where $\mathrm{d}^{m}x\coloneqq\mathrm{d}x^{1}\wedge\cdots\wedge\mathrm{d}x^{m}$ is the standard Euclidean volume form on $U$ and $\mathrm{d}^{m-1}x_{\mu}\coloneqq\iota_{\partial/\partial x^{\mu}}(\mathrm{d}^{m}x)$, where $\iota$ is the interior product (or contraction). Therefore, an equivalent formulation of the multisymplectic conservation law (4) is that

[TABLE]

We call this the integral form of the multisymplectic conservation law. As with (4), this is interpreted as holding when $\omega^{\mu}$ is evaluated on arbitrary variations of a solution to (1).

Note that, by the definition of the $\wedge$ product,

[TABLE]

so (5) may also be written as

[TABLE]

Hence, the multisymplectic conservation law may be interpreted as a symmetry condition on the Poincaré–Steklov operator mapping Dirichlet boundary conditions for $u$ to the corresponding boundary conditions for $\sigma$. See Agoshkov [3], where the symmetry of the Poincaré–Steklov operator is discussed in the context of domain decomposition methods for linear elliptic problems. Similarly, Belishev and Sharafutdinov [6] establish the symmetry of a Poincaré–Steklov operator for harmonic differential forms on a manifold with boundary.

Example \theexample (semilinear elliptic PDE, continued).

Let us revisit the class of semilinear elliptic PDEs we encountered in Section 2.1 and Section 2.2. Let $(u,\sigma)$ be a solution, and consider the linearized problem

[TABLE]

Here, $\partial^{2}F/\partial u^{2}$ is evaluated at \bigl{(}x,u(x)\bigr{)} and hence is a function of $x$ alone. If $(v,\tau)$ and $(v^{\prime},\tau^{\prime})$ are two arbitrary solutions to this problem, then

[TABLE]

By a similar calculation, switching $(v,\tau)$ with $(v^{\prime},\tau^{\prime})$,

[TABLE]

However, since $a$ is symmetric, these integrals are identical, and we conclude

[TABLE]

That this equality holds for every $(v,\tau)$ and $(v^{\prime},\tau^{\prime})$ is precisely the statement (6) of the multisymplectic conservation law.

In the special case $a^{\mu\nu}\equiv\delta^{\mu\nu}$, we have $\tau=\operatorname{grad}v$ and $\tau^{\prime}=\operatorname{grad}v^{\prime}$, so this can be written as

[TABLE]

where $\mathbf{n}$ denotes the outer unit normal to $\partial K$. Hence, in this case, the multisymplectic conservation law expresses the symmetry of the Dirichlet–Neumann operator $v\rvert_{\partial K}\mapsto\operatorname{grad}v\cdot\mathbf{n}\rvert_{\partial K}$ for the linearized problem.

2.4. Multisymplecticity and reciprocity

In many physical systems, the multisymplectic conservation law is closely tied to so-called reciprocity phenomena, such as Green’s reciprocity in electrostatics and Betti reciprocity in elasticity. (See, for example, Abraham and Marsden [1, Section 5.3], Marsden and Hughes [26, Section 5.6], and Lew et al. [25].) These reciprocity phenomena are also exploited, numerically, in formulations of the boundary element method (cf. Partridge et al. [36]). We now briefly discuss the relationship between multisymplecticity and reciprocity, using the language we have developed throughout this section.

Let $(u,\sigma)$ be a solution to (1). The multisymplectic conservation law (4) is just the statement that

[TABLE]

where $(v,\tau)$ and $(v^{\prime},\tau^{\prime})$ are arbitrary variations of $(u,\sigma)$, i.e., solutions to the linearized problem (3). Integrating both sides over $K\Subset U$ and applying the divergence theorem gives

[TABLE]

which is the statement of (5) and (6).

We now generalize the above to the case where the variations $(v,\tau)$ and $(v^{\prime},\tau^{\prime})$ each solve a perturbed version of the linearized problem.

Definition \thedefinition.

Given a solution $(u,\sigma)$ to (1), we say that $(v,\tau)$ solves the linearized problem with incremental sources $\psi$ and $g$ if

[TABLE]

where $\psi=\psi^{i}_{\mu}(x,v,\tau)$ and $g=g_{i}(x,v,\tau)$ are given functions.

Let $(v,\tau)$ be a solution to the linearized problem with incremental sources $\psi$ and $g$, and let $(v^{\prime},\tau^{\prime})$ be a solution to the linearized problem with incremental sources $\psi^{\prime}$ and $g^{\prime}$. By the Leibniz rule,

[TABLE]

which we subtract and rearrange to obtain

[TABLE]

Assuming that (1) is multisymplectic, the $\phi$ and $f$ terms cancel when we substitute (7), leaving the equation

[TABLE]

Integrating over $K\Subset U$ and applying the divergence theorem gives

[TABLE]

The equations (8) and (9) are the differential and integral forms, respectively, of the reciprocity law for a multisymplectic system of PDEs. They may be interpreted as describing a symmetric (or “reciprocal”) relationship between the perturbation of the system by incremental sources and the incremental response of the system to such perturbations. In the special case where the incremental sources vanish, we recover the multisymplectic conservation law.

Example \theexample (semilinear elliptic PDE, continued).

Let us once again examine the semilinear elliptic PDEs considered in Section 2.1, Section 2.2, and Section 2.3. The linearized problem with incremental sources is

[TABLE]

To see how reciprocity arises, we compute

[TABLE]

Subtracting, the terms involving $a$ and $F$ cancel by symmetry, yielding

[TABLE]

which is precisely the statement (8). Integrating this over $K\Subset U$ gives

[TABLE]

which is the statement (9).

As an important special case, which arises in the primal (or Lagrangian) formulation of this system, suppose that $\psi=0$ and that $g=g(\cdot,v)$, so that $\tau=a\operatorname{grad}v$ and $v$ solves the second-order equation

[TABLE]

If the corresponding properties hold for $(v^{\prime},\tau^{\prime})$, then we can write the reciprocity law as

[TABLE]

whose integral form on $K\Subset U$ is

[TABLE]

As a final specialization, let $a^{\mu\nu}\equiv\delta^{\mu\nu}$ and $F(x,u)=f(x)u$, so that $u$ satisfies Poisson’s equation, $-\Delta u=f$, and $v$ and $v^{\prime}$ satisfy

[TABLE]

Then the reciprocity law, in differential form, is

[TABLE]

while the integral form on $K\Subset U$ is

[TABLE]

These last two expressions are two of Green’s identities from vector calculus. If $v$ and $v^{\prime}$ are interpreted as scalar potentials for the electrostatic fields $\tau$ and $\tau^{\prime}$, respectively, then this corresponds to Green’s reciprocity.

3. The flux formulation and multisymplecticity

3.1. Domain decomposition and the flux formulation

In this section, we introduce a weak formulation of the problem (1), called the flux formulation. This decomposes the problem on $U$ into a collection of local solvers for $(u,\sigma)$, coupled through the approximate boundary traces $(\widehat{u},\widehat{\sigma})$. This forms the foundation of the HDG framework of Cockburn et al. [13] and is closely related to the unified DG framework of Arnold et al. [4].

We mention that our presentation of the flux formulation, and of HDG methods, differs from that in Cockburn et al. [13] in a few ways. In particular, Cockburn et al. focus on linear elliptic problems, which allows them to make substantial use of the solution theory for such problems, including well-posedness of the local and global solvers. By contrast, we are interested in obtaining multisymplecticity criteria for the much more general class of systems (1), without assuming anything about the properties of solutions, even their existence and/or uniqueness.

In this section, the function spaces appearing in the flux formulation may be either infinite-dimensional Hilbert spaces or finite-dimensional subspaces (e.g., polynomials up to some degree). This will allow us to prove general multisymplecticity results that apply both to the original, infinite-dimensional problem (Section 3.4) and to finite-element approximation via HDG methods (Section 4).

To begin, observe that if $(u,\sigma)$ is a smooth solution to (1) on $U$, and if $(v,\tau)$ are arbitrary smooth test functions, then

[TABLE]

If $\mathcal{T}_{h}$ is a partition of $U$ into non-overlapping domains $K\in\mathcal{T}_{h}$, then breaking each of the integrals above into a sum over $K\in\mathcal{T}_{h}$ and integrating by parts,

[TABLE]

In a typical finite-element application, $U$ will be polyhedral, and $\mathcal{T}_{h}$ will be a triangulation of $U$ into simplices $K\in\mathcal{T}_{h}$.

Following Cockburn et al. [13] (as well as Arnold et al. [4]), we now relax the regularity and inter-element continuity assumptions on $u,v$ and $\sigma,\tau$ in (10), and we replace the boundary traces of $u$ and $\sigma$ on $\partial K$ by approximate traces $\widehat{u}$ and $\widehat{\sigma}$. Specifically, let

[TABLE]

be specified local function spaces on each $K\in\mathcal{T}_{h}$, and define discontinuous global spaces on $U$ by

[TABLE]

Next, specify a space of approximate traces of functions in $V$,

[TABLE]

where $\mathcal{E}_{h}\coloneqq\bigcup_{K\in\mathcal{T}_{h}}\partial K$, along with the subspace

[TABLE]

of approximate traces vanishing on the domain boundary $\partial U$.

The final ingredient in the HDG framework is the numerical flux $\widehat{\sigma}$, which we define slightly differently to Cockburn et al. [13]. As mentioned in the introduction to this section, our treatment is equivalent to [13] for linear problems, but it extends more naturally to nonlinear problems, even without assuming existence and uniqueness of solutions. Let

[TABLE]

be some space of boundary fluxes on $\partial K$, and define the space of restricted traces \widehat{V}(\partial K)\coloneqq\bigl{\{}\widehat{v}\rvert_{\partial K}:\widehat{v}\in\widehat{V}\bigr{\}}.

Definition \thedefinition.

A local flux function on $K\in\mathcal{T}_{h}$ is a bounded linear map

[TABLE]

Denoting \widehat{\Sigma}\coloneqq\prod_{K\in\mathcal{T}_{h}}\widehat{\Sigma}(\partial K)\subset\prod_{K\in\mathcal{T}_{h}}\bigl{[}L^{2}(\partial K)\bigr{]}^{mn}, this extends naturally to a global flux function,

[TABLE]

Remark \theremark.

Elements of $\widehat{\Sigma}$ and \prod_{K\in\mathcal{T}_{h}}\bigl{[}L^{2}(\partial K)\bigr{]}^{mn} may be interpreted as functions that are double-valued on internal facets of $\mathcal{T}_{h}$ and single-valued on boundary facets in $\partial U$.

We now seek solutions $(u,\sigma,\widehat{u},\widehat{\sigma})\in V\times\Sigma\times\widehat{V}\times\widehat{\Sigma}$ satisfying

[TABLE]

the latter of which states that the normal component of $\widehat{\sigma}$ is single-valued (at least in a weak sense) on the internal facets of $\mathcal{T}_{h}$.

In the language of Cockburn et al. [13], we have “local solvers” (11a)–(11b) on each $K\in\mathcal{T}_{h}$, and these are coupled globally through the numerical flux $\widehat{\sigma}$ by the conservativity condition (11d). A notable distinction between our approach and that of Cockburn et al. [13] is that they assume $\widehat{\sigma}=\widehat{\sigma}(u,\sigma,\widehat{u})$ is a given function of $u$, $\sigma$, and $\widehat{u}$ on each $K\in\mathcal{T}_{h}$, whereas we define it through the flux functions $\Phi_{K}$ by adding (11c) to the flux formulation.

Definition \thedefinition.

The flux formulation of (1) on $\mathcal{T}_{h}$ is given by (11), along with choices of the global function space $\widehat{V}$ and, for each $K\in\mathcal{T}_{h}$, the local function spaces $V(K)$, $\Sigma(K)$, $\widehat{\Sigma}(\partial K)$ and the flux function $\Phi_{K}$. We call this a hybridizable discontinuous Galerkin (HDG) method whenever $\widehat{V}_{0}$, $V$, and $\Sigma$ (but not necssarily $\widehat{V}$ or $\widehat{\Sigma}$) are finite dimensional.

Remark \theremark.

Note that (11) does not impose any particular boundary conditions on $\widehat{u}\rvert_{\partial U}$ or $\widehat{\sigma}\rvert_{\partial U}$. Hence, (11) corresponds to the system of PDEs (1) rather than a particular boundary value problem associated to (1).

We remain agnostic about the choice of boundary conditions for two reasons. First, multisymplecticity is not a statement about a particular solution, but a statement about variations within a general family of solutions. If we pick out an isolated solution (e.g., by the imposition of boundary conditions), then the “family” of solutions becomes zero-dimensional, so any statement about variations is vacuous. Second, the class of PDEs (1) is quite general, including both elliptic and hyperbolic PDEs, among others, depending on $\phi$ and $f$. In the hyperbolic case, when $U$ is a spacetime region, we are not free to impose Dirichlet conditions on all of $\partial U$.

3.2. Local multisymplecticity criteria

In the context of smooth solutions to (1), where $\phi^{i}_{\mu}$ and $f_{i}$ are smooth functions on $U\times\mathbb{R}^{n}\times\mathbb{R}^{mn}$, Section 2.2 states that multisymplecticity holds if and only if the smooth $1$-form $\phi^{i}_{\mu}\,\mathrm{d}\sigma_{i}^{\mu}+f_{i}\,\mathrm{d}u^{i}$ is closed for all $x\in U$. For the flux formulation (11), we wish to relax these smoothness assumptions and express the multisymplecticity condition in terms of function spaces, rather than in a pointwise sense at each $x\in U$.

Observe that (11a)–(11b) still makes sense even if we only have $\phi^{i}_{\mu},f_{i}\in L^{2}(U)$ for $\mu=1,\ldots,m$ and $i=1,\ldots,n$. Therefore, rather than assuming that $\phi$ and $f$ are smooth, let us assume only that

[TABLE]

for each $K\in\mathcal{T}_{h}$, which may be extended naturally to

[TABLE]

These are generally nonlinear maps—and since multisymplecticity is a statement about first variations of solutions, let us assume also that these maps are at least $C^{1}$. When $V\times\Sigma$ is infinite-dimensional, we may interpret this as a variational derivative (either the Gâteaux or Fréchet derivative, which are equivalent for $C^{1}$ maps, cf. Abraham et al. [2, Corollary 2.10]); in the finite-dimensional case, this is just ordinary continuous differentiability. With $\phi$ and $f$ defined in this way, it follows that $\phi^{i}_{\mu}\,\mathrm{d}\sigma_{i}^{\mu}+f_{i}\,\mathrm{d}u^{i}$ is a $C^{1}$ differential $1$-form on $V\times\Sigma$, and we say that this $1$-form is closed if its exterior derivative vanishes as a $2$-form on $V\times\Sigma$.

Definition \thedefinition.

The flux formulation (11) is multisymplectic if solutions satisfy

[TABLE]

for all $K\in\mathcal{T}_{h}$, whenever the $C^{1}$ differential $1$-form $\phi^{i}_{\mu}\,\mathrm{d}\sigma_{i}^{\mu}+f_{i}\,\mathrm{d}u^{i}$ is closed.

Remark \theremark.

The condition (12) is essentially the integral form of the multisymplectic conservation law from Section 2.3, where the approximate traces $\widehat{u}$ and $\widehat{\sigma}$ are used instead of the actual traces of $u$ and $\sigma$. Note that (12) only needs to hold for $K\in\mathcal{T}_{h}$, not for arbitrary subdomains $K\Subset U$ as in (5).

We say that (12) is a local multisymplecticity condition because it is a statement only about the local solvers (11a)–(11b) and numerical flux condition (11c) for each $K\in\mathcal{T}_{h}$. We reserve the global question—whether the multisymplectic conservation law holds for arbitrary unions of elements of $\mathcal{T}_{h}$—for the next section, where the conservativity condition (11d) will also come into play.

To characterize the multisymplecticity of the flux formulation, we first prove a useful lemma, which relates the boundary integral in (12) to the “jumps” $\widehat{u}-u$ and $\widehat{\sigma}-\sigma$ between the approximate and actual traces on $\partial K$.

Lemma \thelemma.

If $\phi^{i}_{\mu}\,\mathrm{d}\sigma_{i}^{\mu}+f_{i}\,\mathrm{d}u^{i}$ is closed and $(u,\sigma,\widehat{u},\widehat{\sigma})\in V\times\Sigma\times\widehat{V}\times\widehat{\Sigma}$ satisfies (11a)–(11b) for $K\in\mathcal{T}_{h}$, then

[TABLE]

Consequently, the local multisymplecticity condition (12) holds if and only if

[TABLE]

Proof.

Since (11a) and (11b) hold for all $v\in V(K)$ and $\tau\in\Sigma(K)$, we have

[TABLE]

so taking exterior derivatives gives

[TABLE]

Subtracting the second equation from the first, the terms involving $\phi$ and $f$ vanish by the closedness assumption, so we are left with

[TABLE]

that is,

[TABLE]

Using this identity, we finally calculate

[TABLE]

which completes the proof. ∎

The equation (13) says that the multisymplecticity of the flux formulation depends entirely on the relationship among $u$, $\sigma$, $\widehat{u}$, and $\widehat{\sigma}$ on $\partial K$ for $K\in\mathcal{T}_{h}$. That is, it depends entirely on the choice of local flux functions $\Phi_{K}$.

Definition \thedefinition.

A local flux function $\Phi_{K}$ is multisymplectic if (13) holds whenever $(u,\sigma,\widehat{u},\widehat{\sigma})\in V\times\Sigma\times\widehat{V}\times\widehat{\Sigma}$ satisfies (11c).

We now prove multisymplecticity for two particular choices of $\Phi_{K}$. The first is used for the hybridized Raviart–Thomas (RT-H), Brezzi–Douglas–Marini (BDM-H), and local discontinuous Galerkin (LDG-H) methods; the second is used for the hybridized continuous Galerkin (CG-H) and nonconforming (NC-H) methods. These methods will be discussed further in Section 4.

Theorem 3.1**.**

Suppose that, for all $v\in V(K)$ and $\widehat{v}\in\widehat{V}$, there exists $\widehat{\tau}\in\widehat{\Sigma}(\partial K)$ such that $\widehat{\tau}_{i}^{\mu}\mathbf{n}_{\mu}=\delta_{ij}(\widehat{v}^{j}-v^{j})\rvert_{\partial K}$ for all $i=1,\ldots,n$. Then, for any $\lambda\in L^{\infty}(\partial K)$ (which is called a “penalty function”), the local flux function

[TABLE]

is multisymplectic.

Proof.

The flux condition (11c) says that

[TABLE]

for all $\widehat{\tau}\in\widehat{\Sigma}(\partial K)$. By assumption, for any $v\in V(K)$ and $\widehat{v}\in\widehat{V}$, we can choose $\widehat{\tau}\in\widehat{\Sigma}(\partial K)$ such that $\widehat{\tau}_{i}^{\mu}\mathbf{n}_{\mu}=\delta_{ij}(\widehat{v}^{j}-v^{j})\rvert_{\partial K}$ for all $i=1,\ldots,n$, and therefore

[TABLE]

Since $v$ and $\widehat{v}$ are arbitrary, this can be written as

[TABLE]

and taking the exterior derivative of both sides yields

[TABLE]

However, this vanishes by the symmetry of $\delta$ and the antisymmetry of $\wedge$, so the multisymplecticity condition (13) holds, as claimed. ∎

Remark \theremark.

More generally, we can replace $\lambda\delta_{ij}\in L^{\infty}(\partial K)$ above with penalty functions $\lambda_{ij}\in L^{\infty}(\partial K)$ such that $\lambda_{ij}=\lambda_{ji}$ for $i,j=1,\ldots,n$, and the argument above still holds. This same generalization applies to the penalty-based HDG methods we will encounter in Section 4.

Theorem 3.2**.**

Suppose that, for all $\tau\in\Sigma(K)$, there exists $\widehat{\tau}\in\widehat{\Sigma}(\partial K)$ such that $\widehat{\tau}_{i}^{\mu}\mathbf{n}_{\mu}=\tau_{i}^{\mu}\mathbf{n}_{\mu}\rvert_{\partial K}$ for $i=1,\ldots,n$. Then the local flux function

[TABLE]

is multisymplectic.

Proof.

The flux condition (11c) says that

[TABLE]

for all $\widehat{\tau}\in\widehat{\Sigma}(\partial K)$. From the assumption on normal traces of elements of $\Sigma(K)$, it follows that

[TABLE]

for any $\tau\in\Sigma(K)$ and $\widehat{\tau}\in\widehat{\Sigma}(\partial K)$. This can be written as

[TABLE]

and taking the exterior derivative yields (13). ∎

3.3. Global multisymplecticity criteria

Whenever (11) is multisymplectic, we may of course sum (12) over an arbitrary collection of elements $\mathcal{K}\subset\mathcal{T}_{h}$ to obtain the global statement

[TABLE]

or equivalently, by Section 3.2

[TABLE]

In the special case $\mathcal{K}=\mathcal{T}_{h}$, for a second-order linear elliptic PDE, (15) is precisely Equation 2.10 from Cockburn et al. [13], which establishes the symmetry of the bilinear form used to solve for $\widehat{u}$ once internal degrees of freedom have been eliminated (i.e., the Schur complement). Hence, (15) can be seen as a generalization of this symmetry condition to nonlinear multisymplectic systems and arbitrary $\mathcal{K}\subset\mathcal{T}_{h}$.

However, (14) is a rather weak global multisymplecticity condition, since it follows trivially from the local condition (12). We now define a stronger version of multisymplecticity, which is more analogous to the classical multisymplectic conservation law (5).

Definition \thedefinition.

The flux formulation (11) is strongly multisymplectic if solutions satisfy

[TABLE]

for all $\mathcal{K}\subset\mathcal{T}_{h}$, whenever the $C^{1}$ differential $1$-form $\phi^{i}_{\mu}\,\mathrm{d}\sigma_{i}^{\mu}+f_{i}\,\mathrm{d}u^{i}$ is closed.

Taking $\mathcal{K}=\{K\}$ immediately implies the local condition (12) for each $K\in\mathcal{T}_{h}$, so the stronger condition (16) indeed implies the weaker condition (14). It follows that (16) holds if and only if the terms of (14) cancel on internal facets. The property that conservation laws “add up” correctly over unions of elements is directly related to the conservativity condition (11d); indeed, this is the reason the term “conservative” is used to describe numerical fluxes with single-valued normal components. (See Equation 3.2 in Arnold et al. [4].)

Let $e=\partial K^{+}\cap\partial K^{-}$ be an internal facet, where $K^{\pm}\in\mathcal{T}_{h}$ are distinct. Recall from Section 3.1 that an element $\widehat{\tau}\in\widehat{\Sigma}$ is generally double-valued on $e$, since the $\widehat{\Sigma}(\partial K^{+})$ and $\widehat{\Sigma}(\partial K^{-})$ components need not agree. As is common in the discontinuous Galerkin literature (including Arnold et al. [4], Cockburn et al. [13]), we define the “normal jump” \llbracket\widehat{\tau}\rrbracket\rvert_{e}\in\bigl{[}L^{2}(e)\bigr{]}^{n} by

[TABLE]

Here, $e^{\pm}$ denotes that $e$ is oriented according to $\partial K^{\pm}$, and we take the corresponding component of $\widehat{\tau}$ and outer normal $\mathbf{n}$ for each term on the right-hand side. Denoting the set of internal facets of $\mathcal{T}_{h}$ by $\mathcal{E}_{h}^{\circ}\coloneqq\mathcal{E}_{h}\setminus\partial U$, we may sum directly over $e\in\mathcal{E}_{h}^{\circ}$ to define \llbracket\widehat{\tau}\rrbracket\in\bigl{[}L^{2}(\mathcal{E}_{h}^{\circ})\bigr{]}^{n}. Hence, the normal component of $\widehat{\tau}$ is single-valued on internal facets if and only if $\llbracket\widehat{\tau}\rrbracket=0$.

With this notation, the conservativity condition (11d) can be rewritten as

[TABLE]

If the extension by zero of $\llbracket\widehat{\sigma}\rrbracket$ to $\mathcal{E}_{h}$ is in $\widehat{V}_{0}$, then applying the conservativity condition with $\widehat{v}=\llbracket\widehat{\sigma}\rrbracket$ immediately implies $\llbracket\widehat{\sigma}\rrbracket=0$, and we say that $\widehat{\sigma}$ is strongly conservative. However, this is not always the case: for example, if $\widehat{\Sigma}(\partial K)$ contains discontinuous traces but $\widehat{V}(\partial K)$ contains only continuous traces, then in general $\llbracket\widehat{\sigma}\rrbracket\notin\widehat{V}_{0}$, so we cannot conclude that $\llbracket\widehat{\sigma}\rrbracket$ vanishes. In this case, we say that $\widehat{\sigma}$ is only weakly conservative. (This terminology is taken from Cockburn et al. [13].)

Theorem 3.3**.**

If $\widehat{\sigma}\in\widehat{\Sigma}$ satisfies the strong conservativity condition $\llbracket\widehat{\sigma}\rrbracket=0$, then for any $\widehat{u}\in\widehat{V}$ and $\mathcal{K}\subset\mathcal{T}_{h}$,

[TABLE]

Consequently, if (11) is multisymplectic and strongly conservative, then it is strongly multisymplectic.

Proof.

Let $e=\partial K^{+}\cap\partial K^{-}$ be an internal facet, where $K^{\pm}\subset\mathcal{K}$. Since $\llbracket\widehat{\sigma}\rrbracket=0$, we have

[TABLE]

Finally, taking the exterior derivative implies that

[TABLE]

so the contributions from internal facets vanish, as claimed. ∎

3.4. The flux formulation for exact solutions

We now apply the theory of the preceding sections to exact solutions of (1), in the sense of distributions, in which case the flux formulation (11) consists of infinite-dimensional function spaces.

Definition \thedefinition.

The exact flux formulation on $\mathcal{T}_{h}$ is the flux formulation (11) associated to the function spaces

[TABLE]

along with the flux functions $\Phi_{K}(u,\sigma,\widehat{u},\widehat{\sigma})=\widehat{\sigma}-\sigma$, for $K\in\mathcal{T}_{h}$.

The next theorem uses a domain-decomposition-type argument to relate the exact flux formulation to solutions of (1), in the sense of distributions, defined over certain global function spaces on $U$.

Theorem 3.4**.**

The element $(u,\sigma,\widehat{u},\widehat{\sigma})\in V\times\Sigma\times\widehat{V}\times\widehat{\Sigma}$ is a solution to the exact flux formulation if and only if

[TABLE]

is a solution to (1) on $U$, in the sense of distributions, and $(\widehat{u},\widehat{\sigma})$ are the exact traces $\widehat{u}=u\rvert_{\mathcal{E}_{h}}$ and $\widehat{\sigma}\rvert_{\partial K}=\sigma\rvert_{\partial K}$ for all $K\in\mathcal{T}_{h}$.

Proof.

Suppose $(u,\sigma,\widehat{u},\widehat{\sigma})\in V\times\Sigma\times\widehat{V}\times\widehat{\Sigma}$ is a solution to the exact flux formulation.

In particular, (11a) holds for all \tau\in\bigl{[}C_{c}^{\infty}(K)\bigr{]}^{mn}, which immediately gives $\partial_{\mu}u^{i}=\phi^{i}_{\mu}$ on $K$, in the sense of distributions. Taking more general test functions \tau\in\bigl{[}C^{\infty}(K)\bigr{]}^{mn}, not necessarily vanishing on $\partial K$, implies that $\widehat{u}\rvert_{\partial K}=u\rvert_{\partial K}$ in the trace sense. Since this holds for all $K\in\mathcal{T}_{h}$, it follows that $\widehat{u}=u\rvert_{\mathcal{E}_{h}}$. Hence, the trace of $u$ is single-valued on $\mathcal{E}_{h}$, and we may therefore conclude (cf. Brezzi and Fortin [7, Proposition III.1.1]) that u\in\bigl{[}H^{1}(U)\bigr{]}^{n}.

Similarly, taking smooth test functions $v$ in (11b) implies that $-\partial_{\mu}\sigma_{i}^{\mu}=f_{i}$ holds in the sense of distributions, and $\widehat{\sigma}_{i}^{\mu}\mathbf{n}_{\mu}\rvert_{\partial K}=\sigma_{i}^{\mu}\mathbf{n}_{\mu}\rvert_{\partial K}$ holds in the trace sense, for all $K\in\mathcal{T}_{h}$. The flux equation (11c) implies further that $\widehat{\sigma}\rvert_{\partial K}=\sigma\rvert_{\partial K}$, and the conservativity condition (11d) gives $\llbracket\widehat{\sigma}\rrbracket=0$ on $\mathcal{E}_{h}^{\circ}$. Hence, the normal trace of $\sigma$ is single-valued on $\mathcal{E}_{h}$, and we may therefore conclude (cf. Brezzi and Fortin [7, Proposition III.1.2]) that \sigma\in\bigl{[}H(\operatorname{div};U)\bigr{]}^{n}.

The converse is a simple verification of (11). If u\in V\cap\bigl{[}H^{1}(U)\bigr{]}^{n}, then it has a (single-valued) trace \widehat{u}\in\bigl{[}L^{2}(\mathcal{E}_{h})\bigr{]}^{n}=\widehat{V}. Likewise, if \sigma\in\Sigma\cap\bigl{[}H(\operatorname{div};U)\bigr{]}^{n}, then it has a trace \widehat{\sigma}\rvert_{\partial K}\in\bigl{[}L^{2}(\partial K)\bigr{]}^{mn}=\widehat{\Sigma}(\partial K) for each $K\in\mathcal{T}_{h}$; this satisfies $\llbracket\widehat{\sigma}\rrbracket=0$ on $\mathcal{E}_{h}^{\circ}$, so (11c)–(11d) hold. Finally, equations (11a)–(11b) hold by the assumption that $(u,\sigma)$ satisfies (1) in the sense of distributions. ∎

Corollary \thecorollary.

The exact flux formulation satisfies

[TABLE]

for all $\mathcal{K}\subset\mathcal{T}_{h}$, whenever $\phi^{i}_{\mu}\,\mathrm{d}\sigma_{i}^{\mu}+f_{i}\,\mathrm{d}u^{i}$ is closed.

Proof.

It follows immediately from Theorem 3.1 (with $\lambda\equiv 0$) and Theorem 3.3 that the exact flux formulation is strongly multisymplectic, so it satisfies (16). Moreover, Theorem 3.4 gives $\widehat{u}=u\rvert_{\mathcal{E}_{h}}$ and $\widehat{\sigma}\rvert_{\partial K}=\sigma\rvert_{\partial K}$ for all $K\in\mathcal{T}_{h}$, so we may “remove the hats” from (16). ∎

4. Multisymplecticity of particular HDG methods

We now apply the results of Section 3 to the particular HDG methods discussed in Cockburn et al. [13]. Although the flux formulation (11) is more general than that for the class of second-order linear elliptic PDEs they consider, the spaces and fluxes used to define the methods are essentially unchanged.

Throughout this section, we assume that $U\subset\mathbb{R}^{m}$ is polyhedral and that $\mathcal{T}_{h}$ is a triangulation of $U$ by $m$-simplices. We denote by $\mathcal{P}_{r}(K)$ the space of degree-$r$ polynomials on $K\in\mathcal{T}_{h}$ and by $\mathcal{P}_{r}(e)$ the space of degree-$r$ polynomials on $e\in\mathcal{E}_{h}$. We also define spaces of discontinuous polynomial boundary traces,

[TABLE]

and continuous polynomial boundary traces,

[TABLE]

which in Cockburn et al. [13] are called $\mathcal{M}_{h,r}$ and $\mathcal{M}_{h,r}^{c}$, respectively.

4.1. The RT-H method

The hybridized Raviart–Thomas (RT-H) method uses the local function spaces

[TABLE]

i.e., degree-$r$ Lagrange finite elements and Raviart–Thomas finite elements, respectively, for each $i=1,\ldots,n$. The trace spaces are taken to be

[TABLE]

and the local flux functions are

[TABLE]

Note that, although $\widehat{\Sigma}(\partial K)$ is infinite-dimensional, the flux condition (11c) simply states that $\widehat{\sigma}\rvert_{\partial K}=\sigma\rvert_{\partial K}$, so we may eliminate this equation and substitute $\sigma$ wherever $\widehat{\sigma}$ appears in the remaining equations.

Theorem 4.1**.**

The RT-H method is strongly multisymplectic.

Proof.

Observe that $\Phi_{K}$ is a special case of the flux in Theorem 3.1, with $\lambda\equiv 0$. Since \widehat{\Sigma}(\partial K)=\bigl{[}L^{2}(\partial K)\bigr{]}^{mn}, it follows that for any $v\in V(K)$ and $\widehat{v}\in\widehat{V}$, we have $(\widehat{v}-v)\mathbf{n}\rvert_{\partial K}\in\widehat{\Sigma}(\partial K)$. Hence, the hypotheses of Theorem 3.1 are satisfied, so the method is multisymplectic.

To show strong multisymplecticity, let $e\in\mathcal{E}_{h}^{\circ}$ be an arbitrary internal facet, and write $e=\partial K^{+}\cap\partial K^{-}$. Since $\sigma\rvert_{K^{\pm}}\in\Sigma(K^{\pm})$, a standard result on Raviart–Thomas elements (cf. Brezzi and Fortin [7, Proposition III.3.2]) implies that $\sigma_{i}^{\mu}\mathbf{n}_{\mu}\rvert_{e^{\pm}}\in\mathcal{P}_{r}(e)$ for $i=1,\ldots,n$.111This result holds since $e$ lies in an affine hyperplane in $\mathbb{R}^{m}$, so $x\cdot\mathbf{n}$ is constant on $e$. It follows that the degree-$(r+1)$ elements of the Raviart–Thomas space nevertheless have degree-$r$ normal traces. The flux condition (11c) implies $\widehat{\sigma}\rvert_{e^{\pm}}=\sigma\rvert_{e^{\pm}}$, so it follows from the above that \llbracket\widehat{\sigma}\rrbracket\rvert_{e}\in\bigl{[}\mathcal{P}_{r}(e)\bigr{]}^{n}. Since this holds for all $e\in\mathcal{E}_{h}^{\circ}$, the extension by zero of $\llbracket\widehat{\sigma}\rrbracket$ to $\mathcal{E}_{h}$ is in $\widehat{V}_{0}$. Therefore, the RT-H method is strongly conservative, so Theorem 3.3 implies that it is strongly multisymplectic. ∎

Remark \theremark.

To prove multisymplecticity, instead of using Theorem 3.1, we could simply have used $(\widehat{\sigma}-\sigma)\rvert_{\partial K}=0$ to see that (13) holds. However, the argument we have developed here is more general, and we will see in Section 4.3 that it also applies to methods where $\lambda\not\equiv 0$.

4.2. The BDM-H method

The hybridized Brezzi–Douglas–Marini (BDM-H) method uses the local function spaces

[TABLE]

i.e., degree-$(r-1)$ Lagrange finite elements and degree-$r$ Brezzi–Douglas–Marini finite elements, respectively, for each $i=1,\ldots,n$. As in the RT-H method, the trace spaces are taken to be

[TABLE]

and the local flux functions are

[TABLE]

As with RT-H, the flux condition (11c) states that $\widehat{\sigma}\rvert_{\partial K}=\sigma\rvert_{\partial K}$, so we may eliminate (11c) and substitute $\sigma$ for $\widehat{\sigma}$ in the remaining equations.

Theorem 4.2**.**

The BDM-H method is strongly multisymplectic.

Proof.

Since the trace spaces and flux function are the same as the RT-H method, multisymplecticity follows exactly as in Theorem 4.1.

Given any internal facet $e=\partial K^{+}\cap\partial K^{-}$, we have \sigma\rvert_{K^{\pm}}\in\bigl{[}\mathcal{P}_{r}(K^{\pm})\bigr{]}^{mn}, so \widehat{\sigma}\rvert_{e^{\pm}}=\sigma\rvert_{e^{\pm}}\in\bigl{[}\mathcal{P}_{r}(e)\bigr{]}^{mn} and therefore \llbracket\widehat{\sigma}\rrbracket\rvert_{e}\in\bigl{[}\mathcal{P}_{r}(e)\bigr{]}^{n}. Since this holds for all $e\in\mathcal{E}_{h}^{\circ}$, the extension by zero of $\llbracket\widehat{\sigma}\rrbracket$ to $\mathcal{E}_{h}$ is in $\widehat{V}_{0}$. Therefore, the BDM-H method is strongly conservative, so Theorem 3.3 implies that it is strongly multisymplectic. ∎

4.3. The LDG-H methods

There are three variants of the hybridized local discontinuous Galerkin method (LDG-H) discussed in Cockburn et al. [13], corresponding to different choices of the local function spaces. In the setting and notation considered here, these three pairs of spaces are:

[TABLE]

Whichever of these we choose, the trace spaces are

[TABLE]

and the local flux functions are

[TABLE]

where the penalty function $\lambda$ is piecewise constant on $\partial K$, i.e., constant on each facet. Note that since $\lambda$ may be different for each $K\in\mathcal{T}_{h}$, on an internal facet $e=\partial K^{+}\cap\partial K^{-}$, the constants $\lambda\rvert_{e^{\pm}}$ need not be equal to one another. The flux condition (11c) states that

[TABLE]

so we may eliminate (11c) and substitute the expression on the right-hand side wherever $\widehat{\sigma}$ appears in the remaining equations.

Theorem 4.3**.**

The LDG-H method is strongly multisymplectic for each of the choices (17a)–(17c) of local function spaces.

Proof.

The flux $\Phi_{K}$ is precisely that of Theorem 3.1. As shown in the proof of Theorem 4.1, the space \widehat{\Sigma}(\partial K)=\bigl{[}L^{2}(\partial K)\bigr{]}^{mn} satisfies the hypotheses of Theorem 3.1, so the LDG-H method is multisymplectic.

Now, for each of (17a)–(17c), we have

[TABLE]

For any internal facet $e=\partial K^{+}\cap\partial K^{-}$, since $\lambda\rvert_{e^{\pm}}$ are constants, it follows that

[TABLE]

so \llbracket\widehat{\sigma}\rrbracket\rvert_{e}\in\bigl{[}\mathcal{P}_{r}(e)\bigr{]}^{n}. (Note that if $\lambda\rvert_{e^{\pm}}$ were arbitrary $L^{\infty}$ penalty functions, rather than constants, this would not necessarily be true; in fact, $\widehat{\sigma}\rvert_{e^{\pm}}$ might not be polynomials at all.) Since this holds for all $e\in\mathcal{E}_{h}^{\circ}$, the extension by zero of $\llbracket\widehat{\sigma}\rrbracket$ to $\mathcal{E}_{h}$ is in $\widehat{V}_{0}$. Therefore, the LDG-H method is strongly conservative, so Theorem 3.3 implies that it is strongly multisymplectic. ∎

4.4. The CG-H method

The hybridized continuous Galerkin (CG-H) method uses the local function spaces

[TABLE]

i.e., degree-$r$ Lagrange finite elements for $u^{i},v^{i}$ and degree-$(r-1)$ Lagrange finite elements for $\sigma_{i}^{\mu},\tau_{i}^{\mu}$, for $\mu=1,\ldots,m$ and $i=1,\ldots,n$. The trace spaces are taken to be

[TABLE]

and the local flux functions are

[TABLE]

Since $u\in V(K)$, we immediately have $u\mathbf{n}\rvert_{\partial K}\in\widehat{\Sigma}(\partial K)$. Moreover, since $\widehat{V}$ consists of continuous polynomials, the degrees of freedom for $\widehat{V}(\partial K)$ are a subset of those for $V(K)$, so $\widehat{u}\mathbf{n}\rvert_{\partial K}\in\widehat{\Sigma}(\partial K)$ as well. Hence, taking the flux condition (11c) with $\widehat{\tau}=(\widehat{u}-u)\mathbf{n}\rvert_{\partial K}$ implies $\widehat{u}\rvert_{\partial K}=u\rvert_{\partial K}$ for all $K\in\mathcal{T}_{h}$.

Remark \theremark.

The CG-H method is so named because it coincides with the classical continuous Galerkin method with Lagrange finite elements when applied to second-order linear elliptic PDEs of the type considered in Section 2.1, Section 2.2, and Section 2.3, as long as $a=a^{\mu\nu}(x)$ is constant on each $K\in\mathcal{T}_{h}$. See Cockburn et al. [14, 13], Cockburn [12], where the relationship between CG-H and the technique of “static condensation” is also discussed.

More generally, this correspondence also holds for the semilinear system

[TABLE]

Indeed, since (11c) implies $\widehat{u}\rvert_{\partial K}=u\rvert_{\partial K}$, substituting $u$ for $\widehat{u}$ in (11a) implies $\sigma=a\operatorname{grad}u$, as long as $a$ is constant on each $K\in\mathcal{T}_{h}$. (Otherwise, $a\operatorname{grad}u$ is generally not in $\Sigma$.) It follows that $u\in C^{0}(U)\cap V$, so for all test functions $v$ in this same space vanishing on $\partial U$, summing (11b) over $K\in\mathcal{T}_{h}$ gives

[TABLE]

The boundary term vanishes by (11d), since the assumption that $v$ is continuous and vanishes on $\partial U$ implies $v\rvert_{\mathcal{E}_{h}}\in\widehat{V}_{0}$. Hence, $u$ is a solution to the continuous Galerkin method for $-\operatorname{div}(a\operatorname{grad}u)=f(\cdot,u)$.

We now prove that the CG-H method is multisymplectic, although—unlike the other HDG methods considered here—it is not strongly multisymplectic except in dimension $m=1$, when multisymplecticity is just ordinary symplecticity.

Theorem 4.4**.**

The CG-H method is multisymplectic. It is strongly multisymplectic (i.e., symplectic) when $m=1$.

Proof.

Since $(\widehat{u}-u)\rvert_{\partial K}=0$, we see directly that (13) holds, so multisymplecticity follows by Section 3.2.

When $m=1$, facets are simply vertices, so $\mathcal{E}_{h}$ is discrete and finite, and the continuity conditions on $\widehat{V}$ are trivial. Hence, \bigl{[}\mathcal{P}_{r}^{c}(\mathcal{E}_{h})\bigr{]}^{n}=\bigl{[}\mathcal{P}_{r}(\mathcal{E}_{h})\bigr{]}^{n}=\mathbb{R}^{\lvert\mathcal{E}_{h}\rvert n}. The result follows immediately from the trivial observation that $\llbracket\widehat{\sigma}\rrbracket\rvert_{e}\in\mathbb{R}^{n}$ at each internal vertex $e$. ∎

Remark \theremark.

Although $\Phi_{K}$ is the same flux function as in Theorem 3.2, the hypotheses of that theorem do not hold for the CG-H method. Here, $\widehat{\Sigma}(\partial K)$ consists only of $\widehat{\tau}$ whose normal traces are continuous on $\partial K$, while this is not necessarily true of $\tau\rvert_{\partial K}$ for arbitrary $\tau\in\Sigma(K)$.

Example \theexample (CG-H for Laplace’s equation in $\mathbb{R}^{2}$).

Let $m=2$, $n=1$, and consider the mixed form of Laplace’s equation,

[TABLE]

Let us apply the lowest-order CG-H method, with $r=1$, so that for each $K\in\mathcal{T}_{h}$, we have $V(K)=\mathcal{P}_{1}(K)$ and \Sigma(K)=\bigl{[}\mathcal{P}_{0}(K)\bigr{]}^{2}.

As discussed in Section 4.4, we have $\widehat{u}\rvert_{\partial K}=u\rvert_{\partial K}$ and $\sigma\rvert_{K}=\operatorname{grad}u\rvert_{K}$. Therefore, $\widehat{\sigma}\rvert_{\partial K}$ is determined by (11b), which in this case states

[TABLE]

Since $\widehat{\sigma}\rvert_{\partial K}=w\mathbf{n}\rvert_{\partial K}$ for some $w\in V(K)$, we may rewrite this as

[TABLE]

However, $w\rvert_{\partial K}$ is generally not equal to $\operatorname{grad}u\cdot\mathbf{n}\rvert_{\partial K}$, since the latter is piecewise constant (and generally discontinuous) on $\partial K$, whereas the former must be continuous and linear. Instead, we must set up a linear system and solve for $w$ in terms of $u$.

For simplicity, let us suppose that $K$ is isometric to the standard, equilateral reference triangle in $\mathbb{R}^{3}$, defined by

[TABLE]

On $T$, the Lagrange basis of linear “hat functions” simply consists of the coordinate functions $x$, $y$, and $z$, and we can write

[TABLE]

Solving the linear system corresponding to (18) yields, after a calculation,

[TABLE]

The multisymplectic form restricted to $\partial T$ is therefore $(\mathrm{d}u\wedge\mathrm{d}w)\mathbf{n}\rvert_{\partial T}$, so the multisymplectic conservation law states that $\int_{\partial T}\mathrm{d}u\wedge\mathrm{d}w=0$.

Let $e_{ij}$ denote the edge in $\partial T$ going from the $i$th standard basis vector to the $j$th standard basis vector in $\mathbb{R}^{3}$. Using the above expressions for $w$ in terms of $u$, another calculation shows that

[TABLE]

from which it is immediately apparent that

[TABLE]

so the method is indeed multisymplectic on $T$, following Theorem 4.4.

We now show that strong multisymplecticity does not hold for the CG-H method when $m>1$, so the result of Theorem 4.4 is the best we can hope for. The proof uses a counterexample based on Section 4.4.

Proposition \theproposition.

The CG-H method is not strongly multisymplectic.

Proof.

Consider the mixed form of Laplace’s equation on the domain $U\subset\mathbb{R}^{2}$ triangulated by two equilateral triangles, as shown in Figure 1. For simplicity, as in Section 4.4, we suppose that each of these is isometric to the standard reference triangle $T$, so that the edge lengths are all $\sqrt{2}$. Note that $\widehat{V}_{0}=\{0\}$, since the degrees of freedom for $\widehat{V}$ all lie on $\partial U$, so the conservativity condition (11d) is trivial. Hence, there are no constraints on $(u,\sigma,\widehat{u},\widehat{\sigma})$ other than the local conditions (11a)–(11c) on each triangle, as discussed in Section 4.4.

From the calculation in (19), we have

[TABLE]

and similarly,

[TABLE]

so adding these gives

[TABLE]

We conclude that this is nonzero since $u_{1},\ldots,u_{4}$ are independent degrees of freedom that may be varied independently. In particular, if we take variations $(v,\operatorname{grad}v,v,\widehat{\tau})$ and $(v^{\prime},\operatorname{grad}v^{\prime},v^{\prime},\widehat{\tau}^{\prime})$ with

[TABLE]

we have

[TABLE]

Hence, the CG-H method is not strongly multisymplectic. ∎

Remark \theremark.

The failure of strong multisymplecticity, in this example, may also be seen via the failure of the bilinear form $(v,v^{\prime})\mapsto\int_{\partial U}v\widehat{\tau}^{\prime}\cdot\mathbf{n}$ to be symmetric. Indeed, if we write $\mathbf{v}=(v_{1},v_{2},v_{3},v_{4})$ and $\mathbf{v}^{\prime}=(v_{1}^{\prime},v_{2}^{\prime},v_{3}^{\prime},v_{4}^{\prime})$, then we may represent this as the quadratic form,

[TABLE]

which is immediately seen not to be symmetric. Neglecting the scalar factor of $\sqrt{3}/6$, the antisymmetrization of this matrix is

[TABLE]

This is precisely the matrix corresponding to

[TABLE]

which agrees with (20).

Remark \theremark.

Cockburn et al. [13] observe that CG-H can be seen as a limiting case of the LDG-H method (17c) where the penalty $\lambda\equiv+\infty$. They also consider more general LDG-H methods where $\lambda$ is infinite on some elements and finite on others [13, Section 3.4]. However, allowing the penalty to be infinite imposes continuity conditions on $\widehat{V}$, meaning that conservativity (and thus multisymplecticity) is only weak rather than strong.

4.5. The NC-H method

The hybridized nonconforming (NC-H) method uses the local function spaces

[TABLE]

just as in the CG-H method. The trace spaces are

[TABLE]

and the local flux functions are

[TABLE]

Unlike CG-H, $(\widehat{u}-u)\mathbf{n}\rvert_{\partial K}$ is generally not in $\widehat{\Sigma}(\partial K)$, so we cannot conclude that $\widehat{u}$ equals $u$ on $\partial K$, except in the weak sense of (11c).

Theorem 4.5**.**

The NC-H method is strongly multisymplectic.

Proof.

The flux $\Phi_{K}$ is that considered in Theorem 3.2, so it suffices to show that the hypothesis of that theorem holds, i.e.: for all $\tau\in\Sigma(K)$, there exists $\widehat{\tau}\in\widehat{\Sigma}(\partial K)$ such that $\widehat{\tau}_{i}^{\mu}\mathbf{n}_{\mu}=\tau_{i}^{\mu}\mathbf{n}_{\mu}\rvert_{\partial K}$ for $i=1,\ldots,n$. Given $\tau\in\Sigma(K)$, this condition is satisfied by $\widehat{\tau}_{i}^{\mu}=\tau_{i}^{\nu}\mathbf{n}_{\nu}\mathbf{n}^{\mu}\rvert_{\partial K}$, i.e., the projection of $\tau\rvert_{\partial K}$ onto the unit normal, so the NC-H method is multisymplectic.

Next, if $e=\partial K^{+}\cap\partial K^{-}$ is an interior facet, then the definition of $\widehat{\Sigma}(\partial K^{\pm})$ implies that \widehat{\sigma}\rvert_{e^{\pm}}\in\bigl{[}\mathcal{P}_{r-1}(e)\bigr{]}^{n}, so \llbracket\widehat{\sigma}\rrbracket\rvert_{e}\in\bigl{[}\mathcal{P}_{r-1}(e)\bigr{]}^{n}. Since this holds for all $e\in\mathcal{E}_{h}^{\circ}$, the extension by zero of $\llbracket\widehat{\sigma}\rrbracket$ to $\mathcal{E}_{h}$ is in $\widehat{V}_{0}$. Therefore, the NC-H method is strongly conservative, so Theorem 3.3 implies that it is strongly multisymplectic. ∎

4.6. The IP-H methods

We finally consider the special case of the hybridized interior penalty method (IP-H), which—unlike the methods considered above—is somewhat idiosyncratic to semilinear elliptic systems. Consider a system of the form

[TABLE]

which generalizes the scalar ($n=1$) semilinear PDEs discussed in Section 2 to $n\geq 1$. Here, $a=a^{\mu\nu}_{ij}(x)$ is a symmetric, positive-definite $mn\times mn$ matrix with inverse a^{ij}_{\mu\nu}(x)\coloneqq\bigl{(}a_{ij}^{\mu\nu}(x)\bigr{)}^{-1}. These are the de Donder–Weyl equations for the Hamiltonian

[TABLE]

so in particular (21) is a canonical multisymplectic system of PDEs.

For such a system, the IP-H method uses the local function spaces

[TABLE]

For both methods, the trace spaces are taken to be

[TABLE]

Based on the observation that the classical solution satisfies $\sigma_{i}^{\mu}=(a\operatorname{grad}u)_{i}^{\mu}=a_{ij}^{\mu\nu}\partial_{\nu}u^{j}$, these methods take the local flux functions to be

[TABLE]

Thus, the IP-H and “IP-H-like” methods are essentially the LDG-H methods (17b) and (17c), respectively, except with $\sigma$ replaced by $a\operatorname{grad}u$ in the local flux functions. As with LDG-H methods, the penalty function $\lambda\rvert_{\partial K}$ is taken to be piecewise constant on $\partial K$ for each $K\in\mathcal{T}_{h}$, and on internal facets $e=\partial K^{+}\cap\partial K^{-}$, the constants $\lambda\rvert_{e^{\pm}}$ need not be equal to one another. The flux condition (11c) states that

[TABLE]

so we may eliminate (11c) and substitute the expression on the right-hand side wherever $\widehat{\sigma}$ appears in the remaining equations. We also assume, following Cockburn et al. [13], that $a$ is constant on each $K\in\mathcal{T}_{h}$.

Theorem 4.6**.**

For the semilinear system (21) where $a$ is constant on each $K\in\mathcal{T}_{h}$, the IP-H and “IP-H-like” methods are strongly multisymplectic.

Proof.

To prove that the methods are multisymplectic, Section 3.2 states that it suffices to show that (13) is satisfied on each $K\in\mathcal{T}_{h}$. Observe that

[TABLE]

We have previously seen that

[TABLE]

by the symmetry of $\delta$ and the antisymmetry of $\wedge$. Therefore, it remains to show that the terms involving $a_{ij}^{\mu\nu}\partial_{\nu}u^{j}-\sigma_{i}^{\mu}$ vanish as well. Integrating (11a) by parts gives

[TABLE]

and since this holds for all $\tau\in\Sigma(K)$, we can write

[TABLE]

Taking the exterior derivative of both sides yields

[TABLE]

where $a^{ij}_{\mu\nu}\,\mathrm{d}\sigma_{j}^{\nu}\wedge\mathrm{d}\sigma_{i}^{\mu}=0$ by the symmetry of $a$ and the antisymmetry of $\wedge$. On the other hand, taking $\tau=a\operatorname{grad}v$ for $v\in V(K)$, we have

[TABLE]

and since this holds for all $v\in V(K)$, we can write

[TABLE]

Note that we need $a$ to be constant on $K$ in order for \tau=a\operatorname{grad}v\in\bigl{[}\mathcal{P}_{r-1}(K)\bigr{]}^{mn}\subset\Sigma(K) to be an admissible test function for both methods (22a)–(22b). Taking the exterior derivative of both sides yields

[TABLE]

where again we have used the symmetry of $a$ and the antisymmetry of $\wedge$. Rearranging this last equality yields

[TABLE]

which is precisely (13). Hence, the IP-H method is multisymplectic.

To show strong multisymplecticity, let $e=\partial K^{+}\cap\partial K^{-}$ be an internal facet. Then

[TABLE]

so \llbracket\widehat{\sigma}\rrbracket\rvert_{e}\in\bigl{[}\mathcal{P}_{r}(e)\bigr{]}^{n}. (Here, we use that $a\rvert_{e^{\pm}}$ and $\lambda\rvert_{e^{\pm}}$ are constant.) Hence, $\llbracket\widehat{\sigma}\rrbracket$ may be extended by zero to an element of $\widehat{V}_{0}$, so both methods are strongly conservative and thus, by Theorem 3.3, strongly multisymplectic. ∎

5. Conclusion

We have generalized the flux formulation and HDG framework of Cockburn et al. [13] to the much larger family of canonical systems of PDEs (1), which includes nonlinear systems. Within this framework, we have established the multisymplecticity of several HDG methods, when such methods are applied to canonical Hamiltonian systems (2). These methods include “hybridized” versions of several widely-used classes of finite element methods, suggesting that—when multisymplectic structure preservation is desired—these general purpose methods may be used instead of the specialized methods constructed in previous work, which are often limited in their order of accuracy and/or use on unstructured meshes.

It is perhaps not surprising that so many finite element methods are multisymplectic. Indeed, like the Ritz–Galerkin method, the multisymplectic conservation law is intimately related to variational principles (cf. Lawruk et al. [24], Kijowski and Tulczyjew [23], Marsden et al. [27], Gotay et al. [19, 18], Vankerschaver et al. [42]). However, any construction of multisymplectic finite element methods must address two difficulties: first, that the multisymplectic conservation law holds for smooth solutions, while finite element spaces are nonsmooth; and second, that finite element basis functions can have support on several elements, posing an obstacle to expressing a localized, per-element conservation law. The flux formulation of Section 3 provides a solution to both of these difficulties. Hybridization (i.e., introducing separate spaces of boundary traces and fluxes) provides a way to express the multisymplectic conservation law on individual elements, while the exact flux formulation of Section 3.4 provides a bridge between the smooth and non-smooth cases (particularly Theorem 3.4 and Section 3.4).

There are two final points we wish to emphasize about these results, by comparison to the previous work discussed in Section 1.

1. (1)

Multisymplectic HDG methods may be applied to systems of the form (1), whether or not the user is aware of any canonical multisymplectic/Hamiltonian structure. Yet, if such a structure is present, it will automatically be preserved. This is analogous to the case for certain classes of symplectic integrators for ODEs (e.g., symplectic partitioned Runge–Kutta methods), but contrasts with previously-studied multisymplectic methods for PDEs on unstructured meshes (e.g., Lagrangian variational methods). 2. (2)

Many of the previously-constructed multisymplectic methods merely satisfy a “discrete version” of the multisymplectic conservation law, e.g., a finite-difference version of (4) on a lattice. This is true even for the previous work on multisymplectic finite element methods [21, 44, 11], which are only “multisymplectic” in a finite-difference sense on the lattice of degrees of freedom. Moreover, this discrete multisymplectic conservation law may differ from method to method, depending on how the divergence operator is discretized.

By contrast, the multisymplectic conservation law (12) satisfied by these HDG methods is exactly the integral multisymplectic conservation law (5), restricted to the numerical traces and fluxes on $\partial K$ for each element $K\in\mathcal{T}_{h}$, while (16) extends this to arbitrary unions of elements. Furthermore, this multisymplectic conservation law does not differ from method to method: it has precisely the same form for weak solutions, in the exact flux formulation of Section 3.4, as it does for each of the HDG methods of Section 4.

One direction for future work is the application of multisymplectic HDG methods to Hamiltonian time-evolution PDEs. There are two natural ways in which such work might proceed. First, a multisymplectic HDG method might be used to semidiscretize the PDE in space, resulting in a finite-dimensional system of ODEs. Sánchez et al. [41] have recently shown that, when an LDG-H method is used to semidiscretize the acoustic wave equation, the resulting system of ODEs is Hamiltonian, and one may then apply a symplectic integrator in time. Second, one might consider the application of spacetime HDG methods, as in Rhebergen and Cockburn [39, 40], Griesmaier and Monk [20]. In this approach, one would simultaneously discretize space and time by applying the flux formulation of Section 3 to the case where $U$ is a spacetime domain and $\mathcal{T}_{h}$ a decomposition into spacetime elements $K\in\mathcal{T}_{h}$. The results of Section 3 are formulated in sufficient generality to include such methods, so one might apply them to investigate the multisymplecticity of specific spacetime HDG methods.

Acknowledgments

This research was supported in part by the Marsden Fund of the Royal Society of New Zealand and by the Simons Foundation (award #279968 to Ari Stern). We also wish to thank the anonymous referees for their helpful comments and suggestions.