Symmetry-preserving finite element schemes: An introductory investigation
Alexander Bihlo Francis Valiquette


fancy




December 30, 2023




Symmetry-preserving finite element schemes:
  An introductory investigation


0.5cm


Alexander Bihlo
    Francis Valiquette

Department of Mathematics and Statistics     Department of Mathematics

Memorial University of Newfoundland    SUNY at New Paltz 

St. John's, NL, Canada  A1C 5S7    New Paltz, NY, USA   12561

abihlo@mun.ca     valiquef@newpaltz.edu 

http://www.math.mun.ca/∼abihlo    http://www2.newpaltz.edu/∼valiquef




0.9
0.25cmKeywords:  Finite elements, geometric numerical integration, invariant discretization, ordinary differential equations, moving frames.
0.25cmMathematics subject classification:   34C14, 65L60


0.5cm


Using the method of equivariant moving frames, we present a procedure for constructing symmetry-preserving finite element methods for second-order ordinary differential equations. Using the method of lines, we then indicate how our constructions can be extended to (1+1)-dimensional evolutionary partial differential equations, using Burgers' equation as an example.
Numerical simulations verify that the symmetry-preserving finite element schemes constructed converge at the expected rate and that these schemes can yield better results than their non-invariant finite element counterparts.




§ INTRODUCTION



Geometric numerical integration is a branch of numerical analysis dedicated to the construction of numerical schemes that preserve intrinsic geometric properties of the differential equations being approximated, <cit.>.  Standard examples include symplectic integrators, <cit.>, Lie–Poisson structure preserving schemes, <cit.>, energy-preserving methods, <cit.>, and general conservative methods, <cit.>. The motivation for considering structure-preserving numerical schemes is that, as a rule of thumb, these integrators provide better global or long term results than their traditional non-geometric counterparts.

In engineering, physics, mathematics, and other mathematical sciences, most  differential equations of interest admit a group of symmetries that encapsulates properties of the equations and their solution spaces.  Over the last 30 years, there has been a considerable amount of work dedicated to the development of finite difference numerical methods that preserve the Lie point symmetries of differential equations, <cit.>. For ordinary differential equations, symmetry-preserving numerical schemes have shown to be very effective, especially when solutions exhibit sharp variations or admit singularities, <cit.>. For partial differential equations, the numerical improvements are not as clear and more work remains to be done, <cit.>. For evolutionary partial differential equations, symmetry-preserving schemes generally require the use of time-evolving meshes which can lead to mesh tangling and other numerical instabilities. To avoid these mesh singularities, various methods have been proposed in recent years, including evolution–projection techniques, invariant r-adaptive methods, and invariant meshless discretizations, <cit.>.

To this day, research on symmetry-preserving numerical schemes has solely focused on finite difference methods. Extending the methodology of symmetry-preserving schemes to other numerical integration techniques such as finite volumes, finite elements, or spectral methods remains to be done.  As such, in this paper we lay out basic ideas for constructing symmetry-preserving finite element methods.  

From a numerical perspective, finite element methods offer several advantages over finite difference methods.  For example, when dealing with complex domains, unstructured grids, or moving boundaries, finite element methods are generally easier to implement than finite difference methods. Also,  the finite element method relies on discretizing a weak form of a system of differential equations and thus has less rigid smoothness requirements than methodologies that rely on the discretization of the strong form of the system. As a result, finite element methods are extensively employed in computational fluid dynamics, structural mechanics, and many other branches of engineering, physics, and applied mathematics.  

As a first attempt to systematically construct symmetry-preserving finite element methods, we restrict our considerations to second-order ordinary differential equations.  As basis functions, we consider piecewise linear (Lagrangian) functions (also called hat functions).  For more accurate schemes, our constructions can easily be applied to higher order Lagrangian interpolants. The ideas developed in this paper can also be applied to higher order ordinary differential equations, with appropriate interpolating functions.  Adapting our results to hierarchical bases, splines and Hermite basis functions, and partial differential equations with multi-dimensional basis functions remains to be considered.

The remainder of this paper is laid out as follows.  In Section <ref> we state the problem we aim to solve in this paper.  Namely, we show via two examples that, in general, the discrete weak formulation of a differential equation will not preserve the symmetries of the original differential equation.  To remedy this situation, we explain how to construct symmetry-preserving finite element schemes using the method of equivariant moving frames.  The basic moving frame constructions, adapted to the problem at hand, are introduced in Section <ref>.  The main results of this paper are found in Section <ref>, where we provide an algorithm for constructing symmetry-preserving finite element schemes.   In Section <ref>, our constructions are illustrated with several examples of ordinary differential equations. Numerical results are presented that verify the convergence of the proposed invariant finite element schemes and show that symmetry-preserving finite element schemes can provide better numerical results than their non-invariant counterparts.  In Section <ref> we explain how to adapt the constructions introduced for ordinary differential equations to (1+1)-dimensional evolutionary partial differential equations using the method of lines. This is illustrated using Burgers' equation as an example. Finally, in Section <ref> we summarize our findings and give some directions for future research.




§ STATEMENT OF THE PROBLEM



Let x ∈ℝ be the independent variable, and u=u(x) a real-valued scalar function. In the following we consider single second-order ordinary differential equations written in the form


    u_xx = Δ(x,u,u_x).


Here and in what follows, we are using the index notation for derivatives, i.e. 

    u_x= u/ x    and     u_xx=^2 u/ x^2.
 
Now, let G be an r-parameter Lie group acting locally on the plane ℝ^2 parametrized by (x,u).  Using capital letters to denote the transformed variables, we have 


    X = g· x    and     U = g· u,  where   g∈ G.


The group action (<ref>) induces a prolonged action on the derivatives given by the chain rule:

    U_X = D_x(U)/D_x(X),     U_XX = D_x(U_X)/D_x (X),

where

    D_x = x + u_xu + u_xxu_x + ⋯

is the total derivative operator with respect to the independent variable x.


A local Lie group of transformations G acting on an open subset of ℝ^2 is said to be a symmetry group of the differential equation (<ref>) if the solution space of the equation in invariant under the given group action.  In other words,


    U_XX = Δ(X,U,U_X)    whenever     u_xx = Δ(x,u,u_x).




The main goal of this paper consists of recasting (<ref>) into its weak form, and to introduce a systematic procedure for constructing a discrete approximation of the weak form that will preserve the symmetries of the original differential equation. To achieve this goal, we introduce the space of real-valued locally integrable functions on ℝ,

    L_1,(ℝ) = {fℝ→ℝ | ∫_K |f| dx < ∞ for all compact subsets  K ⊂ℝ}.

Alternatively, L_1,(ℝ) is defined as the set of functions fℝ→ℝ such that for any compactly supported test function ϕ∈ C_c^∞(ℝ), the integral

    ∫_-∞^∞ |fϕ| dx < ∞

is finite.  We now assume that solutions to (<ref>) and their derivatives are in L_1,(ℝ).  Multiplying the differential equation (<ref>) by a test function ϕ∈ C_c^∞(ℝ), and integrating over ℝ, we obtain, using integration by parts, the weak formulation of equation (<ref>):


    0 = ∫_∞^∞ [-u_xx + Δ(x,u,u_x)]ϕ dx = ∫_-∞^∞ [u_xϕ_x + Δ(x,u,u_x)ϕ] dx.


Let G be the symmetry group of the differential equation (<ref>).  The group G acts on the test function ϕ via the usual group action on functions:

    Φ = g·ϕ = ϕ(g· x) = ϕ(X),     g∈ G.

The induced action on ϕ_x is given by the chain rule

    Φ_X = ϕ_x/D_x X.

The group G also acts on the differential dx, <cit.>.  The action is given by


    ω = g·dx =  (D_xX) dx.


Restricting our attention to local Lie group actions, <cit.>, we assume that g∈ G is near the identity element so that the bounds of integration in (<ref>) remain infinite once an element of the symmetry group acts on the weak form.  The following theorem is essential for our consideration of invariant finite element discretizations.


The invariance of the differential equation (<ref>) implies the invariance of the weak form (<ref>). 



Indeed, 


    ∫_-∞^∞ [U_XΦ_X + Δ(X,U,U_X)Φ] ω 
     = ∫_-∞^∞ [-U_XX + Δ(X,U,U_X)]Φ ω 
     = ∫_-∞^∞ 0·Φ ω = 0,


where (<ref>) was used.


This theorem is essential as it guarantees that for a given differential equation with symmetry group G, the Lie group G remains a symmetry group of its corresponding weak form.  Therefore, when seeking to construct a symmetry-preserving numerical scheme for a particular differential equation, one can either start with the original strong form or work with a suitable weak form. The strong form of a differential equation is the starting point for constructing symmetry-preserving finite difference schemes, which is the route that has been taken so far in the literature, <cit.>. On the other hand, the weak form is the starting point for constructing symmetry-preserving finite element schemes, which is the focus of the present paper.

To approximate (<ref>), we subdivide the real line ℝ into the elements [x_n, x_n+1]. For an introduction to the theory of finite elements, we refer the reader to <cit.>.  In this paper, the space of test functions C_c^∞(ℝ) is replaced by the space of hat functions

    H^ = {ϕ_kℝ→ℝ |  k∈ℤ},

where


0.45


[scale = 0.75]
[
 axis lines=middle,
 xmin=-0.1,xmax=2.5,
 xtick=0.25,1.25,2.25,xticklabels=x_k-1,x_k,x_k+1,
 ymin=-0.1,ymax=1.1,
 ytick=1,yticklabels=1,
 samples=200]
[domain=0.25:1.25,thick] x-0.25;
[domain=1.25:2.25,thick] 2.25-x;

at (3.75,5.5) ϕ_k(x);



0.55

    ϕ_k(x) = x-x_k-1x_k - x_k-1    x ∈ [x_k-1,x_k] 
    x_k+1 - xx_k+1 - x_k    x ∈ [x_k,x_k+1] 
     0.75cm 0     x ∉ [x_k-1,x_k+1]
    .



The solution u(x) to the weak formulation (<ref>) is now approximated by the (infinite) linear combination 


    u(x) ≈ u^(x) = ∑_k=-∞^∞ u_k ϕ_k(x),


where u_k = u(x_k) denotes the value of the function u(x) at the node x_k.  A first order approximation of the weak form (<ref>) is then given by


    0 = ∫_-∞^∞ [u_x^ϕ^'_k + Δ(x,u^,u^_x)ϕ_k] dx,


where ϕ_k^' = D_x(ϕ_k) denotes the derivative of ϕ_k, and 

    u_x^(x) = ∑_k=-∞^∞ u_k ϕ_k^'(x)

approximates the first derivative u_x.

The transformation group (<ref>) induces an action on the discrete weak form (<ref>).  The action on the nodes x_k and the coefficients u_k in the expansion (<ref>) is given by the product action

    X_k = g· x_k,     U_k = g· u_k,     k∈ℤ,

and the action on the hat function ϕ_k(x) is

    Φ_k = g·ϕ_k = ϕ_k(X) = X-X_k-1X_k - X_k-1    X ∈ [X_k-1,X_k] 
    X_k+1 - XX_k+1 - X_k    X ∈ [X_k,X_k+1] 
     0.85cm 0     X ∉ [X_k-1,X_k+1]
    .

We then introduce the transformed interpolated function


    U^(X):= ∑_k=-∞^∞ U_kΦ_k(X).


In the subsequent developments, we require the transformed approximation (<ref>) to be of the same form as the original interpolant (<ref>).  In other words, we require U^ to be a linear combination of basis functions that depend solely on the independent variable x.  This can be achieved by requiring that Φ_k(X) is a function of x only and not of u.  To do so, we require the group action to be projectable, <cit.>. This assumption requires the transformation rule in the independent variable to be a function of x (and the group parameters):

    X = g· x = X(x,g),     U=g· u = U(x,u,g).


The reason for requiring the group action to be projectable comes from the fact that if this were not the case, then the transformed hat function Φ_k would depend on the unknown function u(x), which would make it impossible to evaluate the integral in (<ref>), and therefore make it impossible to obtain the corresponding finite element scheme.  Investigating the possibility of extending the constructions to general, non-projectable group actions remains to be done.  We do stress here though that most symmetry groups of differential equations arising as models in the mathematical sciences are indeed projectable, and thus the projectability assumption captures essentially all equations of practical relevance.

With this in mind, we now recall a theorem due to Lie, <cit.>.


The largest Lie group contained in the diffeomorphism pseudo-group of the real line 𝒟(ℝ) is the special linear group SL(2,ℝ).  Up to a local diffeomorphism, the action of SL(2,ℝ) on the real line ℝ is given by fractional linear transformations:


    X = α x + β/γ x + δ,    αδ - βγ = 1.




Under the fractional linear transformation (<ref>) the hat function ϕ_k(x) transforms according to the formula



    Φ_k(X) = ϕ_k(x) ·γ x_k + δ/γ x + δ,


while the transformation rule for the first derivative is 


    Φ_k^'(X) = (γ x_k+δ)[(γ x+δ) ϕ_k^'(x) - γϕ_k(x)],

where ϕ_k(x) is differentiable.




Formula (<ref>) is obtained by substituting (<ref>) into the definition of the hat function ϕ_k(x) in (<ref>).  As for (<ref>), the chain rule yields


    Φ_k^'(X)    = 1/D_x XD_x[Φ_k(X)] = (γ x+ δ)^2D_x[ϕ_k(x)·γ x_k + δ/γ x + δ]
       = (γ x_k+δ)[(γ x+δ) ϕ_k^'(x) - γϕ_k(x)].



Under the action (<ref>), formula (<ref>)  for the induced action on the differential dx becomes


    ω = (D_xX) dx = dx/(γ x + δ)^2.


Knowing how each constituent of the discrete weak form (<ref>) transforms under the action of the Lie group G, we can now address the main purpose of the paper.  Given a second-order ordinary differential equation with symmetry group G and weak form (<ref>),  we seek to construct, in a systematic fashion, a weak form approximation that will remain invariant under the symmetry group of the differential equation.  In general, the naive discretization (<ref>) will not preserve all the symmetries of the continuous problem.  To construct a symmetry-preserving discrete weak form we will use the method of equivariant moving frames, <cit.>, which is endowed with an invariantization map that can be used to map non-invariant quantities to their invariant counterparts.  In our case, we will use the invariantization map to invariantize the discrete weak form (<ref>), resulting in a symmetry-preserving finite element scheme.


As a simple example, consider the second-order linear ordinary differential equation


    u_xx + p(x)u_x+q(x)u=f(x),


where p, q, and f are arbitrary smooth functions of their argument. Equation (<ref>) admits a two-parameter symmetry group given by


    X=x,     U = u+ϵ_1α(x) + ϵ_2γ(x),


where α(x) and γ(x) are two linearly independent solutions of the homogeneous equation u_xx + p(x)u_x+q(x)u=0.  The corresponding weak form of (<ref>) is



    ∫_-∞^∞ [-u_xϕ_x + (p(x)u_x+q(x)u-f(x))ϕ] dx = 0,



while an approximation to this weak form is given by


    ∫_-∞^∞ [-u_x^ϕ_k^' + (p(x) u_x^ + q(x)u^ - f(x))ϕ_k] dx = 0.


Acting on the latter with the group action (<ref>), we obtain


    0    = ∫_-∞^∞ [-u_x^ϕ_k^' + (p(x) u_x^ + q(x)u^ - f(x))ϕ_k] dx 
        1cm + ϵ_1 ∫_-∞^∞ [-α_x^ϕ_k^' + (p(x) α_x^ + q(x)α^)ϕ_k] dx 
        1cm + ϵ_2 ∫_-∞^∞ [-γ_x^ϕ_k^' + (p(x) γ_x^ + q(x)γ^)ϕ_k] dx,


where

    α^ = ∑_k=-∞^∞α_k ϕ_k(x),    α^_x = ∑_k=-∞^∞α_k ϕ_k^'(x),    α_k = α(x_k),

and similarly for γ^ and γ^_x.  Since the last two integrals in (<ref>) are, in general, nonzero, the discrete weak form (<ref>) does not admit the superposition principle given by (<ref>).



As a less trivial example, consider the  second-order nonlinear ordinary differential equation


    u_xx = 1/u^3.


This equation is invariant under the group action


    X = α x + β/γ x + δ,     U = u/γ x + δ,    αδ - βγ = 1,


and a weak formulation of (<ref>) is given by


    0 = ∫_-∞^∞[u_x ϕ_x + 1/u^3ϕ] dx.


Approximating u^-3 by

    1/u^3≈∑_k=-∞^∞1/u_k^3ϕ_k(x),

we obtain the discrete weak form


    0 = ∫_-∞^∞[∑_ℓ = -∞^∞(u_ℓϕ_ℓ^'ϕ_k^' + 1/u_ℓ^3ϕ_ℓϕ_k)] dx.


Acting on (<ref>) with the symmetry group (<ref>), recalling (<ref>) and (<ref>), we obtain, after simplification,


    0=∫_-∞^∞[∑_ℓ = -∞^∞( u_ℓϕ_ℓ^'ϕ_k^' + 1/u_ℓ^3(γ x_ℓ + δ/γ x + δ)^4 ϕ_ℓϕ_k)] dx.


The extra factor (γ x_ℓ + δγ x + δ)^4 in the second term of (<ref>) shows that the discrete weak form (<ref>) is not invariant under the group action (<ref>).


We conclude this section by observing that all our considerations can be restricted to boundary value problems, which are more standard in the application of the finite element method.  Instead of working on the whole real line ℝ, simply restrict all considerations to an interval [a,b] and impose boundary conditions at x=a and x=b.  The symmetry group G should now consist of all (or a subset of all) transformations that keep the differential equation and its boundary conditions invariant.  As the boundary conditions impose further constraints, the symmetry group of the boundary value problem will usually be smaller than the symmetry group of the differential equation itself, <cit.>.  A slightly less restrictive assumption is to allow symmetry transformations of a given system of differential equations without boundary conditions to act as equivalence transformations preserving a class of boundary value problems containing the problem under consideration, <cit.>.




§ MOVING FRAMES



The theoretical foundations of the discrete equivariant moving frame method have recently been developed in <cit.>. For the sake of completeness of the present exposition, we summarize the theory of moving frames relevant to the construction of symmetry-preserving finite element schemes here. 

After evaluating the discrete weak form (<ref>), the result is a function of the discrete points (x_k-1,u_k-1), (x_k,u_k), and (x_k+1,u_k+1).  In the following, we combine these three points into the second-order discrete jet at k ∈ℤ:

    z^[2]_k = (k,x_k-1,u_k-1,x_k,u_k,x_k+1,u_k+1).

The terminology stems from the fact that z^[2]_k contains sufficiently many points to approximate the function u(x) and its derivatives u_x, u_xx at the node x_k using central differences.  We introduce the second-order discrete jet space 

    ^[2] = ⋃_k=-∞^∞  z_k^[2],

which consists of the union of the second-order discrete jets over the integers k ∈ℤ.  The discrete jet space ^[2] admits the structure of a lattice variety or lattifold, which is a manifold-like object modeled on ℤ rather than ℝ, <cit.>.  Alternatively, ^[2] is a disconnected manifold with fibers isomorphic to the Euclidean space ℝ^6.  In the following, we let π^[2]→ℤ denote the projection onto the discrete index k:

    π(z^[2]_k) = k.


Now, let G be an r-parameter Lie group acting on the plane ℝ^2={(x,u)}.  Extending the action trivially to ℤ,

    g· k = k,

the Lie group G induces an action on the discrete jet z^[2]_k via the product action


    Z^[2]_k    = g· z^[2]_k 
    
    (k,X_k-1,U_k-1,X_k,U_k,X_k+1,U_k+1)     = 
    (k,g· x_k-1,g· u_k-1, g· x_k, g· u_k, g· x_k+1, g· u_k+1).

See <cit.> for further details. In other words, the Lie group G induces an action on each fiber of ^[2] via the product action.  In the following, we assume that the action is (locally) free and regular on each fiber π^-1(k)=^[2]|_k.  This forces G ⩽^[2]|_k = 6.  We recall that the product action is free at z_k^[2] if the isotropy group

    G_z_k^[2] = { g ∈ G |  g· z_k^[2] = z_k^[2]} = {e}

is trivial, and that the action is locally free at z_k^[2] if the isotropy group is discrete.  On the other hand, the action is regular if the group orbits have the same dimension and each point in ^[2]|_k has arbitrarily small neighborhoods whose intersection with each orbit is a connected subset thereof.


Let G act (locally) freely and regularly on (each fiber of) ^[2] by the product action (<ref>).  A discrete (right) moving frame is a G-equivariant map ρ^[2]→ G satisfying 

    ρ(g· z_k^[2]) = ρ(z_k^[2])  g^-1

for all g ∈ G where the product action is defined.


The construction of a discrete moving frame is based on the introduction of a (collection of) cross-section(s) 𝒦⊂^[2] to the group orbits.


A subset 𝒦⊂^[2] is a cross-section to the group orbits if for each k ∈ℤ, the restriction 𝒦|_k ⊂^[2] = π^-1(k) is a submanifold of ^[2]|_k, transverse and of complementary dimension to the group orbits.


In general, a cross-section 𝒦⊂ J^[2] is specified by a system of r= G difference equations



    𝒦 = { E_ℓ(z_k^[2]) = 0  | ℓ = 1,…,r}.



The right moving frame ρ(z_k^[2]) at z_k^[2] is then the unique group element in G that sends z_k^[2] onto 𝒦|_k:

    ρ(z_k^[2])· z_k^[2]∈𝒦|_k.

The coordinate expressions of the moving frame are obtained by solving the normalization equations

    E_ℓ(g· z_k^[2]) = 0,    ℓ = 1,…,r,

for the group parameters g=(g_1,…,g_r).

Given a moving frame, there is a systematic procedure for constructing invariant functions, invariant differential forms, and other invariant quantities, <cit.>.  


Let ρ^[2]→ G be a right moving frame.  The invariantization of the difference function F(k,x_i,u_i,…,x_j,u_j) is the invariant


    ι_k(F)(k,x_i,u_i,…,x_j,u_j) = F(k,ρ_k· x_i,ρ_k· u_i,…,ρ_k· x_j,ρ_k· u_j)


obtained by acting on the arguments of F with the moving frame ρ_k=ρ(z_k^[2]). 


Borrowing the notation from <cit.>, we can rewrite (<ref>) as

    ι_k(F)(k,x_i,u_i,…,x_j,u_j) = F(k,g· x_i,g· u_i,… g· x_j, g· u_j)|_g=ρ_k.

Thus, the invariantization of F is obtained by first acting on its argument by an arbitrary group element g ∈ G, followed by the substitution g=ρ_k.  In particular, the invariantization of the components of a point (x_ℓ,u_ℓ) are the invariants

    ι_k(x_ℓ) = g· x_ℓ|_g=ρ_k,    ι_k(u_ℓ) = g· u_ℓ|_g=ρ_k.

Similarly, we can also invariantize the hat function ϕ_ℓ(x) and its derivative:

    ι_k(ϕ_ℓ) = (g·ϕ_ℓ)|_g=ρ_k,    ι_k(ϕ_ℓ^') = (g·ϕ_ℓ^')|_g=ρ_k.

Therefore, the invariantization of u^(x) and u^_x is

    ι_k(u^) = ∑_ℓ = -∞^∞ι_k(u_ℓ)ι_k(ϕ_ℓ),    ι_k(u^_x) = ∑_ℓ = -∞^∞ι_k(u_ℓ)ι_k(ϕ^'_ℓ).

Finally, according to (<ref>), the invariantization of the one-form dx is the invariant one-form

    ι_k(dx) = (D_x X)|_g=ρ_k dx.

The one-form ι_k(dx) is invariant since we limit our considerations to projectable group actions.  For general group actions, ι_k(dx) would be contact-invariant[A differential form Ω on the jet space ^(n) is said to be contact-invariant if and only if, for every g ∈ G, g^*Ω = Ω + θ_g for some contact form θ_g, <cit.>.].



As an example of the moving frame construction introduced above, let us consider the group action (<ref>).  Introducing the centered difference derivative

    u_x^k = u_k+1-u_k-1/x_k+1-x_k-1,

a cross-section on ^[2] is given by 

    𝒦 = {x_k = 0,  u_k = 1, u_x^k = 0 }.
 
Solving the normalization equations

    0=X_k = α x_k + β/γ x_k + δ,    
    1=U_k = u_k/γ x_k + δ,    
    0=U_X^k = γ(x_k u_x^k - u_k)+δ u_x^k,

for the group parameters, where 

    x_k = x_k+1+x_k-1/2,    u_k = u_k+1+u_k-1/2,

we obtain the discrete moving frame


    α = 1/u_k,    β = -x_k/u_k,    γ = u_k u_x^k/(x_k-x_k)u_x^k + u_k,    δ = u_k[u_k - x_k u_x^k]/(x_k - x_k)u_x^k + u_k.


Invariantizing u_ℓ, we obtain the invariant


    ι_k(u_ℓ) = u_ℓ/γ x_ℓ + δ|_(<ref>) = u_ℓ[(x_k - x_k)u_x^k + u_k]/u_k[(x_ℓ - x_k)u_x^k + u_k],


while the invariantization of the hat function ϕ_ℓ(x) is


    ι_k(ϕ_ℓ(x)) = ϕ_ℓ·γ x_ℓ + δ/γ x +δ|_(<ref>) = ϕ_ℓ·(x_ℓ-x_k)u_x^k + u_k/(x-x_k)u_x^k + u_k.


Combining (<ref>) and (<ref>), we obtain the invariantization of u^:

    ι_k(u^) =  ∑_ℓ=-∞^∞u_ℓ/u_k·(x_k-x_k)u_x^k + u_k/(x-x_k)u_x^k + u_k ϕ_ℓ(x).

Finally, the invariantization of the one-form dx is

    ι_k(dx) = dx/(γ x + δ)^2|_(<ref>) = [(x_k - x_k)u_x^k + u_k]^2/u_k^2[(x-x_k)u_x^k + u_k]^2 dx.







§ SYMMETRY-PRESERVING FINITE ELEMENT SCHEMES



Given a second-order ordinary differential equation of the form (<ref>) with projectable symmetry group G, we now have everything in hand to construct a symmetry-preserving finite element scheme.  First, rewrite the differential equation in its weak form (<ref>).  Then, consider the discrete approximation (<ref>) or any other suitable approximation.  In general, the discrete weak form will not preserve all the symmetries of the differential equation.  To obtain a symmetry-preserving finite element scheme, first construct a discrete moving frame for the symmetry group G as explained in Section <ref>.  Then use the corresponding invariantization map to invariantize the discrete weak form (<ref>) or any suitable approximation.

To guarantee the consistency of the symmetry-preserving finite element scheme, we need to impose certain constraints on the general moving frame constructions introduced in Section <ref>.  Namely, in the continuous limit where the lengths of the elements [x_k-1,x_k] and [x_k,x_k+1] go to zero, all discrete constructions need to converge to their continuous counterparts.  To guarantee this convergence, we have to construct a consistent moving frame compatible with a differential moving frame, <cit.>.  In other words, the discrete moving frame should, in the continuous limit, converge to a moving frames defined for the prolonged action of G on the submanifold jet ^(2)={(x,u,u_x,u_xx)}, <cit.>.  This will be the case if the cross-section 𝒦 used to define the discrete moving frame converges, in the continuous limit, to a cross-section in ^(2).  In practice, this can be accomplished by using solely, for example, the approximations


    x_k,     u_k,     u_x^k = u_k+1-u_k-1/x_k+1-x_k-1,    
    u_xx^k = 2/x_k+1-x_k-1(u_x^k - u_x^k-1)


to define a discrete cross-section as in the continuous limit those quantities converge to x, u, u_x, and u_xx, respectively.



 §.§ Ordinary Differential Equations


In this section we consider several ordinary differential equations to illustrate the construction of symmetry-preserving finite element discretizations. 

We note that all the schemes presented below are implicit and hence require the solution of a (nonlinear) algebraic equation. For this purpose, we used Newton's method with a termination tolerance of 10^-15 in all numerical examples.



  §.§.§ Equation u_xx=exp(-u_x)


As our first example, we consider the equation


    u_xx = exp(-u_x).


This equation admits the three-parameter symmetry group action


    X= e^ϵ x + a,     U = e^ϵ u + ϵ e^ϵ x + b,    ϵ, a, b ∈ℝ.


A weak form formulation of equation (<ref>) is given by


    0 = ∫_-∞^∞(u_x ϕ_x + e^-u_xϕ) dx.


An approximation of (<ref>) is provided by


    0 = ∫_-∞^∞(u^_x ϕ_k^' + e^-u^_xϕ_k) dx.


We now show that the discrete weak form (<ref>) is already invariant under the group action (<ref>).  First, we have

    U_X^ = u_x^ + ϵ,    Φ_k = ϕ_k,    Φ_k^' = 1/e^ϵϕ_k^',    ω = e^ϵ dx.

Therefore

    0    = ∫_-∞^∞ [U_X^Φ_k^' + e^-U_X^Φ_k] ω
    = ∫_-∞^∞[(u_x^ + ϵ)ϕ_k^'/e^ϵ + e^-(u_x^ + ϵ)ϕ_k] e^ϵdx
       = ∫_-∞^∞ [ u_x^ϕ_k^' + e^-u_x^ϕ_k] dx + ϵ∫_-∞^∞ϕ_k^'  x = ∫_-∞^∞ [ u_x^ϕ_k^' + e^-u_x^ϕ_k] dx,


since ∫_-∞^∞ϕ_k^'  x = 0.  Evaluating the integral (<ref>) we obtain the symmetry-preserving finite element scheme


    (Δ x_k+Δ x_k-1) u_xx^k = Δ x_k-1exp[-u_x^k-1] + Δ x_k exp[-u_x^k],


where 

    Δ x_k = x_k+1-x_k,     
    u_x^k = u_k+1-u_k/x_k+1-x_k,    
    u_xx^k = 2/x_k+1-x_k-1[u_x^k - u_x^k-1].

We observe that the finite element scheme (<ref>) differs from the two schemes appearing in <cit.> (equations (4.25) and (4.26)).

We now test the invariant scheme (<ref>) numerically,  by treating equation (<ref>) as an initial value problem. First, we note that the exact solution to the equation (<ref>) is

    u_ a(x)=(x+c_1)ln(x+c_1)-x+c_2,

where c_1 and c_2 are two arbitrary constants.  Using the initial conditions u(0)=1 and u_x(0)=0, the exact solution becomes u_ a(x)=(x+1)ln(x+1)-x+1. Integrating (<ref>)  from x=0 to x=1, the convergence plot of the relative l_∞-error is shown in Figure <ref>. As it can be seen, the scheme converges at first order, in accordance with the derivation of the finite element scheme,  which is based on a first order linear interpolant.






  §.§.§ Equation u_xx+p(x)u_x+q(x)u=f(x)


In Example <ref>, we observed that the discrete weak formulation (<ref>) does not preserve the linear superposition principle (<ref>) for the linear equation (<ref>).  To solve this problem, we now construct a symmetry-preserving finite element scheme.

The first step is to construct a moving frame.  As in (<ref>), let

    α_x^k = α_k+1-α_k-1/x_k+1-x_k-1,    γ_x^k = γ_k+1 - γ_k-1/x_k+1-x_k-1

denote the centered first derivative approximations.  In the following, we assume that

    γ_k = γ(x_k) ≠ 0,    γ_kα_x^k - α_k γ_x^k≠ 0.

The second constraint is a discrete approximation of the Wronskian condition requiring that the solutions α(x) and γ(x) are linearly independent.  We construct a moving frame by choosing the cross-section

    𝒦 = {u_k = u_x^k = 0},

where u_x^k is the centered approximation introduced in (<ref>).
Solving the normalization equations

    0 = U_k =  u_k + ϵ_1α_k + ϵ_2γ_k,    
    0 = U^k_X = u_x^k + ϵ_1 α_x^k + ϵ_2 γ^k_x,

for the group parameters ϵ_1 and ϵ_2, we obtain


    ϵ_1 = u_k γ_x^k - γ_k u_x^k/γ_k α_x^k - α_k γ_x^k,    
    ϵ_2 = α_k u_x^k - u_k α_x^k/γ_k α_x^k - α_k γ_x^k.


Given the moving frame (<ref>), we invariantize the non-invariant discrete weak form (<ref>).  This is done by substituting the group normalizations (<ref>) into (<ref>).  The result is the invariant discrete weak form


    0    = ∫_-∞^∞ [-u_x^ϕ_k^' + (p(x) u_x^ + q(x)u^ - f(x))ϕ_k] dx 
        1cm + u_k γ_x^k - γ_k u_x^k/γ_k α_x^k - α_k γ_x^k∫_-∞^∞ [-α_x^ϕ_k^' + (p(x) α_x^ + q(x)α^)ϕ_k] dx 
        1cm + α_k u_x^k - u_k α_x^k/γ_k α_x^k - α_k γ_x^k∫_-∞^∞ [-γ_x^ϕ_k^' + (p(x) γ_x^ + q(x)γ^)ϕ_k] dx.


For second-order linear homogeneous equations, i.e. when f(x)=0 in (<ref>), we notice that


    u^(x) = c_1 α^(x) + c_2 γ^(x),


where c_1 and c_2 are two arbitrary constants, is an exact solution of the discrete weak form (<ref>).  Indeed, when u^(x) is given by (<ref>), we have that

    u_k γ_x^k - γ_k u_x^k/γ_k α_x^k - α_k γ_x^k =  -c_1,    α_k u_x^k - u_k α_x^k/γ_k α_x^k - α_k γ_x^k = -c_2,

and the right-hand side of (<ref>) is identically zero.





 







  §.§.§ Equation u_xx=u^-3


As a third example, we consider the nonlinear differential equation (<ref>), with discrete weak form (<ref>).  In Example <ref>, we computed a moving frame for the symmetry group (<ref>).  The result is given in equation (<ref>).  Invariantizing the discrete weak form (<ref>), which is obtained by substituting the group parameter normalizations (<ref>) into the transformed discrete weak form (<ref>), we get the symmetry-preserving discrete weak form

    0 = ∫_-∞^∞[∑_ℓ=-∞^∞(u_ℓϕ_ℓ^'ϕ_k^' + 1/u_ℓ^3((x_ℓ - x_k)u_x^k + u_k/(x-x_k)u_x^k + u_k)^4ϕ_ℓϕ_k)] dx.

Integrating this expression yields the symmetry-preserving finite element scheme


    -(u_k+1-u_k/x_k+1-x_k) + (u_k-u_k-1/x_k - x_k-1) + 
    (x_k - x_k-1)[(x_k-1-x_k)u_x^k + u_k]^2/6u_k-1^3[(x_k-x_k)u_x^k + u_k]^2
    
    +(x_k+1-x_k-1)[(x_k-x_k)u_x^k + u_k]^2/3u_k^3[(x_k+1-x_k)u_x^k + u_k][(x_k-1-x_k)u_x^k + u_k]
    +(x_k+1-x_k)[(x_k+1-x_k)u_x^k + u_k]^2/6u_k+1^3[(x_k-x_k)u_x^k + u_k]^2=0.


We now turn to the numerical verification of the resulting invariant scheme. First, we note that the general solution to the differential equation (<ref>) is

    u_ a^2(x)=1/c_1+c_1(x+c_2)^2,

where c_1 and c_2 are arbitrary constants with c_10, <cit.>.

We integrate equation (<ref>) on the interval [0,1] using the invariant finite element scheme (<ref>) and the initial conditions u(0)=1, u_x(0)=0. In this case the exact solution reduces to u_ a(x)=√(1+x^2).   The convergence plot for the scheme (<ref>) is presented in Figure <ref>. As expected, this invariant scheme converges at first order, since it is based on a linear interpolant.




As an alternative discretization of the weak form (<ref>), one could also consider the approximation



    ∫_-∞^∞ [u^_xϕ_k^' + (u^)^-3ϕ_k] dx = 0.



Under the group action (<ref>), this weak form gets dilated to

    0 = ∫_-∞^∞[U_X^Φ_k^' + (U^)^-3Φ_k] ω 
    = (γ x_k + δ) ∫_-∞^∞ [u_x^ϕ_k^' + (u^)^-3ϕ_k] dx,

and therefore preserves the symmetries of the original weak form (<ref>).  The corresponding symmetry-preserving finite element scheme is



    u_k+1-u_k/x_k+1-x_k - u_k - u_k-1/x_k - x_k-1 = x_k+1-x_k/2u_k^2 u_k+1 + x_k-x_k-1/2u_k^2 u_k-1.






  §.§.§ Painlevé Equation u_xx = u^-1u_x^2


As our final example we consider the Painlevé equation


    u_xx = u_x^2/u.


This equation admits a six-parameter symmetry group of projectable transformations given by

    X = α x + β/γ x + δ,    
    U = (u^λ e^ax+b)^1/(γ x + δ),

where αδ - βγ = 1, a, b ∈ℝ, and λ >0.  In the following, we restrict our attention to the two-dimensional symmetry group


    X=x,     U = u e^ax+b,


A weak formulation of the Painlevé equation (<ref>) is given by

    0 = ∫_-∞^∞[ u_x ϕ_x + u_x^2/uϕ]  x.

In the discrete setting, we approximate the weak form by


    0=∫_-∞^∞[ u_x^ϕ_k^' + (u_x^)^2/u^ϕ_k] x.


Integrating (<ref>), we obtain the non-invariant finite element scheme


    -2[(u_k+1-u_k/x_k+1-x_k) - (u_k - u_k-1/x_k-x_k-1)] + u_k-1/x_k - x_k-1ln(u_k-1/u_k) + u_k+1/x_k+1-x_kln(u_k+1/u_k) = 0.


To construct a symmetry-preserving finite element scheme, we construct a moving frame to the group action (<ref>) using the cross-section

    𝒦 = { u_k = 1,  u_x^k = u_k+1-u_k-1/x_k+1-x_k-1 = 0}.

Solving the corresponding normalization equations, we obtain

    a = 1/x_k+1 - x_k-1ln(u_k-1/u_k+1),    
    b = x_k/x_k+1-x_k-1ln(u_k+1/u_k-1) - ln u_k.

Invariantizing the discrete weak form (<ref>) and performing the integration we obtain the symmetry-preserving finite element scheme



    -2[(I_k-1/x_k+1-x_k)-(1-J_k/x_k-x_k-1)]+ J_kln J_k/x_k-x_k-1+I_kln I_k/x_k+1-x_k = 0,


where the invariants I_k and J_k are given by


    I_k = u_k+1/u_kexp[-x_k+1-x_k/x_k+1-x_k-1ln(u_k+1/u_k-1)],    
    J_k = u_k-1/u_kexp[x_k-x_k-1/x_k+1-x_k-1ln(u_k+1/u_k-1)].



We now compare the invariant scheme (<ref>) against the non-invariant scheme (<ref>) numerically. Since the symmetry-preserving scheme is exact, i.e. the only difference between the numerical solution and the exact solution is due to round-off error, we do not need to verify the convergence of the scheme. A straightforward Taylor series analysis reveals that the non-invariant finite element scheme for the Painlevé equation is of second order.

We now solve the initial value problem for the Painlevé equation with initial conditions u(0)=1 and u_x(0)=1, corresponding to the exact solution u=exp(x).  Integrating over the interval [0,1] using a step size of Δ x=0.01, the time series of the relative error between the numerical solutions of the two schemes (<ref>), (<ref>) and the exact solution is depicted in Figure <ref>.  It is obvious that the symmetry-preserving scheme outperforms the non-invariant scheme, with the error of the invariant scheme being several magnitudes smaller and approximately of the size of machine epsilon.




We note that numerically solving the nonlinear algebraic equation (<ref>) for the invariant finite element method is challenging due to the fact that this scheme is exact. Numerically, we observe an accumulation of round-off errors that is growing over the integration interval. The smaller the step size Δ x, the more round-off error can accumulate. To numerically preserve the exactness of the scheme for all step sizes Δ x, variable precision arithmetic may be necessary.





 §.§ Partial Differential Equations



In this section we extend the constructions introduced in the previous sections to the semi-discretization of (1+1)-dimensional evolution equations, where only the spatial variable is discretized.  This allows us to use many of the ideas introduced in the previous sections. To simplify the exposition, we focus on a particular example and consider Burgers' equation


    u_t + u u_x = ν u_xx,    where    ν >0,


which plays an important role in various areas of applied mathematics, such as fluid mechanics, nonlinear acoustics, gas dynamics, and traffic flow. Here ν is the constant viscosity coefficient. Burgers' equation admits a five-parameter maximal Lie symmetry group, see e.g. <cit.>.  One of these admitted symmetry transformations yields an inversion of time, which does not respect the requirement that the time variable t should increase monotonically for a given initial value problem, <cit.>. Thus, we restrict our attention to the four-parameter  subgroup of symmetry transformations


    X= λ (x + v t) + a,     T = λ^2 t + b,     U = u + v/λ,     a,b,v∈ℝ, λ∈ℝ^+.


Multiplying Burgers' equation (<ref>) by a test function ϕ(x) ∈ C_c^∞(ℝ) and integrating over ℝ, we obtain the weak form 


    ∫_-∞^∞ u_t ϕ dx = ∫_-∞^∞(-ν u_x + u^2/2)ϕ_x dx.


In the following, we consider the semi-discretization of Burgers' equation where the spatial variable x is discretized and the time variable t remains a continuous variable. In this setting, the interpolating coefficients in the approximation (<ref>) of the solution now become functions of t:


    u(x,t) ≈ u^(x,t) = ∑_k=-∞^∞ u_k(t) ϕ_k(x).


Substituting (<ref>) into the weak form (<ref>) and replacing the test function by the hat function ϕ_ℓ, we obtain


    ∫_-∞^∞ u_t^ϕ_ℓ dx = ∫_-∞^∞(-ν u_x^ + (u^)^2/2)ϕ^'_ℓ dx,


where

    u_t^ = ∑_k=-∞^∞ u_k/ t (t)ϕ_k(x)    and    
    u_x^ = ∑_k=-∞^∞ u_k(t) ϕ_k^'(x).

Under the group action (<ref>), the differentials dx and dt transform according to


    ω^x = D_x(X) dx + D_t(X) dt = λ(dx + v dt),    ω^t = D_x(T) dx + D_t(T) dt = λ^2 dt,


where D_x and D_t are the total derivative operators in the independent variables x and t, respectively, <cit.>.  Dual to the one-forms (<ref>) are the implicit derivative operators

    D_X = 1/λD_x,    D_T = 1/λ^2(D_t - v D_x).

Therefore, the hat functions and their first derivatives transform according to 

    Φ_ℓ = ϕ_ℓ    and    Φ^'_ℓ = D_X(Φ_ℓ) =  ϕ_ℓ^'/λ.

Finally, we have


    U^   = ∑_k=-∞^∞ U_k(T) Φ_k(X) = ∑_k=-∞^∞u_k+v/λ ϕ_k = u^ + v/λ, 
    
    U^_T    = ∑_k=-∞^∞D_T(U_k)Φ_k = ∑_k=-∞^∞1/λ^3 u_k/ tϕ_k = 1/λ^3  u_t^,
    
    U^_X    = ∑_k=-∞^∞ U_k D_X(Φ_k) = ∑_k=-∞^∞u_k+v/λ^2 ϕ_k^' = 1/λ^2  u_x^,


where we used the fact that ∑_k=-∞^∞ϕ_k = 1 and ∑_k=-∞^∞ϕ_k^' = 0, where the sums are defined.  

We now act on the discrete weak form (<ref>) with the symmetry group (<ref>).  Since the weak form is evaluated at a fixed time, we substitute

    ω^x = λ(dx + v dt) ≡λ dx

into the transformed weak form.  After simplification, we obtain


    ∫_-∞^∞ u_t^ϕ_ℓ  x = ∫[-ν u_x^ + (u^ + v)^2/2]ϕ_ℓ^'  x 
    = ∫_-∞^∞[-ν u_x^ + (u^)^2/2 + v u^]ϕ_ℓ^'  x.


Due to the occurrence of the Galilean boost parameter v, we conclude that the discrete weak form (<ref>) is not invariant under the symmetry subgroup (<ref>).




  §.§.§ Symmetry-Preserving Lagrangian Scheme



In this section we introduce a discrete weak form of Burgers' equation that preserves the symmetry subgroup (<ref>).  This is done by using the Lagrangian form of Burgers' equation given by


    du/dt = ν u_xx,    dx/dt = u,


where

    d/dt = D_t + dx/dtD_x = D_t + u D_x.

In this setting, x is now a function of the time variable t.  Therefore, the nodes x_ℓ are functions of t and the element [x_ℓ,x_ℓ+1] varies as a function of time.

Following the general procedure introduced in the previous sections, the first step in constructing a symmetry-preserving weak form consists of computing a discrete moving frame.  Assuming, for simplicity, that

    u_x^ℓ = u_ℓ+1-u_ℓ-1/x_ℓ+1-x_ℓ-1 > 0    for all    ℓ∈ℤ,

we introduce the cross-section

    𝒦 = { x_ℓ = t = u_ℓ = 0,  u_x^ℓ = 1}.

Solving the normalization equations

    0 = λ (x_ℓ + vt)+a,     0=λ^2t+b,     0 = λ^-1(u_ℓ + v),     1=λ^-2 u_x^ℓ,

for the group parameters, we obtain the moving frame


    a = -√(u_x^ℓ) (x_ℓ - t  u_ℓ),     b = -t  u_x^ℓ,     v = -u_ℓ,    λ = √(u_x^ℓ).


Since u_k and x_k are functions of t, we now wish to invariantize u_k/ t and x_k/ t.  Under the symmetry group action (<ref>), we have


    g· u_k/ t   =  U_k/ T = 1/λ^2/ t[u_k+v/λ] = 1/λ^3 u_k/ t,
    
    g· x_k/ t   =  X_k/ T = 1/λ^2/ t [λ(x_k+vt)+a] = 1/λ( x_k/ t + v),


where

    / T = 1/λ^2/ t

is the derivative operator dual to the one-form ω^t = λ^2 dt.  Using the moving frame (<ref>), we have

    ι_ℓ(du_kdt) = 1/(u_x^ℓ)^3/2du_k/dt,    ι_ℓ(dx_kdt) = 1/√(u_x^ℓ)(dx_k/dt - u_ℓ).

Next, invariantizing the discrete weak form (<ref>), which is obtained by substituting the group normalizations (<ref>) into the transformed discrete weak form (<ref>), we obtain the symmetry-preserving discrete weak form

    ∫_-∞^∞ u_t^ϕ_ℓ  x = ∫_-∞^∞(-ν u_x^ + (u^)^2/2-u_ℓ  u^)ϕ_ℓ^'  x.

Evaluating the integrals, and simplifying the expressions, we obtain the symmetry-preserving finite element scheme



    1/3[x_ℓ-x_ℓ-1/x_ℓ+1-x_ℓ-1·du_ℓ-1/dt + 2 ·du_ℓ/dt + x_ℓ+1-x_ℓ/x_ℓ+1-x_ℓ-1·du_ℓ+1/dt] = ν u_xx^ℓ - (u_ℓ+1-2u_ℓ+u_ℓ-1/3)u_x^ℓ,


where

    u_xx^ℓ = 2/x_ℓ+1-x_ℓ-1[(u_ℓ+1-u_ℓ/x_ℓ+1-x_ℓ) - (u_ℓ-u_ℓ-1/x_ℓ-x_ℓ-1) ].

In the Lagrangian formalism, we need to supplement (<ref>) with a mesh equation that will describe how the node x_ℓ will evolve as a function of time.  This can be achieved, in a symmetry-preserving fashion, by setting ι_ℓ(dx_ℓ/dt)=0, which yields the invariant differential equation


    dx_ℓ/dt = u_ℓ.



In the continuous limit, the invariant scheme (<ref>) converges to (<ref>).


In equation (<ref>) there is no built-in term that would allow to control the evolution of the mesh.  The limit only holds provided the mesh points satisfy equation (<ref>).  This is to be expected as we have invariantized the discrete weak form (<ref>), defined on a fixed mesh together with the mesh equation dx_ℓ/dt = 0, which forces the nodes to stay fixed as the time variable evolves.  


Due to the Lagrangian mesh equation (<ref>), the invariant scheme (<ref>) will in general suffer from mesh tangling and singularities <cit.>.  To avoid these problems, we now construct a symmetry-preserving finite element scheme that will hold on any moving mesh.




  §.§.§ Symmetry-Preserving r-Adaptive Scheme



In this section we construct a symmetry-preserving finite element scheme with a built-in term that takes into account the evolution of the mesh.  This is achieved by invariantizing


    ∫_-∞^∞(u_t^ - u^_x x_t)ϕ_ℓ dx = 
    ∫_-∞^∞(-ν u_x^ + (u^)^2/2)ϕ_ℓ^' dx,


where the extra term on the left-hand side of equation (<ref>) takes into account the movement of the mesh, and where

    x_t = ∑_k=-∞^∞dx_k/dtϕ_k.

In particular, we note that when x_k/ t = 0 for all k, then we recover the discrete weak form (<ref>) of Burgers's equation on a fixed mesh.  

Under the group action (<ref>),

    X_T = ∑_k=-∞^∞ X_k/ TΦ_k = ∑_k=-∞^∞1/λ( x_k/ t + v)ϕ_k = x_t+v/λ.

Therefore, acting by the symmetry group (<ref>) on the discrete weak form (<ref>) we obtain the transformed weak form


    ∫_-∞^∞[u_t^ - u_x^(x_t+v)]ϕ_ℓ dx = ∫_-∞^∞[-ν u_x^ + (u^)^2/2 + v u^]ϕ_ℓ^' dx.


The invariantization of (<ref>) is obtained by substituting the moving frame expressions (<ref>) into (<ref>), which yields the symmetry-preserving discrete weak form

    ∫_-∞^∞[u_t^ - u_x^(x_t-u_ℓ)]ϕ_ℓ dx = ∫_-∞^∞[-ν u_x^ + (u^)^2/2 - u_ℓ  u^]ϕ_ℓ^' dx.







Evaluating the integrals, we obtain the symmetry-preserving finite element scheme


    x_ℓ-x_ℓ-1/3(x_ℓ+1-x_ℓ-1)·du_ℓ-1/dt
    + 2/3·du_ℓ/dt
    + x_ℓ+1-x_ℓ/3(x_ℓ+1-x_ℓ-1)·du_ℓ+1/dt
    -u_ℓ-u_ℓ-1/3(x_ℓ+1-x_ℓ-1)·dx_ℓ-1/dt
    
    - 2/3u_x^ℓ·dx_ℓ/dt 
    - u_ℓ+1-u_ℓ/3(x_ℓ+1-x_ℓ-1)·dx_ℓ+1/dt = ν  u_xx^ℓ - (u_ℓ+1 + u_ℓ+u_ℓ-1/3)u_x^ℓ.


The remaining ingredient is to prescribe dx_ℓ/dt using an invariant mesh equation.  For example, when using the mesh equation (<ref>), the scheme reduces to 

    x_ℓ-x_ℓ-1/3(x_ℓ+1-x_ℓ-1)·du_ℓ-1/dt
    + 2/3·du_ℓ/dt
    + x_ℓ+1-x_ℓ/3(x_ℓ+1-x_ℓ-1)·du_ℓ+1/dt = ν  u_xx^ℓ.

Other invariant mesh equations for Burgers' equation were proposed e.g. in <cit.>, and will not be discussed further here.





§ CONCLUSIONS AND OUTLOOK


In this paper we have, for the first time, laid out a partial theory for constructing symmetry-preserving finite element schemes. This contribution is timely given the large body of literature that exists nowadays regarding the construction of symmetry-preserving finite difference schemes, and due to the obvious importance that finite element discretizations play in mathematical sciences.
 
While we have primarily restricted our attention to second-order differential equations, the principles introduced in this paper are applicable to higher-order differential equations, boundary value problems (though this usually reduces the size of the admitted symmetry group) as well as to multi-dimensional systems of partial differential equations.  A main complication when tackling higher-order differential equations is the necessity to use higher-order basis functions. Conceptually, these higher-order basis functions can readily be included in the theory laid out in the present paper.  Since the resulting computations substantially grow in complexity, we have however abstained from including them here for the sake of clarity of this first exposition on invariant finite element methods.
 
Invariant discretization schemes are a particular class of geometric numerical integrators that are designed to preserve at the discrete level (a subgroup of) the maximal Lie symmetry group of a system of differential equations. The motivation for the development of geometric numerical integrators is that, in general, maintaining the intrinsic geometric properties of a system of differential equations improves the long-term behavior of a numerical integration scheme. In the case of symmetries, it has been shown that invariant integrators play an essential role for blow-up problems, where they have been shown to outperform standard non-invariant integrators. The preservation of symmetries in finite element schemes opens up the possibility to compare invariant finite element schemes against non-invariant finite element discretizations, which has been done for a single example in the present work.  We reserve a more detailed comparison for future work.



 §.§ Acknowledgements


This research was undertaken, in part, thanks to funding from the Canada Research Chairs program, the NSERC Discovery Grant program and the LeverageR&D program of the Research and Development Corporation of Newfoundland and Labrador.  FV would like to thank Memorial University of Newfoundland, where this research was initiated, for the hospitality during his stay.
 

99



BDK-1997
Bakirova, M.I., Dorodnitsyn, V.A., and Kozlov, R.V., Symmetry-preserving difference schemes for some heat transfer equations, J. Phys. A 30 (1997), 8139–8155.

B-2013
Bihlo, A., Invariant meshless discretization schemes, J. Phys. A 46  (2013), 062001, 12 pp.

BCW-2015
Bihlo, A., Coiteux, X., and Winternitz, P., The Korteweg-de Vries equation and its symmetry-preserving discretization, J. Phys. A 48 (2015), 055201, 25 pp.

BN-2013
Bihlo, A., and Nave, J.-C., Invariant discretization scheme using evolution-projection techniques, SIGMA 9 (2013), 052, 23 pp.

BN-2014
Bihlo, A., and Nave,  J.-C., Convecting reference frames and invariant numerical models, J. Comput. Phys.  271 (2014), 656–663.

BP-2012
Bihlo, A., and Popovych, R.O., Invariant discretization schemes for the shallow water equations, SIAM J. Sci. Comput.  34 (2012), B810–B839.

BV-2017
Bihlo, A., and Valiquette, F., Symmetry-preserving numerical schemes, in Symmetries and Integrability of Difference Equations (2017), 261–324, Springer.

BC-2016
Blanes, S., and Casas, F., A Concise Introduction to Geometric Numerical Integration, Monographs and Research Notes in Mathematics, Vol. 23, CRC Press, Boca Raton, 2016.

BA-2002
Bluman, G., and Anco, S., Symmetry and Integration Methods for Differential Equations, Applied Mathematical Sciences, Vol. 154, Springer–Verlag, New York, 2002.

BCW-2006
Bourlioux, A., Cyr-Gagnon, and Winternitz, P., Difference scheme with point symmetries and their numerical tests, J. Phys. A: Math. Gen.  39 (2006), 6877–6906.

BRW-2008
Bourlioux, A., Rebelo, R., and Winternitz, P., Symmetry preserving discretization of SL(2,ℝ) invariant equations, J. Nonlin. Math. Phys.  15 (2008), 362–372.

BS-2007
Brenner, S. and Scott, R., The mathematical theory of finite element methods, Texts in Applied Mathematics, Vol. 15, Springer, New York, 2007.

BD-2001
Budd, C., and Dorodnitsyn, V., Symmetry-adapted moving mesh schemes for the nonlinear Schrödinger equation.  Symmetry and integrability of difference equations, J. Phys. A 34 (2001), 10387–10400.

CRW-2016
Campoamor-Stursberg, R., Rodríguez, M.A., and Winternitz, P., Symmetry preserving discretization of ordinary differential equations.  Large symmetry groups and higher order equations, J. Phys. A: Math. Gen.  49 (2016), 035201, 21 pp.

D-1989
Dorodnitsy, V.A., Transformation groups in net space, J.  Sov. Math.  55 (1991), 1490–1517; translated from Itogi Nauki i Tekhniki, Seriya Sovremennye Problemy Mathematiki, Noveishie Dostizheniya 43 (1989), 149–191. 

DKW-2000
Dorodnitsyn, V., Kozlov, R., and Winternitz, P.,  Lie group classification of second-order ordinary difference equations, J. Math. Phys. 41 (2000), 480–504. 

DW-2000
Dorodnitsyn, V., and Winternitz, P., Lie point symmetry preserving discretization for variable coefficient Korteweg-de Vries equations. Modern group analysis, Nonlinear Dynam.  22 (2000), 49–59.

FO-1999
Fels, M., and Olver, P.J., Moving coframes.  II.  Regularization and theoretical foundations, Acta Appl. Math.  55 (1999), 127–208.

HLW-2006
Hairer, E., Lubich, C., and Wanner, G., Geometric Numerical Integration:  Structure-Preserving Algorithms for Ordinary Differential Equations, Springer Series in Computational Mathematics, Vol. 31, 2nd Edition, Springer, New York, 2006.

HA-1975
Hermann, R., and Ackerman, M., Sophus Lie's 1880 transformation group paper, Math. Sci. Press, Brookline, Massachusetts, 1975.

HLW-2000
Hernándes Heredero, R., Levi, D., and Winternitz, P.,  Symmetry preserving discretization of the Burgers equation.  SIDE III – Symmetries and Integrability of Difference Equations, 197–208, CRM Proc.  Lecture Notes, Vol. 25, Amer. Math. Soc.  Providence, RI, 2000.

HR-2011
Huang, W., and Russell, R.D., Adaptive Moving Mesh Method, Applied Mathematical Sciences, Springer, New York, 2011.

K-2008
Kim, P., Invariantization of the Crank–Nicolson method for Burgers' equation, Phys. D 237 (2008), 243–254.

KO-2004
Kim, P., and Olver, P.J., Geometric integration via multi-space, Regular Chaotic Dyn.  9 (2004), 213–226.

LR-2004
Leimkuhler, B., and Reich, S., Simulating Hamiltonian Dynamics, Cambridge Monographs on Applied and Computational Mathematics, Vol. 14, Cambridge University Press, Cambridge, 2004.

LMW-2015
Levi, D., Martina, L., and Winternitz, P., Structure preserving discretization of the Liouville equation and their numerical tests, SIGMA 11 (2015), 080, 20 pp.

L-1880
Lie, S., Theorie der Transformationsgruppen, Math. Ann.  16 (1880), 441–528; Gesammelte Abhandlungen, Vol. 6, B.G. Teubner, Leipzig, 1927, pp. 1–94, See reference <cit.> for an English translation.

M-2010
Mansfield, E.L., A Practical Guide to the Invariant Calculus, Cambridge University Press, Cambridge, 2010.

MM-2016
Marí Beffa, G., and Mansfield, E.L., Discrete moving frames on lattice varieties and lattice-based multispaces, Found. Comp. Math. 18 (2018), 181–247.

O-1993
Olver, P.J., Applications of Lie Groups to Differential Equations, Second Edition, Vol. 107, Springer, New York, 1993.

O-1995
Olver, P.J., Equivalence, Invariants, and Symmetry, Cambridge University Press, Cambridge, 1995.

O-2001
Olver, P.J., Geometric foundations of numerical algorithms and symmetry, Appl. Alg. Engin. Comp. Commun.  11 (2001), 417–436.

P-1950
Pinney, E., The nonlinear differential equation y”+p(x)y+cy^-3=0, Proc. Amer. Math. Soc. 1 (1950), 681.

QM-2008
Quispel, G.R.W., and McLaren, D.I., A new class of energy-preserving numerical integration methods, J. Phys. A: Math. Theor.  41 (2008), 045206.

RV-2013
Rebelo, R., and Valiquette, F., Symmetry preserving numerical schemes for partial differential equations and their numerical tests, J. Difference Eq. Appl.  19 (2013), 738–757.

SC-1994
Sanz-Serna, J.M., and Calvo, M.P., Numerical Hamiltonian problems,
 Applied Mathematics and Mathematical Computation, Vol. 7, Chapman & Hall,
  London, 1994.

WBN-2016
Wan, A.T.S., Bihlo, A., and Nave, J.C., The multiplier method to construct
  conservative finite difference schemes for ordinary and partial differential
  equations, SIAM J. Numer. Anal. 54 (2016), 86–119.

WBN-2017
Wan, A.T.S., Bihlo, A., and Nave, J.C., Conservative methods for dynamical systems,
  SIAM J. Numer. Anal. 55 (2017), 2255–2285.

ZM-1988
Zhong, G., and Marsden, J.E., Lie–Poisson, Hamilton–Jacobi theory and Lie–Poisson integrators, Phys. Lett. A 133 (1988), 134–139.



