Discrete gradient methods have an energy conservation law

We show for a variety of classes of conservative PDEs that discrete gradient methods designed to have a conserved quantity (here called energy) also have a time-discrete conservation law. The discrete conservation law has the same conserved density as the continuous conservation law, while its flux is found by replacing all derivatives of the conserved density appearing in the continuous flux by discrete gradients.

The discrete gradient method [3,5,7,9] is a powerful and general technique for constructing integral-preserving integrators. It includes projection methods and the average vector field method [6] as special cases. In particular, it is useful for constructing energy-preserving integrators.
For PDEs, a conservation law is more fundamental than a conserved quantity, for the former implies the latter only in the presence of suitable boundary conditions, and the flux is (usually) a local quantity. Discrete gradient methods have been widely used to construct energy-preserving integrators for PDEs and their semidiscretizations; see [2] for a survey. In this paper we show that discrete gradient methods also preserve energy conservation laws (ECLs) when the symplectic structure is constant. Discrete gradient methods are not usually symplectic or variational and the resulting conservation laws developed in this paper are not instances of Noether's theorem. Despite their nonvariational nature, the use of 'conservative' discretizations of systems of conservation laws of mass, momentum and so on-that is, of discretizations with discrete conservation laws-is widespread in numerical methods for PDEs. The typical case is that the conserved density is (one of) the dependent variable(s): for example, is a conservation law with conserved density u and flux f (u). Euler's method is a time discretization in the form of a (time-)discrete conservation law. (The full discretization where ∆ is any difference operator in x, is a fully discrete conservation law, but the spatial discretization is covered by discrete versions of Noether's theorem (see, e.g., [4] and references therein) and is not the subject of this paper.) Another known case is when a conservation law arises in a system of Euler-Lagrange equations in which the Lagrangian admits a symmetry acting only on the dependent variables.
If the symmetry is fairly simple, e.g., if it is linear, then there may exist symplectic or discrete Lagrangian methods that preserve the symmetry and which will have a discrete conservation law. See [8] where this situation is illustrated for the phase symmetry of the nonlinear Schrödinger equation. Discrete gradient methods are a way of ensuring that a time discretization has a conservation law when the conserved density, that we call energy in this paper, is not or cannot be one of the dependent variables, in a way that is independent of any variational structure or symmetry.
A discrete gradient on a vector space V with inner product a ⊤ b is a function satisfying the discrete gradient axiom for all smooth functions H: V → R and all z 0 , z 1 ∈ V . For consistency in the continuous limit, generally one also requires although for conservation properties only the discrete gradient axiom is needed. A common discrete gradient is the Average Value discrete gradient given by In coordinates we write H y for the y-component of ∇H. If z = (x, y), then the discrete gradient axiom becomes We first give an explicit calculation of the ECL for first-order canonical Hamiltonian PDEs in one space dimension and their discrete gradient time discretizations. Proposition 1. Let z(x, t) ∈ R n , let K be an n × n antisymmetric matrix and let H be a Hamiltonian of the form H = H(z, z x ) dx, and consider the equations of motion The discrete gradient method has a discrete ECL. Its conserved energy density, H, is the same as that of the PDE. Its flux is the same bilinear form as that of the continuous ECL, but with the gradients of the energy density replaced by their discrete gradients.
Proof. The continuous ECL can be found as follows.
Applying a discrete gradient method in (z, z x ) to the PDE gives the time discretization The change in energy density over one time step is which establishes the result.
Note that the discrete calculation exactly follows the continuous one, the key step being the discrete gradient axiom. The calculation is entirely local and does not require integration by parts.

Proposition 2.
Let K be any constant (i.e., independent of z) skew-adjoint differential operator and let H = H(z, z x , z xx , . . .) be any differential function of z (and possibly x, although we suppress the x). Let K and H determine the Poisson system is the Euler operator applied to H. The flux associated to the conserved density H is a quadratic form in the derivatives of H. A discrete gradient method applied to this system has a discrete ECL with conserved density H and flux given by the same quadratic form with derivatives of H replaced by discrete gradients.
Proof. From skew-adjointness of K we can write ) Then the flux for the energy conservation law is the quadratic form [10] S(E(H), E(H)) − A (KE(H), H) where The calculation of the discrete ECL now follows as in Proposition 1.
In the literature there has been much attention paid to multisymplectic discretizations of variational PDEs [1]. These often take the multi-Hamiltonian form which has a multisymplectic conservation law (dz ∧ Kdz) t + (dz ∧ Ldz) x = 0. Multisymplectic methods preserve a discretization of this (differential) conservation law, but do not have a energy conservation law. (This is the PDE analogue of the fact that a integrator for a general Hamiltonian ODE cannot be both symplectic and energy-preserving). Our next result covers this case. As in Proposition 2, it could be extended to differential operators K and arbitrary differential polynomials H, but the case considered here shows the basic structure. Proposition 3. Let H = H(z, z x ) and let K be a (possibly singular) antisymmetric matrix. Then PDE Kz t = δH δz has an ECL with conserved density H and flux given by a bilinear form in z t and the derivatives of H. The discrete gradient method has a discrete ECL with conserved density H and flux given by the same bilinear form with z t replaced by (z 1 − z 0 )/∆t and the derivatives of H replaced by discrete gradients.
Proof. First we note that multi-Hamiltonian PDEs take the given form with H = S(z) − 1 2 z ⊤ Lz x and δH δz = ∇S(z) − Lz x . The continuous ECL is found as follows: The key step is the use of the product rule for the derivative of the flux H ⊤ zx z t . Note that there is no requirement for K to be invertible. If K is invertible then one can solve and substitute for z t to give the same flux as in Proposition 1. For the discrete gradient method we have establishing the result.
We now consider the fully discrete case. The existence and derivation of conservation laws for variational and Hamiltonian discretizations is the subject of discrete versions of Noether's theorem, see e.g. [4]. We restrict ourselves here to the observation that if the semidiscretization has a semidiscrete ECL, then a discrete gradient method applied to this semidiscretization will have a fully discrete ECL. In the following proposition, the extension to H = H(z, z x ) is routine, we omit the z-dependence for clarity. Proposition 4. Let z(x, t) ∈ R n , let K be an n × n antisymmetric matrix and let H be a Hamiltonian of the form H = H(z x ) dx with equations of motion has a semidiscrete ECL and the discrete gradient methods applied to this semidiscretization has a fully discrete ECL.
Proof. The semidiscretization is