Noether-Type Discrete Conserved Quantities Arising from a Finite Element Approximation of a Variational Problem

In this work, we prove a weak Noether-type Theorem for a class of variational problems that admit broken extremals. We use this result to prove discrete Noether-type conservation laws for a conforming finite element discretisation of a model elliptic problem. In addition, we study how well the finite element scheme satisfies the continuous conservation laws arising from the application of Noether’s first theorem (1918). We summarise extensive numerical tests, illustrating the conservation of the discrete Noether law using the p-Laplacian as an example and derive a geometric-based adaptive algorithm where an appropriate Noether quantity is the goal functional.


Introduction
Symmetries are an extremely important and continually occurring feature of differential equations arising from many applicable areas, including mathematical physics, meteorology and differential geometry, that was first developed by Sophus Lie for the purpose of studying solutions of differential equations in the late 19th century [Lie71,c.f.].
Noether's (1st) Theorem [Noe71] is a striking result which in the continuous setting connects these symmetries with conservation laws associated to the Euler-Lagrange equations of a variational problem. Roughly, the theorem states that given a variational problem with an underlying symmetry, there exists a natural conservation law associated to it. For example, a symmetry of translation with respect to the spatial coordinates results in a conservation of linear momentum, a symmetry of rotation results in conservation of angular momentum and a symmetry of translation with respect to the temporal coordinate gives a conservation of energy. A famous example from meteorology is that of potential vorticity. This is a conservation law arising from a particle relabelling pseudo-group symmetry. This quantity is extremely important in studying the evolution of a cyclone.
The work of Noether has gained public attention recently with the publication of an article in the New York Times [Ang12] where the result is "consider [ed] . . . as important as Einstein's theory of relativity; it undergirds much of today's vanguard research in physics, including the hunt for the almighty Higgs boson." In the discrete setting, Noether's Theorem has been studied in terms of difference equations [Dor01,HM04], where it was shown that a discrete equivalent of the conservation law holds when a discrete symmetry was built into the discrete Lagrangian. In this work we turn our attention to the finite element method (FEM). FEMs form one of the most successful numerical methods for approximating the solution to partial differential equations (PDEs) [Bra01,BS94,Cia78,c.f.]. A topic which has been the subject of much ongoing research is that of constructing FEMs which inherit some property of the continuous problem. The notion of discretisations inheriting some geometric property from the continuous problem can be seen as a generalisation of geometric integration [HLW06,c.f.] to the case of PDEs and is a rapidly developing area of research. Some of the properties studied in the discretisation of PDEs are the same as in the geometric integration of the ODE, for example the Hamiltonian structure of a given problem. Others are based on completely new notions, for example, the recent development of the discrete exterior calculus [AFW10], which, as the name suggests, is a discrete equivalent to the Cartan based exterior calculus. This has allowed for a rigorous description of discrete differential forms and the associated discrete function spaces as a discrete differential complex. This provides a framework which may be used as a first step in the construction of a variational complex in a similar light to that developed in [HM04] for difference equations. A first step in this direction was taken in [Man06]. A review of some of the huge quantity of topics arising from this area, including Lie group integrators, discrete gradient methods as well as FEMs for differential forms is given in [CMKO11].
As opposed to geometric integrators, the term used for numerical methods with some geometric property of an ODE, the methods for PDEs are generally called mimetic methods. The class of FEMs which fall under the mimetic framework are the mixed methods, for example the Raviart-Thomas scheme [RT77]. It is not only FEMs which fall under the category of mimetic methods, in fact there are finite difference (FD) [BBL09] and finite volume (FV) schemes which are characterised as mimetic. Note that also there is an intrinsic relationship between each of them. For example, an appropriate choice of quadrature for the Raviart-Thomas finite element scheme results in the mimetic FD scheme [CMR09].
The Lagrangian piecewise polynomial FEM is not a mimetic method. Most standard methods cannot inherit geometric properties of the continuous PDE. There is an underlying algebraic condition which must be satisfied for these properties to be inherited by the approximation scheme [HM04,Man06].
The purpose of this paper is to show that variational numerical problems have their own conservation laws which derive from the same principle as that discovered by Noether, giving rise to discrete (numerical) forms of conservation laws which are automatically preserved by the scheme.
The classical Noether theorem is only applicable to classical solutions of the variational problem. As such we derive weaker versions of the theorem applicable to a wider class of solutions to the problem, including the broken extremals. We will discuss how these laws are naturally passed down to the Lagrangian finite element scheme and hence quantify the discrete Noether quantity associated to this FEM. That is, we write the exact Noether quantity for this discretisation, in the same spirit as [HM04].
We will also study how well the Lagrangian finite element scheme satisfies the strong conserved quantities arising from Noether's Theorem measured in an appropriate weak norm. That is, we consider how well this finite element scheme approximates the Noether conservation law for the continuous problem (when one exists). We will also present some interesting numerical results, quantifying the deviation of the approximation in terms of a computable estimator which we are able to use to construct an h-adaptive scheme (local mesh refinement) aimed at minimising the violation of the smooth conservation law to a user specified tolerance.
The paper is set out as follows: In §2 we introduce some fundamental notation and the model problem we consider. In §3 we briefly describe Noether's Theorem and the background material needed. To illustrate its application we apply the Theorem to a simple model problem. In §4 we weaken the invariance criterion on which the classical Noether Theorem is based, ultimately allowing us prove two versions of the theorem applicable to weaker solutions of the problem. In §5 we discuss how the results of §4 can be passed down to give discrete counterparts to our weak Noether's Theorem. We perform numerical experiments to demonstrate that the quantities derived are indeed conserved at the discrete level. We also discuss trivial Lie group actions (those of translation with respect to the dependent variable) and how the mimetic methods relate to this case. Finally, in §6 we study the properties of the finite element solution with respect to the original (strong) Noether Theorem. We also detail an interesting numerical result by constructing a computable estimator, aimed at measuring the violation of the strong Noether Theorem in a weak norm for the Lagrangian finite element scheme. We perform some numerical experiments demonstrating that there is a superconvergence of the estimator over the finite element approximation of the solution to the Euler-Lagrange equations. We then proceed to test an adaptive scheme based on the estimate allowing us to minimise the discrete violation of the continuous conserved quantity up to user specified tolerance.

Notation
Let Ω ⊂ R d be a bounded domain with boundary ∂Ω. We begin by introducing the Sobolev spaces [Cia78,Eva98] L p (Ω) = where α = {α 1 , . . . , α d } is a multi-index, |α| = d i=1 α i and derivatives D α are understood in a weak sense. We pay particular attention to the cases k = 1, 2 and Let L = L(x, u, ∇u) be the Lagrangian. We will let be known as the action functional. The problem arising from the calculus of variations is to seek a function extremising the action functional. For simplicity we will consider the minimisation problem, that is, to to find u ∈ (2.9) Note that we are implicitly coupling the minimisation problem with homogenous Dirichlet boundary conditions. We will use the notation that

Noether's First Theorem
For the reader's benefit we will briefly describe Noether's first theorem in the continuous, smooth case and necessary background material. We assume, in this section, that L is smooth and the minimisation problem (2.9) has a solution (not necessarily unique) which is at least C 2 (Ω), i.e., smooth enough to satisfy the Euler-Lagrange equations (2.12) 3.1. Definition (one-parameter group). The transformation is said to be a one-parameter group if the following conditions hold (1) The parameter choice of = 0 yields the identity, i.e., (2) The inverse is given by the parameter − , i.e., (3) The transformation is closed under composition, i.e., if then ( x, u) = (Ξ(x, u; + δ), Φ(x, u; + δ)) .

3.2.
Definition (infinitesimal). The infinitesimals, ξ(x, u) and φ(x, u) of the one parameter group are defined as 3.3. Definition (characteristics). We define the characteristics, which are given in terms of the infinitesimals of the group, to be x ∈ Υ} be the graph of u over a subdomain such that Υ ⊂ Ω. Also let Υ Ξ = Ξ(Γ; ), then the transformation (3.1) is said to be a variational symmetry if holds for any smooth subdomain Υ of Ω.
3.5. Theorem (infinitesimal invariance [Olv93, Thm 4.12]). A variational symmetry group with infinitesimals ξ, φ and characteristics Q of the action functional For the problem we consider in this work, that of a first order Lagrangian, the conservation law, C , takes the form Proof Using the result of Theorem 3.5 we have that Noting by the product rule that then it holds that (3.13) Since the quantity are the Euler-Lagrange equations associated with the action functional (2.8) they vanish, concluding the proof with as required.
3.7. Remark (the form of C ). It is clear from the identity (3.8) that for our model problem C = C (u, ∇u).

3.8.
Remark (the beauty of the theorem). What makes Theorem 3.6 truly remarkable is its constructive nature. For completeness we will give an example of the construction of C for the Laplacian.
3.9. Example (Laplace's problem). Let us consider the case f = f (|x|) then the Lagrangian, is invariant under the rotational group SO(d). For simplicity we restrict d = 2, we calculate the infinitesimals from the group of rotations, note that in this case Φ ≡ 0 and (3.17) It then holds that In this case the characteristic of the group of rotations is Making use of Theorem 3.6 we see is a conservation law over solutions of L [u] = 0.
3.10. Remark (trivial Lie group actions). For any variational problem, the Euler-Lagrange equations, as already mentioned, are given in variational (or divergence) form. As such, if we assume that L does not depend on u, that is L = L(x, ∇u), then the Euler-Lagrange equations themselves can be a Noether conservation law. Indeed, consider the case of Example 3.9 with f ≡ 0. It is clear by definition that ∆u = div(∇u) = 0 is a conservation law. It arises from Noether's Theorem under the trivial Lie group action, that of translation in the dependent variable For this action, the infinitesimals are ξ = 0 and φ = 1.

Noether's theorem for weak solutions
Noether's Theorem (Theorem 3.6) as it is stated in §3 only makes sense for classical solutions of the Euler-Lagrange equations (2.12). We wish to "weaken" the theorem such that it is applicable to extremals which are W 1 ∞ (Ω) with jump discontinuities, the so called broken extremals [GH96]. We begin by defining a weaker invariance condition than that of Definition 3.4. 4.1. Definition (weak variational symmetry). The transformation (3.1) is said to be a weak variational symmetry if holds over the domain Ω.

4.2.
Remark (strong symmetry ⇒ weak symmetry). We note that any strong variational symmetry is also a weak symmetry but the converse is not true.
4.3. Theorem (Noether type conserved quantities for weak variational symmetries). Suppose that the variational problem (2.9) has a weak variational symmetry and that the minimiser is smooth, i.e., u ∈ C 2 (Ω). Let φ and ξ be the infinitesimal generators of the symmetry as in Definition 3.2. Then 4.4. Remark (structure of (4.2)). The weak conservation law given in Theorem 4.3 has a very clear structure. The first two terms represent the weak Euler-Lagrange equations, which do not vanish in general since the infinitesimal φ does not necessarily belong to the correct space of test functions. The last two terms essentially represent the weak conservation law itself.
Since u is an extremal we have that o (0) = 0 and we may compute the weak variation to be For the inner variation we vary the independent variables, to that end let i : R + → R be another function defined such that Again we have that i (0) = 0, we then look to compute We have that where ‹ ∇ denotes the spatial gradient operator with respect to x. After a change of variable and using the invariance criterion (4.1) we see that Subsituting (4.8) into (4.6) and taking the sum of the inner and outer variations concludes the proof.
4.5. Theorem (strong conservation law ⇒ weak conservation law). Let the variational problem (2.9) have a variational symmetry in the sense of Definition 3.4 and that the minimiser to the variational problem is smooth u ∈ C 2 (Ω), then (4.2) holds.
Proof Recall from (3.13) we have that Using Definition 3.3 of the characteristics, Q we see (4.10) which after some simplification yields Å ∂L Integrating over the domain concludes the proof. Now we have developed the framework sufficiently to state our main result in this section. Here we are concerned with broken extremals, that is, functions whose derivatives have finitely many jump discontinuities. 4.6. Definition (broken extremal). An extremal, u ∈ C 0 (Ω), to the problem (2.9) is said to be a broken extremal if it is piecewise C 2 (Ω) over the domain Ω with bounded derivative. That is, Ω can be decomposed into finitely many open subsets, (1) the subsets make up the entire domain, i.e., Ω = i Ω i , (2) they are non-overlapping, i.e., Ω i ∩ Ω j = ∅ and (3) the solution is smooth over each of the subsets, i.e., u ∈ C 2 (Ω i ) ∩ W 1 ∞ (Ω). 4.7. Definition (skeleton and jumps). We define to be the skeleton of the decomposition. Let n i be the outward pointing normal to Ω i , we then define jumps of scalars and vector valued functions as v := v| Ω1 n 1 + v| Ω2 n 2 (4.13) v := (v| Ω1 ) n 1 +(v| Ω2 ) n 2 , (4.14) respectively.
4.8. Definition (piecewise variational symmetry). Let u be a broken extremal to the variational problem (2.9) and let {Ω i } N i=1 be the decomposed domain of u. Then (3.1) is a piecewise variational symmetry if it is a variational symmetry over Ω i for each i = 1, . . . N . 4.9. Theorem (conserved quantities for broken extremals of the variational problem). Suppose u is a broken extremal to (2.9) and the problem has a piecewise variational symmetry. Then (4.16) The other jump term follows similarly by splitting the inner variation from the proof of Theorem 4.3 into subdomains.

Finite element conservation laws
5.1. Discretisation. In this section we calculate the discrete counterpart to Theorem 4.9 in the finite element context. To that end, let T be a conforming triangulation of Ω, namely, T is a finite family of sets such that (1) K ∈ T implies K is an open simplex (segment for d = 1, triangle for d = 2, tetrahedron for d = 3), (2) for any K, J ∈ T we have that K ∩ J is a full subsimplex (i.e., it is either ∅, a vertex, an edge, a face, or the whole of K and J) of both K and J and (3) K∈T K = Ω. We let E be the skeleton (set of common interfaces) of the triangulation T and say e ∈ E if e is on the interior of Ω and e ∈ ∂Ω if e lies on the boundary ∂Ω.
The shape regularity of T is defined as where ρ K is the radius of the largest ball contained inside K and h K is the diameter of K. We use the convention where h : Ω → R denotes the meshsize function of T , i.e., where h K is the diameter of an element K. We introduce the finite element space where P k denotes the linear space of polynomials in d variables of degree no higher than a positive integer k. We consider k ≥ 1 to be fixed and denote by N := dim V.
The Galerkin approximation to the variational problem (2.9) is to seek U ∈ V ⊂ H 1 0 (Ω) such that The finite element scheme defined by (5.4) is guaranteed to be well posed under some assumptions on L that allow us to invoke the Lax-Milgram Theorem or the more generally applicable inf-sup condition [EG04]. Henceforth we now assume the continuous minimisation problem admits a unique solution.
We may now proceed in deriving a finite element Noether type conservation law. As already seen, the conservation law arises after taking inner and outer variations of the variational problem. The outer variation can be characterised by the following Lemma.
5.2. Lemma (discrete Euler-Lagrange equations). The discrete Euler-Lagrange equations associated to the variational minimisation problem (5.4) are to seek U ∈ V such that where V ∈ V ⊂ H 1 0 is the discrete variation. Since U ∈ V is the discrete minimiser of the energy functional we certainly have that o (0) = 0 and we may explicitly compute this quantity, the first variation.
(5.7) Note that since ∇U is not continuous over the skeleton of the triangulation T we have that ∂L ∂(∇U ) is, in general, not well defined. But ∇U is smooth over the interior of each element. We thus split the integral into elementwise contributions and integrate by parts elementwise. For brevity we note that the Lagrangian L = L(x, U, ∇U ) and drop the dependency.
where div K denotes an elementwise divergence. We now use the identity and hence (5.10) as required.
5.3. Example (discrete Laplace's problem). For example the discrete Euler-Lagrange equations associated to Laplace's problem (Example 3.9) are to find U ∈ V such that where ∆ K is an elementwise Laplacian. Note that if U is a piecewise linear function, the first term of (5.11) is zero. Hence the discrete Laplacian can be completely characterised in terms of the jump of the gradient of U over the internal skeleton. 5.4. Definition (L 2 (Ω) projection operator). We define P V : L 2 (Ω) → V such that for each w ∈ L 2 (Ω) we have (5.12) 5.5. Theorem (conserved quantities over C 0 (Ω)-finite element spaces). Let u be the unique weak extrema to the minimisation problem (2.9) and U be its finite element approximation. Suppose that this problem satisfies a weak variational symmetry. Then the finite element solution satisfies the following (5.14) We have that a finite element minimiser of this energy functional is continuous over the domain Ω but its derivative is not.
Using the discrete Euler-Lagrange equations from Lemma 5.2 testing with φ we need only calculate the inner variation.
The inner variation can be regarded as a change of variables on the independent variable [GH96, §3.3]. In a similar calculation to that of the Proof of Theorem 4.3 we let i( ) = J [U ( x)] with x = x+ ξ. We again split the integral into subdomains to obtain (5.15) as required.
5.6. Applications to the p-Laplacian. In this section we give a numerical verification to Theorem 5.5 for a simple test problem, that of the p-Laplacian where we will restrict p ∈ (1, ∞) for degeneracy issues. The p-Laplacian is the Euler-Lagrange equation of the following minimisation problem: Find u ∈ (5.17) with the (parametrised) action functional J p given by Note that for p = 2 this problem coincides with the standard Laplace's problem (see Example 3.9). For general p it is well known that the problem is uniquely solvable.
The discrete weak formulation associated to the minimisation problem (5.17) is to find U ∈ V such that In this test we choose f such that solves the p-Laplace equation (5.16). We have that f can be written as f = f (|x|) and hence the Lagrangian is invariant under SO(d) group actions. We fix d = 2 and take T to be a structured triangulation of Ω, the unit circle, as given in Figure 1   It is well known [BL94,c.f.] that the finite element approximation (5.19) is well posed and has optimal convergence properties. In Tables 1-3 we show errors, convergence rates and the values of the finite element Noether quantity as written in Theorem 5.5 for various cases of p. The tables also study the experimental order of convergence of the numerical approximation which we now define. 5.7. Definition (experimental order of convergence). Given two sequences a(i) and h(i) 0, i = l, . . . ,, we define experimental order of convergence (EOC) to be the local slope of the log a(i) vs. log h(i) curve, i.e., EOC(a, h; i) := log(a(i + 1)/a(i)) log(h(i + 1)/h(i)) .
(5.22) Table 1. In this test we computationally study the behaviour of the finite element conserved quantity, N [U ], given in Theorem 5.5. To that end we consider the p-Laplacian with p = 3. We fix f such that u is known and is given by (5.20). We compute the piecewise linear (k = 1) finite element approximation given by (5.19). To do this we formulate (5.19) as a system of nonlinear equations, the solution to this is then approximated by a Newton method with tolerance set at 10 −10 . At each Newton step the solution to the linear system of equations is approximated using a stabilised conjugate gradient iterative solver with an algebraic multigrid preconditioner, also set at a tolerance of 10 −10 . Note that N [U ] is below the tolerance of the solvers when the mesh is sufficiently resolved.  Table 2. This is the same problem as in Table 1 with the exception that we take p = 4.  Table 3. This is the same problem as in Table 1 with the exception that we take p = 5. The correct function space setting is to seek u ∈ L 2 (Ω) and p ∈ H div (Ω) := {Ψ : div Ψ ∈ L 2 (Ω)}. A conformal approximation of this problem can be sought using the Raviart-Thomas and piecewise constant finite element pair [RT77], for example. A sufficient condition for the construction of a conformal finite element space of H div (Ω) is that the jumps of the discrete functions vanish over the skeleton of the domain. Taking a generic P ∈ RT 0 we see for our model problem This gives a "weak" Stokes theorem. Recall Remark 3.10 concerned itself with the trivial Lie group action of translation in the dependent variable. For our model problem we have that ∇u is a conservation law. The mimetic scheme weakly enforces this conservation law.

Conservative properties of Lagrangian FEs for strong solutions
In this section we present results concerning the approximability of the strong continuous conservation laws arising from Theorem 3.6. We examine numerically the behaviour of the Lagrangian finite element method. In this sense we wish to measure the quantity div(C [U ]) and evaluate how far it deviates from zero. For clarity of exposition, we will assume henceforth that the continuous minimisation problem takes the form: Find u ∈ H 1 0 such that (6.1) 6.1. Theorem (Bound on the finite element approximation of Noether's laws.). Let u ∈ H 2 (Ω) ∩ H 1 0 (Ω) be a strong extrema to the variational problem (2.9) (and hence a strong solution to the Euler-Lagrange equations (2.12)). Suppose we have that Theorem 3.6 holds under a variational symmetry group with infinitesimals ξ and φ.
In addition assume that we have Then if U ∈ V is the finite element approximation to u there exists a constant C such that

(6.4)
Proof We begin by noting that since u is a strong extremal, Theorem 3.6 holds and we have that div Now, we may use the fact that for a generic ϕ ∈ H 1 Since ϕ was generic we may divide through by ∇ϕ and take the supremum over ϕ then by the definition of the H −1 (Ω) norm By the definition of the Noether quantity C from (3.10) we have that where for clarity we have written the dependencies explicitly. Now for each of the I i we add and subtract appropriate quantities and make use of the triangle inequality. We thus have the following bounds: ∂L(x, u, ∇u) ∂(∇u) (6.9) ã (6.10) and Under assumptions (6.2)-(6.3) and using the fact that φ and ξ are infinitesimals of a smooth Lie group action we have that ã (6.12) ã (6.13) (6.14) Taking the sum of the I i gives the desired result.
6.2. Corollary. Let the conditions of Theorem 6.1 hold under the same variational symmetry group with infinitesimals ξ and φ. In addition assume that the Lagrangian is sufficiently smooth such that both L and ∂L ∂(∇u) are (locally) Lipschitz with respect to the second and third variable then the bound (6.4) can be simplified to (6.15) 6.3. Remark (relating to the p-Laplacian). We may relate Theorem 6.1 and Corollary 6.2 to the p-Laplacian studied in §5.6. We were considering that L was invariant under rotations in the independent variable. In that case we have that φ ≡ 0 and that ξ(x, U ) − ξ(x, u) = 0 (6.16) hence we may infer that The leading term here is ∇U − ∇u which is well known to be = O(h k ).
6.4. Lemma (a computable upper bound for the conservation law). Let u ∈ H 2 (Ω)∩ H 1 0 be a strong extrema to the variational problem (2.9). Suppose that Theorem 3.6 holds and that U is the finite element approximation to u. Then there exists a constant C dependent on the shape regularity of T such that  Proof The proof is a standard aposteriori argument where the quantity of interest is split into regular and singular parts. Recall for the model elliptic problem C [U ] = C (x, U, ∇U ) and hence div C [U ] / ∈ L 2 (Ω). Let · | · denote the H −1 (Ω) -H 1 0 duality pairing then for any ϕ ∈ H 1 0 it holds that Note that for generic p ∈ [L 2 (Ω)] d and v ∈ H 1 (Ω) that Applying a Cauchy-Schwarz inequality followed by a Poincaré inequality together with a trace inequality yields Noting ϕ was a generic function in H 1 0 , dividing through by ∇ϕ and taking the supremum over ϕ yields the desired result. 6.5. Remark. It is worth noting that since the quantity div C [U ] is not orthogonal to V there are no powers of h appearing in the estimate given in Lemma 6.4. What is interesting is that the estimate still converges to zero at the same rate as that of the residual of the problem. 6.6. Numerical experiments. In Figure 2 we show numerically that the estimate given in Lemma 6.4 converges at the same rate as the apriori bound given in Remark 6.3. All numerical experiments are conducted for simplicity on the 2-Laplacian taking Ω = [−1, 1] 2 which is discretised using an unstructured triangulation. 6.7. Adaptive methods. We conclude our numerical experiments in Figure 3 by using the observed convergence of the computable estimator E(U, f ) to construct an adaptive scheme aimed at minimising E(U, f ) (and hence div C [U ] −1 ). The adaptive algorithm we make use of is of standard type (SOLVE → ESTIMATE → MARK → REFINE [SS05, c.f.]) utilising the maximum strategy marking and newest vertex bisection refinement.

Conclusions and outlook
In this work we have given a consise statement of Noether's first theorem, applied to a set of model problems.
We have proved a Noether type theorem for a specific class of weak extremum to a model variational problem, that is, those that are W 1 ∞ (Ω) with finitely many jump discontinuities.   Figure 3. In this experiment we fix f such that u is as in Figure 2. We solve the FE problem on concurrently refined meshes and compute the L 2 (Ω)-error, the H 1 0 -error and the estimate of Noether's conservation law, E(U, f ), from Lemma 6.4. We construct an h-adaptive approximation to the problem with the aim of minimising the violation of the Noether quantity. Note that the estimator decreases far quicker than O(N −k/2 )