In numerical ordinary differential equations (ODEs), it is known that all B-series methods (including Runge–Kutta methods) preserve linear invariants, while only certain ones preserve quadratic invariants. Linear invariants arising in physical systems include mass, charge, and linear momentum; quadratic invariants include angular momentum and other momentum maps, as well as the canonical symplectic form for Hamiltonian systems. SeeHairer, Lubich, and Wanner  and references therein.
However, for partial differential equations (PDEs) describing time evolution, it is desirable for a numerical integrator to preserve not only global invariants but also local conservation laws. For instance, the evolution may preserve total mass (a global invariant), but the mass in a particular region may change by flowing through the boundary of the region (a local conservation law). Another example is the canonical multisymplectic conservation law for Hamiltonian PDEs, which is a quadratic local conservation law for the variational equation. Focusing only on global invariants overlooks this more granular, local form of conservativity.
This paper develops a new framework for the preservation of such properties by numerical integrators. We do so by answering a much more general question: When does a numerical integrator preserve the evolution of certain classes of observables (e.g., linear, quadratic), even when those observables are not invariants? This includes not only global invariants, as previously studied, but also local conservation laws and other balance laws encountered in both conservative and dissipative dynamical systems.
The main idea of our approach is summarized as follows. Suppose evolves in a (finite- or infinite-dimensional) Banach space according to . Given a functional
, the chain rule implies thatevolves according to . We say that a numerical integrator is -functionally equivariant if applying it to the augmented system
preserves the relation , i.e., the integrator maps . In other words, the following diagram commutes:
This is weaker than equivariance in the usual sense, since the diagram need only commute for (1), not arbitrary
-related vector fields. Preserving invariants becomes the special case where the augmented equation readsand the integrator leaves constant.
We develop a theory of functional equivariance and show that it provides a useful tool kit for understanding the behavior of (especially affine and quadratic) observables, including local conservation laws and multisymplecticity. The paper is organized as follows:
Section 2 characterizes the functional equivariance of B-series methods and explores some consequences for both conservative and non-conservative dynamical systems. The main result, Theorem 2.8, shows that a method is functionally equivariant for a class of observables if and only if it preserves invariants in that class. In particular, all B-series methods are affine functionally equivariant, and those preserving quadratic invariants are quadratic functionally equivariant.
Section 3 applies this framework to local conservation laws for PDEs and spatially semidiscretized PDEs. In particular, affine/quadratic functionally equivariant numerical integrators are seen to preserve discrete-time local conservation laws for affine/quadratic observables.
Section 4 applies this framework to the multisymplectic conservation law for canonical Hamiltonian PDEs and spatially semidiscretized PDEs. We show that multisymplectic semidiscretization in space, followed by a symplectic integrator in time, yields a multisymplectic method in spacetime. We also show that hybrid finite elements may be used for multisymplectic semidiscretization, generalizing the results of McLachlan and Stern  to time-evolution problems.
Finally, Section 5 extends these results from B-series methods to additive and partitioned methods, including additive/partitioned Runge–Kutta methods.
2. Functional equivariance of B-series methods
We begin by considering the case of affine equivariant numerical integrators, which McLachlan, Modin, Munthe-Kaas, and Verdier  showed are precisely B-series methods. In Section 5, we will discuss how these results generalize to additive and partitioned methods, including splitting/composition methods and additive/partitioned Runge–Kutta methods.
2.1. Basic definitions and results
Let be a numerical integrator, whose application to a vector field with time-step size gives a map , . All the methods we will consider have , so it suffices to consider the integrator map with unit time step.
Given an affine map , a pair of vector fields and is -related if . A numerical integrator is affine equivariant if for all -related and , all affine maps , and all Banach spaces and .
Given a Gâteaux differentiable map and , define by . We say that a numerical integrator is -functionally equivariant if for all . That is, if , then . Given a class of maps , we say that is -functionally equivariant if this holds for all and all Banach spaces and .
This is a slight generalization of the situation considered in the introduction: may now be any Banach space rather than , and is only required to be Gâteaux differentiable rather than . Note that is precisely the vector field corresponding to the augmented system (1).
Example 2.4 (Runge–Kutta methods).
An -stage Runge–Kutta method has the form
where and are given coefficients defining the method. When this method is applied to the augmented system (1), we augment the method by
Note that the internal stages are not needed, since the augmented vector field depends only on . Hence, for a Runge–Kutta method, -functional equivariance says that
In particular, if is an invariant of , so that for all , then the terms of this sum vanish, and we get .
Every affine equivariant method is affine functionally equivariant.
If is an affine map, then so is . Since the vector fields and in Definition 2.3 are -related, the conclusion follows by affine equivariance. ∎
The converse is generally not true. For instance, the map , which is defined for finite-dimensional , is seen to be affine functionally equivariant but is not affine equivariant except in the weaker sense mentioned in Remark 2.2. This is an example of an “aromatic” series that is not a B-series, cf. Munthe-Kaas and Verdier .
Since we are concerned for the time being with affine equivariant numerical integrators, it is natural to make the following assumptions on . This includes cases where contains not only affine maps but also quadratic or higher-degree polynomial maps.
The class of maps satisfies the following:
contains the identity map for all ;
is a vector space for all and ;
is invariant under composition with affine maps, in the following sense: If and are affine and , then .
As noted in the introduction, preservation of invariants may be seen as a special case of functional equivariance, so one might expect the latter property to be strictly stronger than the former. Perhaps surprisingly, our first main result shows that they are equivalent.
Let be affine equivariant and satisfy creftype 2.7. Then preserves -invariants if and only if is -functionally equivariant.
Suppose preserves -invariants. Given , it follows from creftype 2.7 that is in . This is an invariant of the augmented vector field , since . Hence, preservation of -invariants implies satisfies , i.e., . In particular, implies .
Conversely, suppose is -functionally equivariant. If is an invariant of , then the augmented vector field is , and -functional equivariance implies . However, any constant functional is also in and has the same augmented vector field , so for all . Hence, . ∎
For B-series methods, the following statements hold:
Every B-series method is affine functionally equivariant.
B-series methods preserving quadratic invariants (e.g., Gauss–Legendre collocation methods) are quadratic functionally equivariant.
No B-series method is cubic functionally equivariant.
2.2. Strong equivariance vs. functional equivariance
There is a stronger notion of -equivariance, based on a straightforward generalization of Definition 2.1 to nonlinear maps . Two vector fields and are -related if for all , and is -equivariant if whenever this is the case.
To illustrate the distinction with functional equivariance, we now show that the implicit midpoint method is not quadratic equivariant in this stronger sense, although Corollary 2.92 tells us that it is quadratic functionally equivariant. Let , and observe that the vector fields
are -related. Applying the implicit midpoint method with time step size gives
Since for , the method is not -equivariant. On the other hand, applying the method to the augmented equation with gives
which illustrates that the method is -functionally equivariant.
Essentially, functional equivariance requires only that commute with particular pairs of related vector fields, while strong equivariance requires that it commute with all such pairs.
2.3. Affine equivariance and closure under differentiation
In addition to invariants and observables that depend on itself, we are often interested in those that depend on variations of . We say that is a variation of if satisfy
whose flow is the derivative of the flow of . An especially important example is the canonical symplectic form for Hamiltonian ODEs, which can be understood as a quadratic invariant depending on two variations of .
A numerical method is said to be closed under differentiation if the method applied to (2) is the derivative of , i.e., .
Bochev and Scovel  showed that Runge–Kutta methods are closed under differentiation, from which it follows that those preserving quadratic invariants are symplectic integrators. The same argument can be applied to B-series methods, where closure under differentiation can be established by showing that it holds for all trees [15, Theorem VI.7.1]. Here, we present a new, tree-free proof that uses only affine equivariance, and which will readily generalize to additive and partitioned methods in Section 5.
Affine equivariant numerical integrators are closed under differentiation.
Given , consider the system
corresponding to . Since is -related to , where is either of the projections or , it follows that . Now, let and take , giving the augmented system
By Proposition 2.5, applying to this system yields
Finally, let and take the limit as . ∎
Let be an affine equivariant numerical integrator preserving -invariants. Given , define by
Then , where and .
It is trivial to extend Corollary 2.12 to the case where depends on two or more variations of , e.g., where and .
Before discussing applications to conservation laws for PDEs, which will be the subject of Section 3, we first illustrate some examples of functional equivariance for numerical ODEs.
2.4.1. Hamiltonian systems
Suppose is equipped with a Poisson bracket . Given , the corresponding Hamiltonian vector field is determined by the condition for . That is, the augmented system (1) can be written
Hence, if is -functionally equivariant, then applying to this system gives a “discrete version” of for . For a Runge–Kutta method, this has the form
This holds for any
Poisson bracket, not just the canonical bracket or those for which the Poisson tensor is constant. Preservation of-invariants is the special case where .
2.4.2. Canonical Hamiltonian systems with damping/forcing
Let with canonical coordinates . Consider the system
where is a Hamiltonian and a constant parameter. If has the special form , where is a positive definite mass matrix and is a potential energy function, then energy is dissipated according to , and the parameter dictates the rate of dissipation.
If is also quadratic, then so is , and hence any quadratic functionally equivariant method yields a discrete version of this dissipation law. If is a Runge–Kutta method preserving quadratic invariants, then this has the form
There is no reason to restrict to linear damping: if we replace the damping term in (3) by an arbitrary forcing term , then we obtain
where the sum on the right-hand side approximates the work done by . When is not quadratic, the identities above generally do not hold. However, since the kinetic energy functional is quadratic, we still get the weaker identity
where the sum approximates work done by both conservative and non-conservative forces.
2.4.3. Monotone observables
Suppose is such that , so is monotone decreasing. If is an -functionally equivariant Runge–Kutta method with , then
so is also monotone decreasing along the numerical solution. Conversely, any method with this monotonicity property also preserves -invariants, and is thus -functionally equivariant, since is an invariant when are both monotone decreasing.
For Runge–Kutta methods, the additional condition is needed to get monotonicity. Functional equivariance alone is not sufficient. We are not aware of a more general version of this condition for arbitrary B-series methods.
An immediate consequence is the known B-stability of Runge–Kutta methods preserving quadratic invariants with . If is a Hilbert space with inner product , consider
2.4.4. Symplectic and conformal symplectic systems
Suppose that is a continuous bilinear form on . Let and each be variations of , so that satisfy
Then is a quadratic functional of this augmented system, evolving according to
where is the Lie derivative along of at [1, Theorem 6.4.1].
If preserves quadratic invariants, then we may apply quadratic functional equivariance to describe the numerical evolution of . Taking to be
where , , and . Furthermore, this implies
where is the pullback of by at . For a Runge–Kutta method preserving quadratic invariants, (5) takes the form
and the sum on the right-hand side expresses the difference between and at .
In particular, suppose that is antisymmetric and nondegenerate, so that is a symplectic vector space. If is a symplectic vector field, satisfying , then we recover the result of Bochev and Scovel  that if preserves quadratic invariants, then , i.e., is a symplectic integrator. An interesting generalization is the case of conformal symplectic vector fields, satisfying for some constant , of which (3) is a canonical example; see McLachlan and Perlmutter . In this case, (6) becomes
which can be seen as an approximate conformal symplecticity relation. However, generally does not equal exactly unless ; see McLachlan and Quispel [19, Example 7] for a counterexample when is the implicit midpoint method.
The arguments above apply without modification if is a vector-valued bilinear form, i.e., a continuous bilinear map for some Banach space .
3. Application to conservation laws in PDEs
In this section, we apply the general results of Section 2 to local conservation laws in time-evolution PDEs. We also consider discrete conservation laws in numerical PDEs, when semidiscretization in space is combined with numerical integration in time.
3.1. General approach and examples
Let correspond to a time-dependent system of PDEs on a domain , where the Banach space is a function space (or product of function spaces) on . Suppose that solutions satisfy a local conservation law, in the form
where and depend on . The notation is deliberately suggestive of Maxwell’s equations, where is charge density, is current density, and (7) is local conservation of charge.
From Theorem 2.8 and Corollary 2.9, we immediately obtain a powerful general statement about preservation of local conservation laws under numerical integration. If , where and is an appropriate space of densities, then , and thus an -functionally equivariant integrator satisfies a discrete-time version of (7). For instance, a Runge–Kutta method preserving -invariants satisfies
We note that, while is required to be related to by a functional in (e.g., is affine or quadratic in ), no such restriction is placed on . In particular, all B-series methods preserve affine local conservation laws, while those preserving quadratic invariants also preserve quadratic local conservation laws. In the case of symplectic Runge–Kutta methods, Frasca-Caccia and Hydon [12, Section 3.1] recently proved this by a direct computation, whereas here it is seen as a particular instance of quadratic functional equivariance.
In addition to the differential form of the conservation law (7), we may also integrate over a compact subdomain and apply the divergence theorem to get
where denotes the outer unit normal to . This may be seen as an integral form of the conservation law (7). In this case, if with , then an -functionally equivariant method satisfies a discrete-time version of (8). In the case of a Runge–Kutta method preserving -invariants, this has the form
Maxwell’s equations in consist of the vector evolution equations
along with the scalar constraint equations
and the constitutive relations and . Here, and are the electric and magnetic fields, and are the electric and magnetic flux densities, and are the electric permittivity and magnetic permeability tensors, and and are charge and current density.
Taking the divergence of the first evolution equation, we see that is a local invariant, so the constraint is preserved by the evolution. Next, interpreting the constraint to define as a function of , we see that taking the divergence of the second evolution equation gives the local conservation law . Since and are both linear in , any B-series method will preserve the constraint , together with a discrete version of the conservation law relating and .
Consider the nonlinear Schrödinger (NLS) equation,
A direct computation shows that solutions satisfy
which is a local conservation law for the quadratic functional . Since the implicit midpoint method is quadratic functionally equivariant, it follows that this conservation law is preserved, in the sense that
More generally, for any quadratic functionally equivariant Runge–Kutta method,
Consider the wave equation , written as the first-order system
If is a solution, then
which is a local conservation law for the quadratic functional . Hence, applying any quadratic functionally equivariant Runge–Kutta method gives
3.2. Discrete conservation laws in numerical PDEs
In practice, of course, we will not be applying numerical integrators to infinite-dimensional function spaces. Rather, we typically first semidiscretize in space (e.g., using a finite difference, finite volume, or finite element scheme), yielding a system on a finite-dimensional vector space to which we can apply a numerical integrator.
Suppose the spatial discretization scheme is such that it preserves a semidiscrete conservation law , where is some “discrete divergence” operator. Then it follows from the previous section that, if for some , then applying an -functionally equivariant numerical integrator yields a fully discrete conservation law corresponding to (7). We illustrate this with a few examples, which are semidiscretized versions of those considered in the previous section.
Nédélec  introduced a finite element semidiscretization of Maxwell’s equations, in which and are approximated by piecewise-polynomial vector fields and . This method may be written as
where , , and is any vector field from the same space as .
Taking the divergence of the first equation gives , so the constraint is preserved by the evolution. For the second, when for a piecewise-polynomial scalar field , we get
which we may write as . Thus, taking implies the semidiscrete charge conservation law . (See Berchenko-Kogan and Stern  for a hybridization of Nédélec’s method that preserves a stronger form of this conservation law, using rather than .) Since and are linear in , any B-series method will preserve exactly and give a discrete-time version of the charge conservation law relating and .
For the one-dimensional NLS equation, the finite-difference semidiscretization
satisfies the semidiscrete local conservation law
where the right-hand side is a difference of midpoint approximations to . Hence, a discrete-time version of this conservation law is preserved by any B-series method that preserves quadratic invariants.
For the one-dimensional wave equation, consider the finite-difference semidiscretization
If we define
which is a finite-difference approximation to , then a short calculation gives the semidiscrete conservation law
where the right-hand side is a difference of midpoint approximations to . As in the previous example, a discrete-time version of this conservation law is therefore preserved by any B-series method that preserves quadratic invariants.
3.3. Remarks on quadratic conservation laws arising from point symmetries
Conservation laws with quadratic densities are common in partial differential and differential-difference equations because of their association with linear symmetries of Hamiltonian PDEs. (See, e.g., Olver .) However, not all such symmetries are easily preserved under semidiscretization. We focus here on affine point symmetries, those arising from actions on the field variables.
For example, the one-dimensional NLS equation may be written in the form
where and is the variational derivative with respect to . The integrand is called the Hamiltonian density. Observe that is invariant under the diagonal action , where . This point symmetry leads to the local conservation law for in Example 3.2. More generally, any Hamiltonian density of the form has the same point symmetry, and hence has a local conservation law for .
Similarly, the one-dimensional semidiscretized NLS equation in Example 3.6 can be written
where the summand can be viewed as a discrete Hamiltonian density . The invariance of under the point symmetry , , yields the semidiscrete local conservation law for obtained in Example 3.6. More generally, we get such a local conservation law whenever the discrete Hamiltonian density has the form .
A related example involves orthogonal (rather than unitary) point symmetry. Suppose and its conjugate momentum both take values in , and let act by . Then any Hamiltonian density that depends only on the 10 invariants , , is invariant and thus has a local conservation law for . Like the point symmetry discussed above, this point symmetry is preserved under a wide class of lattice semidiscretizations, which have corresponding semidiscrete quadratic conservation laws.
By contrast with point symmetries, symmetries that involve spatial translations are typically broken by semidiscretization. However, special semidiscretizations can be constructed that preserve versions of the associated conservation laws, although these are generally not symplectic. An example is provided by the Korteweg–de Vries equation
which has a local conservation law with . The semidiscretization
has a semidiscrete conservation law with density only for the parameter (Ascher and McLachlan ).
Frasca-Caccia and Hydon  give general techniques for constructing finite-difference semidiscretizations preserving several local conservation laws—linear, quadratic, or otherwise—with many examples. When such methods are used in conjunction with B-series methods for time integration, it follows from Theorem 2.8 and Corollary 2.9 that affine local conservation laws are always preserved in a discrete sense, while quadratic local conservation laws are preserved by any B-series method that preserves quadratic invariants.
4. Multisymplectic integrators
In this section, we apply the foregoing theory to the multisymplectic conservation law for canonical Hamiltonian PDEs and its preservation by numerical integrators. Since this is a quadratic local conservation law depending on variations of solutions, it follows that B-series methods preserving quadratic invariants also preserve a discrete-time version of the multisymplectic conservation law. Furthermore, we discuss techniques for spatial semidiscretization that preserve a semidiscrete multisymplectic conservation law, reviewing some known results for finite-difference semidiscretization and introducing new results for finite-element semidiscretization. Consequently, when such methods are used in conjunction with B-series methods preserving quadratic invariants, the resulting method will satisfy a fully discrete multisymplectic conservation law.
4.1. Canonical Hamiltonian PDEs
Before discussing the canonical Hamiltonian formalism for time-evolution PDEs, we first quickly recall the stationary (time-independent) case, following the treatment in McLachlan and Stern .
where and . Here and henceforth, we use the Einstein index convention of summing over repeated indices; for instance, has an implied sum over and therefore corresponds to the divergence of the vector field .
Now, for time-dependent problems, we let and depend on , and we introduce an additional unknown field . The de Donder–Weyl equations for , , are then given by
For , the de Donder–Weyl equations are not in the form , since we have expressions for and but not . To deal with this, we assume that the second equation of (10) defines as an implicit function of , , , , and . By the implicit function theorem, this is true (at least locally) if the matrix is nondegenerate. Therefore, we may eliminate the second equation and substitute this expression for into the other two equations. Assuming the Hamiltonian does not depend on , this gives a system of the form with .
Let , so that and are scalar fields and is a vector field on , and take . Then the de Donder–Weyl equations are
Eliminating the second equation and substituting into the third, we obtain the first-order form of the wave equation with , as in Example 3.3.
4.2. The multisymplectic conservation law
For Hamiltonian ODEs, the symplectic conservation law is a statement about variations of solutions to Hamilton’s equations. Similarly, for Hamiltonian PDEs, the multisymplectic conservation law is a statement about variations of solutions to the de Donder–Weyl equations.
Let be a solution to (10). A (first) variation of is a solution to the linearized problem