On the development of symmetry-preserving finite element schemes for ordinary differential equations

TL;DR

This paper presents a new finite element method for ordinary differential equations that preserves their Lie point symmetries, leading to improved long-term accuracy and applicability to various symmetry groups and polynomial degrees.

Contribution

It introduces a symmetry-preserving finite element scheme using equivariant moving frames, applicable to diverse ODEs and polynomial degrees, with demonstrated numerical benefits.

Findings

Symmetry-preserving schemes show better long-term accuracy.

Applicable to arbitrary order ODEs and polynomial degrees.

Numerical experiments confirm improved performance.

Abstract

In this paper we introduce a procedure, based on the method of equivariant moving frames, for formulating continuous Galerkin finite element schemes that preserve the Lie point symmetries of initial value problems for ordinary differential equations. Our methodology applies to projectable and non-projectable symmetry group actions, to ordinary differential equations of arbitrary order, and finite element approximations of arbitrary polynomial degree. Several examples are included to illustrate various features of the symmetry-preserving process. We summarise extensive numerical experiments showing that symmetry-preserving finite element schemes may provide better long term accuracy than their non-invariant counterparts and can be implemented on larger elements.

Tables

Table 1. Table 1. The standard finite element approximation ( 18 ) where ( 89 ) and ( 90 ) hold with T = 10 𝑇 10 T=10 .

q	$\tau$	Maximal	$L_2$ error	EOC
		nodal error
	1.56e-01	7.49e-04	1.70e-03	-
0	7.81e-02	1.87e-04	4.25e-04	2.00
	3.91e-02	4.68e-05	1.06e-04	2.00
	1.95e-02	1.17e-05	2.66e-05	2.00
	1.56e-01	3.04e-07	2.19e-05	-
1	7.81e-02	1.90e-08	2.74e-06	3.00
	3.91e-02	1.19e-09	3.43e-07	3.00
	1.95e-02	7.43e-11	4.28e-08	3.00
	1.56e-01	5.31e-11	1.58e-07	-
2	7.81e-02	8.30e-13	9.91e-09	4.00
	3.91e-02	1.39e-14	6.20e-10	4.00
	1.95e-02	4.75e-15	3.87e-11	4.00

Table 2. Table 3. The standard finite element approximation ( 54 ) where ( 93 ), ( 94 ) and ( 95 ) hold with T = 1000 𝑇 1000 T=1000 .

q	$\tau$	Maximal	$L_2$ error	EOC
		nodal error
	1.56e-01	1.45e-01	1.27e-01	-
0	7.81e-02	7.52e-02	3.17e-02	2.00
	3.91e-02	3.83e-02	7.91e-03	2.00
	1.95e-02	1.93e-02	1.98e-03	2.00
	1.56e-01	1.45e-01	7.79e-05	-
1	7.81e-02	7.52e-02	9.81e-06	2.99
	3.91e-02	3.83e-02	1.23e-06	3.00
	1.95e-02	1.93e-02	1.54e-07	3.00
	1.56e-01	1.45e-01	1.48e-06	-
2	7.81e-02	7.52e-02	9.38e-08	3.98
	3.91e-02	3.83e-02	5.88e-09	4.00
	1.95e-02	1.93e-02	3.68e-10	4.00

Table 3. Table 5. The standard finite element approximation ( 65 ) where ( 96 ) and ( 97 ) hold.

q	$\tau$	Maximal	$L_2$ error	EOC
		nodal error
	1.56e-01	2.41e-01	2.48e-02	-
0	7.81e-02	1.35e-01	6.30e-03	1.98
	3.91e-02	7.25e-02	1.58e-03	1.99
	1.95e-02	3.76e-02	3.96e-04	2.00
	1.56e-01	2.34e-01	1.22e-03	-
1	7.81e-02	1.34e-01	1.58e-04	2.94
	3.91e-02	7.23e-02	2.00e-05	2.99
	1.95e-02	3.75e-02	2.50e-06	3.00
	1.56e-01	2.34e-01	6.22e-05	-
2	7.81e-02	1.34e-01	4.11e-06	3.92
	3.91e-02	7.23e-02	2.60e-07	3.98
	1.95e-02	3.75e-02	1.64e-08	3.99

Table 4. Table 7. A table confirming whether the standard finite element approximation ( 73 ) and the invariant approximation ( 76 ) may be successfully solved for various step sizes τ 𝜏 \tau{} when approximating the exact solution ( 99 ) with C = 1 , y 0 = 0.5 formulae-sequence 𝐶 1 subscript 𝑦 0 0.5 C=1,y_{0}=0.5 .

$\tau$	Standard scheme	Invariant scheme
0.390625	✓	✓
0.78125	✓	✓
1.5625	✗	✓
3.125	✗	✓
6.25	✗	✗

Equations

$$\mathcal{F}(t,y,y_t,y_{tt},...,y_{t^{m+1}})=0,$$

$$y(0)=y_0,\dots ,y_{t^m}(0)=y_m.$$

$$y^{(m+1)}=(y,y_t,...,y_{t^{m+1}}),$$

$${\mathcal{F}(g\cdot t,g\cdot y^{(m+1)})|}_{\mathcal{F}(t,y^{(m+1)})=0}=0$$

$$\hat{t}:=g\cdot t,\hat{y}:=g\cdot y,$$

$$\widehat{y_{t^i}}=g\cdot y_{t^i}=\frac{{\mathrm{d}}^i\hat{y}}{\mathrm{d}{\hat{t}}^i}.$$

$$\mathbf{u}=(u_0,u_1,u_2,\dots ,u_m)=(y,y_t,\dots ,y_{t^m}),$$

$$\begin{array}{rl}\mathcal{F}(t,\mathbf{u},{u_m}_t)& =0,\\ {u_0}_t-u_1& =0,\\ ⋮& \\ {u_{m-1}}_t-u_m& =0,\end{array}$$

$$\mathbf{u}(0)=(y_0,y_1,\dots ,y_m).$$

$$y_{tt}-y^{-1}{y}_t^2=0,$$

$$y(0)=y_0,y_t(0)=y_1,$$

$$\hat{t}=\frac{\alpha t+\beta }{\gamma t+\delta },\hat{y}=(y^{\lambda }{e}^{ax+b}{)}^{1/(\gamma x+\delta )},$$

$$\hat{t}=t,\hat{y}=\exp (at+b)y.$$

$$\widehat{y_t}$$

$$\widehat{y_{tt}}$$

$$\widehat{y_{tt}}-(\widehat{y})^{-1}(\widehat{y_{t}})^{2}=\big{(}y_{tt}-y^{-1}y_{t}^{2}\big{)}\mathrm{exp}(at+b)=0,$$

$$u_0=y,u_1=y_t,$$

$$\begin{array}{rl}{u_1}_t-{u_0}^{-1}{u_1}^2& =0,\\ {u_0}_t-u_1& =0,\\ u_0(0)=y_0,u_1& (0)=y_1.\end{array}$$

$$\hat{t}$$

$$\widehat{u_1}$$

$${\mathbb{V}}_q([0,T])=\{Y:{Y|}_{I_n}\in {\mathbb{P}}_q(I_n),n=0,...,N-1\}.$$

$${\mathbb{V}}_q^C([0,T])={\mathbb{V}}_q([0,T])\cap {\mathcal{C}}^0([0,T]),$$

$$\begin{array}{rl}{\int }_{I_n}\mathcal{F}(t,\mathbf{U},{U_m}_t)\varphi \mathrm{d}t& =0\forall \varphi \in {\mathbb{V}}_q(I_n),\\ {\int }_{I_n}({U_0}_t-U_1)\psi \mathrm{d}t& =0\forall \psi \in {\mathbb{V}}_q(I_n),\\ ⋮& \\ {\int }_{I_n}({U_{m-1}}_t-U_m)\chi \mathrm{d}t& =0\forall \chi \in {\mathbb{V}}_q(I_n),\end{array}$$

$$\begin{array}{rl}{\int }_{I_n}({U_1}_t-{U_0}^{-1}{U_1}^2)\varphi \mathrm{d}t& =0\forall \varphi \in {\mathbb{V}}_q(I_n),\\ {\int }_{I_n}({U_0}_t-U_1)\psi \mathrm{d}t& =0\forall \psi \in {\mathbb{V}}_q(I_n).\end{array}$$

$$\hat{t}$$

$$\widehat{U_1}$$

$$\begin{array}{rl}{\int }_{I_n}\exp (at+b)({U_1}_t-{U_0}^{-1}{U_1}^2+a{U_0}_t-a{U}_1)\varphi \mathrm{d}t& =0,\\ {\int }_{I_n}\exp (at+b)({U_0}_t-U_1)\psi \mathrm{d}t& =0,\end{array}$$

$$\begin{array}{rl}{\int }_{I_n}\exp (at+b)({U_1}_t-{U_0}^{-1}{U_1}^2)\varphi \mathrm{d}t& =0,\\ {\int }_{I_n}\exp (at+b)({U_0}_t-U_1)\psi \mathrm{d}t& =0.\end{array}$$

$$\hat{\mathbf{u}}=g\cdot \mathbf{u}$$

$${\mathbf{U}}^{(1)}=(U_0,\dots ,U_m,{U_0}_t,\dots ,{U_m}_t)$$

$$\widehat{{U_i}_t}=\frac{d{\hat{U}}_i}{d\hat{t}},i=0,\dots ,m.$$

$$\mathbf{\lambda }(H({\mathbf{Z}}^{(1)}))={\mathbf{\tau }}^{\ast }H({\mathbf{Z}}^{(1)})=H(\hat{\mathbf{Z}}{}^{(1)})$$

Full text

Department of Mathematics and Statistics, Memorial University of Newfoundland, St. John’s, NL, A1C 5S7, Canada abihlo@mun.ca.

Department of Mathematics and Statistics, Memorial University of Newfoundland, St. John’s, NL, A1C 5S7, Canada jjackaman@mun.ca.

Department of Mathematics, Monmouth University, West Long Branch, NJ, 07764, USA fvalique@monmouth.edu.

This research was supported, in part, thanks to the Canada Research Chairs, the InnovateNL LeverageR&D and NSERC Discovery grant programs

1. Introduction

For the accurate long time simulation of ordinary differential equations (ODEs), it is paramount that intrinsic properties of the underlying problem be preserved, giving rise to the field of geometric numerical integration, [24]. Well-known geometric numerical methods include symplectic integrators such as collocation and Runge-Kutta methods, [9, 24, 33, 49], Lie–Poisson structure preserving schemes, [52], energy-preserving methods, [46], general conservative methods, [50, 51], and more. By preserving certain characteristics of the differential equations, these geometric numerical schemes typically provide better global and long term results than their non-geometric analogues.

Over the last 30 years there has been a considerable amount of work focusing on the elaboration of finite difference numerical schemes that preserve the Lie point symmetries of differential equations, [1, 14, 15, 16, 17, 30, 31]. Using either infinitesimal generators or the method of equivariant moving frames, the algorithms for constructing symmetry-preserving finite difference schemes are now well-established. For a review and comparison of the two approaches, we refer the reader to [7]. For ordinary differential equations, symmetry-preserving numerical schemes have shown to give good numerical results, especially when solutions exhibit sharp variations or admit singularities, [11, 12, 31]. For partial differential equations, the numerical improvements are not as clear, [3, 30, 34, 48], and there remains a lot more to be done. Most of the work thus far has focused on evolutionary partial differential equations and it is now accepted that in order to preserve symmetries, numerical schemes have to be defined on time-evolving meshes. To avoid mesh tangling and other numerical instabilities, the basic invariant numerical schemes have to be combined with evolution–projection techniques, invariant r-adaptive methods, or invariant meshless discretisations, [2, 4, 5, 6].

Recently, the first and third authors began to adapt some of the techniques used to construct symmetry-preserving finite difference schemes to finite element methods, [8], focusing on second order ODEs and Burgers’ equation. From a numerical perspective, finite element methods offer several advantages over finite difference methods. Firstly, when considering complex domains, unstructured grids or moving boundaries, finite element methods are generally easier to implement than finite difference methods. Secondly, since finite element methods are based on the discretisation of the weak form of a differential equation, such numerical schemes have less rigid smoothness requirements than finite difference schemes. From an application standpoint, finite element methods are used in a wide range of fields such as fluid dynamics, engineering, physics, and applied mathematics.

Given a group of Lie point transformations acting freely and regularly on a manifold, the method of equivariant moving frames, originally developed in [21], and extended to Lie pseudo-groups, discrete groups, and finite differences in [37, 41, 42, 43], is a powerful tool for constructing invariant quantities of the group action. Indeed, the method of equivariant moving frames is equipped with a systematic process for constructing invariant quantities using the so-called invariantisation map, which sends non-invariant quantities to their invariant counterparts. Therefore, given a finite element method, which will in general not preserve the symmetries of the original differential equation, a symmetry-preserving scheme can be easily constructed by invariantising the non-invariant discrete weak formulation. This is the same key idea used to construct symmetry-preserving finite difference schemes in [7, 48, 30, 31],

In [8], the authors limited themselves to second order ODEs and only considered projectable group actions. In the present paper we generalise the ideas set out in [8] to ODEs of arbitrary orders and to general Lie point transformation groups. In particular, we utilise the conforming nature of the finite element approximation to exploit the invariantisation procedure within a continuous framework, whereas in the prequel a finite difference like approach was considered. We note that viewing high degree finite elements in the difference setting leads to complex calculations, which adds computational challenges to the methodology considered in [8]. The approach introduced in this paper has several benefits, not least of which is the ability to construct invariant finite element approximations of arbitrarily high order simultaneously. Though, we note that for some group actions it will be necessary to consider certain aspects of the invariantisation procedure for each polynomial degree independently. In these cases, these degree specific modifications are significantly simpler than those introduced in [8].

While finite element methods were originally introduced as a methodology for spatial discretisation, they are also ideal for temporal discretisations. In particular, by rewriting an arbitrarily high order ODE as a system of first order equations we may introduce both continuous and discontinuous time-stepping Galerkin approximations which are well studied, see [28, 22, 19, 18]. It is known that for a large class of geometric ODEs, the continuous Galerkin method preserves the underlying energy of the problem, [25], while this is not the case for the discontinuous Galerkin method, [20]. Therefore, in the sequel we will disregard the discontinuous method and focus on its conservative continuous counterpart. Furthermore, the continuous Galerkin approximation is known to be “close” to a family of symplectic collocation methods, see [27, Chapter 3.1].

The remainder of this work is set out as follows: In Section 2, after introducing notation and giving necessary definitions, we outline a general methodology for constructing a discrete finite element method that preserves the Lie point symmetries of an ODE of arbitrary order. In Section 3 we examine several examples highlighting the intricacies of the methodology. In Section 4 we provide numerical experiments that highlight favourable aspects of symmetry-preserving finite element schemes. In particular, our simulations show that for certain equations the invariant schemes are more accurate long term and can be implemented on larger elements than their non-invariant counterparts. Finally, in Section 5 we contrast the ideas introduced in this paper with those first considered in [8]. In particular we notice that our methodology extends the work of [8] to non-projectable symmetry actions and is much easier to generalise to higher order ODEs.

Before diving into the subject, we note that the present paper should be regarded as an exploratory study in assessing the feasibility of systematically incorporating symmetries into the finite element method. Since this paper shows that this is possible, the next natural step will be to consider applications to partial differential equations, which offers a much richer playing ground.

2. General methodology

Let

[TABLE]

be an $(m{+}1)$-th order ODE with initial conditions

[TABLE]

To simplify the notation, we introduce the $(m{+}1)$-th jet notation

[TABLE]

which collects the dependent variable $y$ and its time derivatives up to order $m+1$. The tuple $(t,y^{(m+1)})$ provides local coordinates for the $(m{+}1)$-th order extended jet bundle $\mathrm{J}^{(m+1)}(\mathbb{R}^{2},1)$ of curves in the plane, [40].

We assume that equation (1) admits a certain Lie group of point symmetries, which are transformations depending on the independent and dependent variables that map solutions of (1) to themselves. There are many excellent textbooks on this topic, such as [10, 40, 26].

Definition 2.1**.**

A Lie group $G$ is said to be a symmetry group of equation (1) if and only if for all $g\in G$ near the identity, the equality

[TABLE]

holds.

The symmetry group (or at least the infinitesimal symmetry generators) of a given differential equation can easily be computed using symbolic software packages such as Maple, Mathematica, or Sage. Given the action of $G$ on the independent and dependent variables

[TABLE]

the prolonged action on the time derivatives is obtained by implicit differentiation

[TABLE]

In the sequel, we rewrite the $(m{+}1)$-th order ODE (1) as a system of $m{+}1$ first order differential equations by introducing the variables

[TABLE]

to obtain the system

[TABLE]

subject to the initial data

[TABLE]

In light of how $\bm{u}$ is defined in (6), the induced action of $G$ on $\bm{u}$ is given by the prolonged action (5) after a change of variables.

Remark 2.2**.**

In general, the system of equations (7) may possess more symmetries than the original equation (1). In the following, when working with system (7), we restrict ourselves to the symmetry group of the original equation (1).

Example 2.3** **(A working example).

To illustrate the notions introduced throughout this section, we will consider the nonlinear differential equation

[TABLE]

subject to the initial conditions

[TABLE]

for some prescribed constants $y_{0},y_{1}$, and with the assumption that $y(t)\neq 0$ for all $t\in[0,T]$, for some $T>0$. The differential equation (8) admits a six-parameter symmetry group of projectable transformations given by

[TABLE]

where $\alpha\delta-\beta\gamma=1$, $a,b\in\mathbb{R}$, and $\lambda>0$. In the following we will restrict ourselves to the two-parameter subgroup

[TABLE]

This makes the exposition easier to follow by keeping the formulas simple, and furthermore we will see in Section 4 that the finite element method that preserves this two-parameter symmetry group yields desirable numerical results. Up to order two, the prolonged action of (10) is simply obtained by computing the product rule for differentiation:

[TABLE]

Then

[TABLE]

provided that $y$ is a solution of (8), which confirms that (10) is a symmetry (sub-)group of the equation.

Equation (8) is written as a system of first order ODE by introducing the variables

[TABLE]

so that (8) becomes

[TABLE]

Making the change of variables (12) into the transformation formulas (10), we conclude that the system of equations (13) is invariant under the transformation group

[TABLE]

2.1. Temporal finite element methods

Let $T>0$, and consider a partition of the temporal interval $[0,T]$ into $N$ sub-intervals given by $0=t_{0}<t_{1}<...<t_{N}=T$. For $n\in\{0,\ldots,N-1\}$, let $I_{n}=(t_{n},t_{n+1}]$ denote the finite element of length $\tau_{n}=t_{n+1}-t_{n}$. Given a continuous or differentiable function $y(t)$, we systematically use capital letters to denote the functions finite element approximation. For example, $Y=Y(t)$ represents the finite element function approximating $y(t)$. When there is no ambiguity, we shall not explicitly write the time dependency of functions.

Definition 2.4**.**

Let $\mathbb{P}_{q}(I_{n})$ denote the space of polynomials of degree $q$ on the element $I_{n}\subset[0,T]$. Then the discontinuous finite element space is

[TABLE]

The continuous finite element space is defined analogously with global continuity enforced, i.e.

[TABLE]

where $\mathcal{C}^{0}([0,T])$ is the space of continuous functions over the interval $[0,T]$.

We let $\mathbb{V}_{q}^{C}(I_{n})$ represent the localisation of the continuous finite element space to a single element. Here, the values of functions at time $t_{n}$ in $\mathbb{V}_{q}^{C}(I_{n})$ are fixed by the right endpoint of the corresponding functions in $\mathbb{V}_{q}^{C}(I_{n-1})$, for $\mathbb{V}_{q}^{C}(I_{0})$ the initial function values are enforced by the initial conditions (2). We note that $\mathbb{V}_{q}(I_{n})\subset\mathcal{C}^{q}(I_{n})$ and $\mathbb{V}_{q}^{C}(I_{n})\subset\mathcal{C}^{q}(I_{n})$, i.e., both discontinuous and continuous finite element approximations are smooth over a single element. This property will be crucial in the sequel.

As we are restricting ourselves to initial value problems, we wish for our finite element approximations to boast a time stepping implementation. Moreover, as we have introduced the auxiliary variables (6), we rewrite the ODE (1) as the first order system (7) with dependent variable $\bm{u}$. This allows for the design of finite element approximations which encompass all initial value problems for ODEs. As mentioned in the introduction, we consider the continuous Galerkin method, which is well studied for non-degenerate temporal ODEs, [19]. For such problems the method is described by seeking for $\bm{U}\in(\mathbb{V}_{q+1}^{C}(I_{n}))^{m+1}$, where $n=\{0,\ldots,N-1\}$, such that

[TABLE]

where $\bm{U}(0)$ is fixed by the initial condition $\bm{u}(0)$, and $\bm{U}(t_{n})$ is determined by the solution at the right endpoint on the previous element $I_{n-1}$. For a detailed discussion and analysis of this method see [19, 22]. The functions $\bm{U}$ are typically referred to as trial functions, and $\phi,\psi,..,\chi$ are known as test functions. We note that the upwind discontinuous Galerkin approximation, [18], is another widely used finite element method for time evolving ODEs. However, this method does not conserve any underlying structures associated to the continuous problem, [20]. In this work, we will refer to the continuous Galerkin formulation as “standard” in view of the following remark.

Remark 2.5** **(**Preservation of geometric structures by

(17)).**

The standard continuous Galerkin formulation (17) is well studied and understood. In fact, for a large class of geometrically interesting ODEs, namely Hamiltonian ODEs, it exactly preserves the associated energy at the points $t_{n}$, see [25, 27]. For an introduction to Hamiltonian ODEs and their geometric properties see [24, 33]. In the lowest order case, i.e., when $q=0$, the standard finite element method is equivalent to the average vector field (AVF) discrete gradient method [38], and may therefore be viewed as a generalisation of this method. Furthermore, we note the close relation between the standard method and well known symplectic integrators. Through an appropriate choice of Gaussian quadrature, i.e., choosing the number of quadrature points equal to the polynomial degree of the numerical solution, we obtain the $q+1$ point Gauss collocation method. This collocation method is, in turn, equivalent to a member of the Gauss family of Runge-Kutta methods which are known to be symplectic. For more information see [27, §3.1]. With this in mind, we expect the “standard” formulation to perform well numerically.

We note that throughout this work we shall not study this method under quadrature approximation, instead evaluating integrals exactly where possible, or employing a quadrature which has a negligible error. This choice will be fully justified in Remark 4.2.

As the next example shows, the approximate weak formulation (17) will, in general, not preserve the symmetries of the original system of differential equations (6).

Example 2.6**.**

Continuing Example 2.3, the “standard” finite element formulation of the first order system (13) is obtained by seeking for ${U_{0}},{U_{1}}\in\mathbb{V}_{q+1}^{C}(I_{n})$ such that

[TABLE]

Observe that if ${U_{0}}(t)=0$ for some $t\in I_{n}$, then this finite element is not well defined, so we must assume that this is not the case numerically. However, we do not have an analytic assurance of this requirement. Due to the local smoothness of the approximation, the group action on $(t,U_{0},U_{1},{U_{1}}_{t})$ is obtain by simply replacing the lower case $u$ in (14) by its capital case counterpart $U$. Thus,

[TABLE]

Therefore, the transformed weak form is

[TABLE]

observing that $\phi$ and $\psi$ are both arbitrary functions in the same function space. The transformed finite element method reduces to

[TABLE]

As the exponential function does not commute with the integral, we conclude that the standard finite element method is not invariant under the group action (19).

To construct finite element methods that preserve Lie point symmetries, we will use the method of equivariant moving frames, which we now summarise in the particular context of our problem.

2.2. Moving frames and invariantisation

For an introduction to the method of equivariant moving frames, we refer the reader to the original papers [21, 37, 43, 42] and the textbook [36].

Let $G$ be the symmetry group of the differential equation (1), or possibly a symmetry subgroup. Via the prolonged action, it induces an action

[TABLE]

on the variables (6), and remains a symmetry group of (7). Since the function $\bm{U}=(U_{0},\ldots,U_{m})$ is smooth in the interior of any element $I_{n}$, the action (22) also applies to $\bm{U}$ inside $I_{n}$, to give $\widehat{\bm{U}}=g\cdot\bm{U}$. It is prolonged to

[TABLE]

using the usual prolongation formula

[TABLE]

Therefore, in the following, we consider the action of the Lie group $G$ on the first order jet space \mathrm{J}^{(1)}_{I_{n}}=\mathrm{J}^{(1)}(\mathbb{R}^{m+1},1)\big{|}_{I_{n}} with local coordinates $\bm{Z}^{(1)}=(t,\bm{U}^{(1)})$. We combine the Lie group and the first order jet space by introducing the first order lifted bundle $\mathcal{B}^{(1)}_{I_{n}}=G\times\mathrm{J}^{(1)}_{I_{n}}$. This bundle admits a Lie groupoid structure, [35], with source map $\bm{\sigma}(g,\bm{Z}^{(1)})=\bm{Z}^{(1)}$ corresponding to the projection onto $\mathrm{J}^{(1)}_{I_{n}}$ and the target map $\bm{\tau}(g,\bm{Z}^{(1)})=\widehat{\bm{Z}}{}^{(1)}=g\cdot\bm{Z}^{(1)}$ given by the prolonged group action.

Definition 2.7**.**

Let $H\colon\mathrm{J}^{(1)}_{I_{n}}\to\mathbb{R}$ be a function. The lift of $H$ is the function

[TABLE]

obtained by substituting the arguments $\bm{Z}^{(1)}$ of $H$ by their prolonged action expressions $\widehat{\bm{Z}}{}^{(1)}$.

The lift of a function is a new function defined on the lifted bundle $\mathcal{B}^{(1)}_{I_{n}}$. The lift map $\bm{\lambda}$ extends to differential forms, see [32] for more detail. The (horizontal) lift of $\thinspace\mathrm{d}{t}$ is the one-form

[TABLE]

obtained by taking the (horizontal) differential of $\widehat{t}=g\cdot t$, where the group parameters $g=(g^{1},\ldots,g^{r})$ are treated as constants. In order to have a well defined action of the symmetry group $G$ on the system of equations (17), there remains to define the action of $G$ on the basis functions $\phi\in\mathbb{V}_{q}(I_{n})$. By definition, the function $\phi$ is a degree $q$ Lagrangian interpolating function depending on the nodes of the element $I_{n}$, i.e., $t_{n}$, $t_{n+1}$, a finite number of interior points $t_{i}\in(t_{n},t_{n+1})$, and obviously the continuous time variable $t$. We write these dependencies explicitly as

[TABLE]

where in $t_{n},\ldots,t_{n+1}$, we collect the dependency of $\phi$ on the nodes $t_{n}$, $t_{n+1}$ and the interior points $t_{i}$. Then, the lift of $\phi$ is defined as the function

[TABLE]

where the group $G$ acts via the product action on $t_{n}$, $\ldots$, $t_{n+1}$, and $t$, simultaneously.

To simplify the exposition, we now introduce the notation

[TABLE]

to refer to the system of equations (17), where $\boldsymbol{\phi}=(\phi,\psi,\ldots,\chi)$ denotes the test functions occurring in (17). Then the lift of (17), in other words, the action of $G$ on (17) is defined as

[TABLE]

where $\omega$ is given in (23).

Definition 2.8**.**

The system of equations (17) is said to be invariant under the action of a Lie group $G$ if and only if

[TABLE]

As mentioned in the previous section, the approximate weak form (17) will not, in general, preserve the symmetries of the original system (6). To obtain a weak formulation that preserves the symmetries, we follow the methodology of invariantisation, which is based on the introduction of a moving frame.

Definition 2.9**.**

Let $G$ be a Lie group acting on $\mathrm{J}^{(1)}_{I_{n}}$. A right moving frame is a $G$-equivariant map $\rho\colon\mathrm{J}^{(1)}_{I_{n}}\to G$ satisfying the equality

[TABLE]

To guarantee the existence of a moving frame, we need to impose two regularity assumptions on the group action, [21].

Definition 2.10**.**

The Lie group $G$ is said to act freely at $\bm{Z}^{(1)}\in\mathrm{J}^{(1)}_{I_{n}}$ if the isotropy group

[TABLE]

is trivial. The group action is said to be locally free at $\bm{Z}^{(1)}$ if $G_{\bm{Z}^{(1)}}$ is a discrete subgroup. We say that $G$ acts (locally) freely on $\mathrm{J}^{(1)}_{I_{n}}$ if the action of $G$ is (locally) free at every point $\bm{Z}^{(1)}\in\mathrm{J}^{(1)}_{I_{n}}$.

Definition 2.11**.**

The action of $G$ on $\mathrm{J}^{(1)}_{I_{n}}$ is said to be regular if the orbits form a regular foliation, that is to say that the orbits are submanifolds of dimension $r=\dim G$.

Theorem 2.12**.**

Let $G$ act (locally) freely and regularly on $\mathrm{J}^{(1)}_{I_{n}}$, then a moving frame exists in a neighbourhood of every point $\bm{Z}^{(1)}\in\mathrm{J}^{(1)}_{I_{n}}$.

The construction of a moving frame is based on the introduction of a cross-section.

Definition 2.13**.**

A cross-section $\mathcal{K}\subset\mathrm{J}^{(1)}_{I_{n}}$ is a submanifold of co-dimension $r=\dim G$ transverse to the group orbits.

Given a cross-section $\mathcal{K}\subset\mathrm{J}^{(1)}_{I_{n}}$, the right moving frame at the point \bm{Z}^{(1)}=\big{(}t,\bm{U}^{(1)}\big{)} is defined as unique group element $\rho\in G$ such that

[TABLE]

In practice, a moving frame is frequently obtained by choosing a coordinate cross-section, which consists of setting certain coordinates of the jet $\bm{Z}^{(1)}=(t,\bm{U}^{(1)})$ equal to constant values. Due to the definition of the variables $\bm{u}$ in (6), and the fact that the transformation group $G$ comes from the symmetry group of the differential equation (1), we work under the assumption that it is possible to determine a coordinate cross-section by setting $r=\dim G$ coordinates from $\bm{Z}^{(0)}=(Z_{0},\ldots,Z_{m+1})=(t,\bm{U}^{(0)})$ to constant values. Thus, let

[TABLE]

where $0\leq i_{1}<i_{2}<\cdots<i_{r}\leq m+1$. The requirement (29) leads to the normalisation equations

[TABLE]

Solving for the group parameters $g$ yields the moving frame $g=\rho(\bm{Z}^{(1)})=\rho(\bm{Z}^{(0)})$.

Example 2.14**.**

Continuing Example 2.6, we now construct a moving frame for the action (19). Working on the set where $U_{0}\neq 0$, the action (19) is free and regular, therefore a moving frame exists. A coordinate cross-section to the group orbits is given by

[TABLE]

The corresponding normalisation equations are

[TABLE]

Solving these two equations for the group parameters $a,b$ leads to the right moving frame

[TABLE]

Given a moving frame, there is a systematic procedure for mapping non-invariant quantities (e.g. functions, differential forms, functionals) to their invariant analogue, known as invariantisation. To introduce the invariantisation map, we note that a moving frame $\rho\colon\mathrm{J}^{(1)}_{I_{n}}\to G$ induces a moving frame section $\varrho\colon\mathrm{J}^{(1)}_{I_{n}}\to\mathcal{B}^{(1)}_{I_{n}}$ in the lifted bundle $\mathcal{B}^{(1)}_{I_{n}}$ defined by

[TABLE]

Then, given a function $H\colon\mathcal{B}^{(1)}_{I_{n}}\to\mathbb{R}$, or more generally a differential form, we can consider the moving frame pull-back

[TABLE]

obtained by replacing the group element $g$ by the moving frame $\rho(\bm{Z}^{(1)})$.

Definition 2.15**.**

Let $\rho\colon\mathrm{J}^{(1)}_{I_{n}}\to G$ be a moving frame. The invariantisation of the system of equations (17) is defined as

[TABLE]

where we define $\varrho^{*}[\bm{\lambda}(\boldsymbol{\phi})]=\boldsymbol{\phi}(\rho\cdot t_{n},\ldots,\rho\cdot t_{n+1};\rho\cdot t)$.

Therefore, the system (35) is obtained by first lifting its arguments, i.e. by acting on the arguments with the symmetry group $G$, and then replacing the group element $g$ by the moving frame $\rho$. The map $\iota$ is called the invariantisation map. The fact that (35) is invariant follows from the $G$-equivariance of the moving frame $\rho$. We note that if the finite element approximation is already invariant, then the invariantisation of (35) will simply yield back the original system of equations.

Example 2.16**.**

As observed in Example 2.6, the discrete weak formulation (18) is not invariant under the group of transformations (19). To obtain a weak formulation that will remain invariant under the group action, it suffices to invariantise (18). The lift of the discrete weak form has already been computed in (21). Using the moving frame (34) computed in Example 2.14, the moving frame pull-back of (21) is simply obtained by substituting the group normalisations (34) into (21). The result is the symmetry-preserving finite element scheme

[TABLE]

The invariance of the resulting scheme may be verified simply by applying the group action (19).

To summarise the previous two sections, we now carefully discuss the invariantisation procedure of a finite element approximation and its various intricacies. More examples of the invariantisation procedure are considered in Section 3.

2.3. Outline of the methodology for constructing symmetry-preserving finite element methods

To simplify the presentation, we restrict our exposition to third order ODEs, which encapsulates all examples considered in Section 3. Following the exposition introduced in the previous sections, it is straightforward to extend the methodology to higher order equations.

Our starting point is a third order ODE

[TABLE]

with initial condition $y(0)=y_{0}$, $y_{t}(0)=y_{1}$, $y_{tt}(0)=y_{2}$, which admits a symmetry group $G$. The steps for constructing a symmetry-preserving continuous Galerkin finite element scheme are as follows:

Step 1

Introduce the auxiliary variables

[TABLE]

and recast the third order equation (37) as a system of first order ODEs

[TABLE]

with initial conditions $u_{0}(0)=y_{0}$, $u_{1}(0)=y_{1}$, $u_{2}(0)=y_{2}$. We note that while $\mathcal{F}(t,u_{0},u_{1},u_{2},{u_{2}}_{t})$ is uniquely prescribed here, there may be other, equally valid, ways to describe the system of equations (39). Such an example is provided in Example 3.3.

From the definition of the variable $u$ in (38), the symmetry group $G$ induces the transformations

[TABLE]

via the prolonged action. Furthermore,

[TABLE]

Step 2

Formulate a “standard” finite element approximation in the spirit of (17) by seeking $U_{0},U_{1},U_{2}\in\mathbb{V}_{q+1}^{C}(I_{n})$ such that

[TABLE]

where $U_{0}(t_{n}),U_{1}(t_{n}),U_{2}(t_{n})$ are fixed by either the solution on the previous element or the initial data.

Step 2.1

This step is not necessary, but if possible we wish for our space of test functions $\mathbb{V}_{q}(I_{n})$ to be invariant under the group action $\widehat{t}=g\cdot t$. If this is not the case, it may be possible to choose a more appropriate function space for a given symmetry group action, however, no such examples are presented in the sequel.

Step 3

Construct a moving frame. Assuming the action of $G$ on

[TABLE]

is free and regular, choose a cross-section

[TABLE]

and solve the normalisation equations $g\cdot Z_{i_{1}}=c_{i_{1}},\ldots,g\cdot Z_{i_{r}}=c_{i_{r}}$ for the group parameters to obtain the moving frame $\rho(\bm{Z}^{(0)})$.

Step 4

Invariantise (40) according to the invariantisation formula (35), which is obtained by first lifting $(t,\bm{U}^{(1)})$, $\boldsymbol{\phi}=(\phi,\psi,\chi)$, and $\mathrm{d}t$ using the group action $G$ and then substituting the group elements my their moving frame expressions. The result is a discrete weak formulation of (39) that is invariant under the action of the symmetry group of the differential equation (37).

Remark 2.17** **(Normalisation equations without exact solutions).

Steps 3 and 4 are based on the assumption that the normalisation equations can be solved for the group parameters to obtain an explicit moving frame $\rho(\bm{Z}^{(0)})$, which is then used to invariantised (40). If the normalisation equations cannot be solved analytically, it is still possible to obtain a symmetry-preserving finite element method as follows. Instead of steps 3 and 4, consider the enlarged system of equations consisting of the normalisations with the lift of (17):

[TABLE]

Solving these equations numerically to a high precision will then constitute the desired invariant finite element method.

Remark 2.18** **(Generalisation to systems of ODEs).

In the above exposition, we have restricted our attention to a single ODE. The procedure for constructing a symmetry-preserving finite element method can easily be extended to systems of ODEs. From a notational standpoint, simply replace the original ODE (1) by a system of $\ell$ ODEs and let $y(t)=(y^{1}(t),\ldots,y^{\ell}(t))$ denote a vector of dependent variables. Then the rest of the exposition remains essentially unchanged, except that $u_{i}=y_{t^{i}}$ now denotes the time derivative of an $\ell$-dimensional vector and the discontinuous and continuous finite element spaces should be replaced by the Cartesian products

[TABLE]

3. Examples

In this section we provide several examples highlighting different aspects of the methodology for constructing symmetry-preserving continuous Galerkin finite element methods. Our first example is concerned with second order linear equations which occur in classical mechanics, electric circuits, and many other branches of science. Our second example involves the Schwarzian differential equation which occurs in geometry, the theory of Sturm–Liouville equations, [45], and in the study of gravity–dilation–antisymmetric tensor systems, [53, 54]. Finally, our last two examples have been selected so as to admit interesting symmetry groups and to illustrate certain aspects of our methodology. In Section 4, numerical simulations using the obtained schemes are performed and compared to standard non-invariant finite element methods.

Example 3.1** **(Second order linear ODE).

Our first example is not an illustration of the invariantisation procedure per say, but rather an illustration of the fact that for certain differential equations the continuous Galerkin method is naturally invariant. To this end, let $y=y(t)$ satisfy the initial value problem

[TABLE]

for some prescribed constants $y_{0}$ and $y_{1}$. This equation possesses a two-parameter symmetry group given by

[TABLE]

where $\epsilon_{1},\epsilon_{2}\in\mathbb{R}$, and $\alpha(t)$, $\gamma(t)$ are linearly independent solutions of the corresponding homogeneous equation

[TABLE]

This symmetry group reflects the fact that the differential equation is linear. Introducing the auxiliary variables ${u_{0}}=y$ and ${u_{1}}=y_{t}$, we may rewrite (43) as the system of first order equations

[TABLE]

which is invariant under the two-parameter group of transformations

[TABLE]

Applying the finite element discretisation (17) results in us seeking for ${U_{0}},{U_{1}}\in\mathbb{V}_{q+1}^{C}(I_{n})$ such that

[TABLE]

Since the independent variable $t$ is an invariant of the group action, it follows that $\thinspace\mathrm{d}{t}$ and the test functions $\phi$ and $\psi$ remain unchanged under the group action (47). Thus, the transformed finite element scheme is simply obtained by substituting the transformation rules

[TABLE]

and

[TABLE]

into (47). The result is

[TABLE]

Since $\alpha$ and $\gamma$ are solutions of the homogeneous problem (45), it follows that the terms in $\epsilon_{1}$ and $\epsilon_{2}$ vanish, and therefore, the finite element scheme (48) is invariant under the group action (47) without any need to invariantise the functional. In other words, since the scheme is linear it preserves the superposition principle.

Example 3.2** **(Schwarzian differential equation).

Consider the third order ODE

[TABLE]

subject to the initial data $y(0)=y_{0},y_{t}(0)=y_{1},y_{tt}(0)=y_{2}$, and where we assume that $y_{t}(t)\neq 0$ for all $t\in[0,T]$. The differential equation is known, [36], to be invariant under the linear fractional group action

[TABLE]

By letting ${u_{0}}=y$, we can rewrite (50) as the system of first order differential equations

[TABLE]

which is invariant under the extended group of transformations

[TABLE]

The “standard” finite element approximation is obtained by seeking for ${U_{0}},{U_{1}},{U_{2}}\in\mathbb{V}_{q+1}^{C}(I_{n})$ such that

[TABLE]

subject to appropriate initial data. Applying the group action (53) to (54), we find the transformed functionals

[TABLE]

which shows that the weak formulation is not invariant. We observe that this finite element approximation is consistent, i.e., substituting sufficiently globally smooth trial and test functions into the numerical approximation gives back the original system of ODEs (52).

A moving frame is obtained by choosing the cross-section

[TABLE]

This yields the normalisation equations

[TABLE]

the solution of which gives the moving frame

[TABLE]

The invariantisation of the weak form (54) is then obtained by substituting the group normalisations (58) into the transformed functionals (55). The result is the symmetry-preserving finite element method

[TABLE]

Example 3.3** **(A second order quasi-linear ODE).

We now move our attention to the second order ODE

[TABLE]

for some constants $y_{0},y_{1}$. This differential equation admits a two-parameter symmetry group with nontrivial temporal action given by

[TABLE]

where $a,b\in\mathbb{R}$. When considering a finite domain in both the continuous and discrete setting, it is important to note that due to dilation and translation in time, the action will accordingly dilate and shift the domain of consideration. With this in mind, the “invariant” scheme we will consider for this problem will be permitted to dilate and shift the elements through time. Introducing the auxiliary variable ${u_{0}}=y$ we may write (60) as

[TABLE]

which is invariant under the extended group action

[TABLE]

We note that the first equation in (62) is not of the form discussed in Section 2.3, as it contains the term ${u_{0}}_{t}$. Here we include terms in ${u_{0}}_{t}$ to simplify the computations. The first order system

[TABLE]

would be equally valid. However, since the group action on (64) leads to significantly more complex expressions, we do not consider this system here. In the spirit of (17), we introduce a standard finite element approximation by seeking for ${U_{0}},{U_{1}}\in\mathbb{V}_{q+1}^{C}(I_{n})$ such that

[TABLE]

For any function $f(t)\in\mathbb{V}_{q}(I_{n})$ we observe that $f(g\cdot t)$ remains in $\mathbb{V}_{q}$ up to translation in the time variable $t$. Therefore, our space of test functions is preserved up to temporal translations and we conclude that our space of test functions is “invariant” under this group action as discussed in Step 2.1 of Section 2.3. Acting on the finite element approximation (65) with (63) we find, after simplification, the transformed functionals

[TABLE]

We observe consistency of the first equation by equating ${U_{0}}_{t}$ and ${U_{1}}$. However, the ${U_{1}}$ term scales differently in $t$ from the second equation and so cannot be substituted.

To construct a moving frame, we choose the cross-section

[TABLE]

The resulting normalisation equations are

[TABLE]

Solving for the group parameters $a,b$, under the assumptions that ${U_{0}}\neq 0$ and ${U_{0}}\neq-\cfrac{t{U_{1}}}{2}$, we obtain the moving frame

[TABLE]

Substituting these normalised group parameters into (66) yields the invariant finite element approximation

[TABLE]

where ${U_{0}},{U_{1}}\in\mathbb{V}_{q+1}^{C}(I_{n})$.

Example 3.4** **(A non-projectable action).

We now draw our attention to the first order ODE

[TABLE]

for a known fixed constant $C$, subject to the initial $y(0)=y_{0}$. This ODE is invariant under the two-parameter non-projectable symmetry group action

[TABLE]

where $\alpha,\beta\in\mathbb{R}$. As the ODE is a first order system we do not need to introduce auxiliary variables to propose a finite element discretisation. However, for consistency with the remainder of this work we shall write $y(t):={u_{0}}(t)$. A standard finite element approximation for (71) is given by seeking for ${U_{0}}\in\mathbb{V}_{q+1}^{C}(I_{n})$ such that

[TABLE]

Applying the prolonged group action

[TABLE]

to the finite element approximation (73) yields

[TABLE]

The induced action of the transformation $\widehat{t}=t+\alpha{U_{0}}$ on the test functions $\phi$ will depend on the degree of the functions considered. First, let us consider the case where the test functions are naturally invariant, which occurs when $q=0$ and the test functions are time independent. In this case $\widehat{\phi}=\phi$ and we obtain the transformed finite element scheme

[TABLE]

where ${U_{0}}\in\mathbb{V}_{1}^{C}(I_{n})$. To construct a moving frame, we choose the cross-section

[TABLE]

The corresponding normalisation equations are

[TABLE]

Solving for the group parameters yields the moving frame

[TABLE]

Substituting (79) into (76) gives the invariant scheme where one seeks for ${U_{0}}\in\mathbb{V}_{1}^{C}(I_{n})$ such that

[TABLE]

For $q>0$, the invariantisation of the finite element scheme is more complicated, as the space of test functions is not invariant. In fact, as the group action is not projectable there is no known conventional choice of function space that is invariant under the symmetry group. We must instead invariantise every basis function that spans our function space, which is encapsulated into the general invariantisation formula (35). We note that after the invariantisation of the basis functions we no longer expect them to span a function space. Nevertheless, implementing the invariantisation procedure will produce a consistent symmetry-preserving finite element method. As an illustrative example, consider the case where $q=1$. Here the space of test functions $\mathbb{V}_{1}(I_{n})$ may be given in terms of the Lagrange basis functions

[TABLE]

where

[TABLE]

Letting

[TABLE]

with $a,b\in\mathbb{R}$, the transformed finite element method (75) becomes

[TABLE]

where

[TABLE]

and ${U_{0}}\in\mathbb{V}_{2}^{C}(I_{n})$. Substituting the moving frame expressions (79) into (83) yields the following invariant finite element approximation: Seek ${U_{0}}\in\mathbb{V}_{2}^{C}(I_{n})$ such that for all $a,b\in\mathbb{R}$,

[TABLE]

where

[TABLE]

We note that at the endpoints of the element, i.e. at $t_{n}$ and $t_{n+1}$, the value of modified basis functions $\mathcal{M}_{i}$ is the same as with the Lagrange basis functions. However, we do not expect them to form a partition of unity. This procedure may be replicated for arbitrarily high order Lagrange basis functions, leading to arbitrarily high order invariant finite element approximations.

Remark 3.5** **(**The connection between Lie point symmetry preserving

methods and other geometric numerical integrators).**

We observe that, in view of Remark 2.5, under a $q+1$ point quadrature approximation the standard and invariant finite element approximations are equivalent, following the methodology outlined in the proof of [27, Theorem 3.1.9]. This suggests that the perturbation we are making to our numerical method to preserve the Lie symmetry is small. In fact, within the context of Hamiltonian ODEs the standard finite element method is energy preserving, indicating that Lie point symmetry preserving integrators are “close” to energy preserving integrators. Furthermore, we note both finite element approximations under this quadrature yield a symplectic collocation method [27, §3.1], suggesting that Lie point symmetry preserving integrators are “close” to symplectic integrators. We conjecture that, within the field of geometric numerical integration, Lie point symmetry preserving schemes are competitive alternatives to energy preserving and symplectic methods. Significant theoretical work beyond the scope of this paper is required to verify this conjecture.

4. Numerical experiments

In this section, we consider select numerical experiments showcasing the impact of considering symmetry-preserving finite element methods. In particular, we assume that our element size $\tau_{n}=t_{n+1}-t_{n}=\tau$ is uniform. We note that there are no theoretical requirements for considering a uniform element size. This simplification is made solely to add clarity to the results of our numerical experiments.

We measure all errors in the $L_{2}$ norm, which for a vectorial solution is defined as

[TABLE]

where $\bm{U}=(U_{0},...,U_{m})$ is the numerical solution and $\bm{u}=(u_{0},...,u_{m})$ is the corresponding exact solution. We also consider the maximal error at the temporal nodes, which we shall use to investigate nodal super-convergence. It is important to remark that the maximal nodal error only induces a norm when $q=0$, however, we shall not use it to investigate convergence rates in the sequel instead relying on the $L_{2}$ norm (87).

We use the $L_{2}$ errors to compute an experimental order of convergence (EOC), which is defined as follows.

Definition 4.1** **(Experimental order of convergence).

Given two finite sequences $\{a_{k}\}_{k=0}^{n}$ and $\{b_{k}\}_{k=0}^{n}$, we define the experimental order of convergence (EOC) as

[TABLE]

In the sequel $a_{k}$ will represent errors and $b_{k}$ will represent element sizes.

In practice, our finite element approximation is solved through the following steps:

1. (1)

The problem is linearised, allowing us to solve an underlying linear problem which converges to the solution of the nonlinear problem. In our numerical experiments we employ a Newton solver, which is solved up to a tolerance of $10^{-12}$. 2. (2)

Our finite element functions are decomposed into basis functions, for example $U(t)=U^{1}\mathcal{L}_{1}(t)+\cdots+U^{q+2}\mathcal{L}_{q+2}(t)$, where $\mathcal{L}_{i}(t)$ are the degree $q{+}1$ Lagrange basis functions over the element $I_{n}$ and $U^{i}$ are the values of the finite element function at the Lagrange points, otherwise known as the degrees of freedom. Additionally, decomposing the test functions into their basis functions we may assemble a linear system of equations $A\bm{U}=b$, where $\bm{U}=(U^{1},...,U^{q+2})$, $A$ is a matrix that represents terms depending on both $U$ and the test functions, and $b$ is a vector involving terms depending solely on the test functions and time. We observe that for finite element formulations of the form (17), the linear system $A\bm{U}=b$ is square after the enforcement of initial data. 3. (3)

To solve the finite element approximation we must evaluate the integral, while it often may be computed exactly, in practice it is typically evaluated through a quadrature approximation. In the sequel we employ an order $16$ Gauss quadrature, the error of which we expect to be below machine precision when the quadrature is not exact. 4. (4)

Finally, we solve the linear system iteratively to obtain the finite element solution over one element.

For more information on the practical implementation of finite element methods see [13, §0]. In the sequel we shall use Firedrake, [47], to conduct steps (2) and (3), and utilise the NumPy linalg routine solve, [39], for (4).

Remark 4.2** **(Computational costs).

In step (3), we note that the quadrature approximation computed is of significantly higher order than the numerical method. While the quadrature approximation is of higher order, this does not significantly contribute to the computational complexity of the method. We also note that steps (1)–(3), up to the inclusion of initial data, may be conducted entirely independently of the particular time step under consideration. In fact, the linear system solved iteratively in step (4) is identical on each time step up to the enforcement of initial data. The linearisation of the method and cost of assembling basis functions and evaluating the integrals are an initial set-up cost of the method and may be conducted offline. In practice, we are only required to enforce different initial data on the linear system at each time step. Finally, we note that the Newton solver does have a leading order computation cost as it is employed at every time step, which is typical of nonlinear numerical approximations.

Example 4.3** **(Working example).

Here we compare standard finite element approximation (18) of the working example (13) with the symmetry-preserving finite element scheme (36), which is investigated as an illustrative example throughout Section 2. We initialise both finite element schemes using

[TABLE]

which approximates the exact solution

[TABLE]

i.e., we wish to numerically simulate exponential decay, with this in mind, we consider a relatively short time simulation with $T=10$, as over significantly longer time the exact solution will, up to numerical precision, be zero. Simulating the standard finite element scheme (18) for various polynomial degrees we obtain Table 2. Similarly simulating the invariant finite element scheme (36) we obtain Table 2.

First, we observe optimal convergence in each polynomial degree. In addition, we notice that, while the $L_{2}$ errors in both the standard and invariant scheme are comparable, the standard scheme has slightly smaller errors in the $L_{2}$ norm. Regardless, we observe that the invariant scheme is exact at the nodes, indicating that by preserving the Lie point symmetries we exactly preserve the flow of the solution at the nodes. This phenomenon originates from the fact that the transformation group (19) is of the same form as the exact solution, and thereby allows for the exact preservation of the solution through exactly preserving the Lie point symmetries.

Now, consider the case where (18) models an exponential growth problem with initial conditions

[TABLE]

and exact solution

[TABLE]

As both the invariant and non-invariant schemes are nonlinear, we note that this situation is difficult to simulate, as the exponential growth of the solution may prevent our nonlinear solver from converging for large time $t$. This issue may be overcome by decreasing the step size $\tau{}$. However, as $T$ increases linearly $\tau{}$ must decrease exponentially. Regardless, when $q=0$ and $\tau{}=0.25$ we obtain Figure 1. We observe that the invariant scheme is not only exact at the nodes, but that the error does not grow significantly over time. Conversely, we notice that the error of the standard scheme increases over time.

Example 4.4** **(Schwarzian differential equation).

We now compare the standard finite element scheme (54) against the invariantised approximation (59) for the third order ODE (50) considered in Example 3.2. In our numerical study we simulate the case where

[TABLE]

subject to the initial conditions

[TABLE]

up to the end time $T=1000$. Such numerical simulation approximates the exact solution

[TABLE]

The standard finite element approximation (54) leads to Table 4 while the invariant approximation (59) gives the results in Table 4.

For this problem, we observe that the maximal nodal error is comparable, however, the $L_{2}$ errors of the invariantised scheme are significantly smaller in the $q=0$ case highlighting a significant improvement through the invariantisation procedure. On the other hand, for higher polynomial degree, the $L_{2}$ errors are comparable for both schemes. Also, we note that in the $L_{2}$ norm, both schemes convergence optimally in each polynomial degree.

Example 4.5** **(A second order quasi-linear ODE).

We now consider the standard and invariant finite element schemes (65) and (70), which approximate the second order quasi-linear ODE (60). We enforce the initial data

[TABLE]

so that the exact solution is

[TABLE]

and simulate the solution over the domain $t\in[1,1000]$, where $T=1000$. Simulating the standard finite element approximation (65) for various polynomial degree we obtain Table 6, and for the invariant approximation (70), we obtain Table 6.

We observe that the errors for the standard and invariant scheme are comparable, with a slight improvement in the errors in the case $q=0$. These comparable results are a consequence of the exact solution behaving linearly for large $t$, as linear behaviour may be captured exactly by both of the schemes.

Example 4.6** **(**An ODE with non-projectable

action).**

The ODE

[TABLE]

has the simple solution given by

[TABLE]

Both the standard and invariant finite element approximations experimentally achieve best approximability, i.e., they exactly reconstruct polynomials of order $q{+}1$, as has been verified by the authors numerically. This means that, in practice, even for the case $q=0$ both finite element approximations are exact. As such, making a comparison of their respective errors is an exercise in futility. While the solution we are approximating is linear, both finite element approximations are nonlinear. That is to say that the approximations must be linearised and iteratively solved. In our experiments we utilised a Newton solver, which for sufficiently large time steps may not converge. In Table 7 we have tabulated a list of both standard and invariant schemes highlighting at which time step size their respective Newton solvers fail to converge. We observe that the symmetry-preserving scheme allows for the problem to be solved with much larger time steps compared to its non-invariant counterpart.

Example 4.7** **(**An illustrative example highlighting the benefits of

preserving Lie point symmetries for a naive finite element discretisation).**

Excluding Example 4.3, where the benefits of invariantisation are clear, we note that the numerical experiments conducted thus far do not highlight overwhelming benefits of the invariantisation procedure. We conjecture that this is due to the natural preservation of geometric structures by the standard finite element method, as discussed in Remark 2.5. Invariantising a less natural finite element approximation would indeed highlight the benefits of preserving Lie point symmetries, however, the authors feel this would make for an unfair comparison. For the sake of completeness, we now consider a poor finite element discretisation and show that its invariantisation recovers accurate long term simulations. We shall construct our naive approximation through the linearisation of an ODE model. Such an approximation may appear to be pathological, however, within the study of PDEs, linearisation is a common tool to obtain reduced models. This numerical example is motivated by the invariantisation of PDEs, which this work aims to inform.

Consider the initial value problem

[TABLE]

for some given constants $y_{0},y_{1}$. The ODE is invariant under the special linear group action

[TABLE]

Defining ${u_{0}}=y$, we may rewrite (100) as the first order system

[TABLE]

This system of ODEs is invariant under the extended group action

[TABLE]

Now, instead of employing the finite element formulation (17), we now consider the linearisation of (100), which of course for a fundamentally nonlinear problem we expect will yield poor results. Thus, we formulate a finite element approximation of (102) by seeking ${U_{0}},{U_{1}}\in\mathbb{V}_{1}^{C}(I_{n})$ such that

[TABLE]

This linearisation has been designed intentionally poorly, and is inconsistent with the solution of (102). Note that here we have fixed the degree of the finite element approximation for simplicity of exposition. Higher order finite element approximations would yield similar results. In particular, as our test space is $\mathbb{V}_{0}(I_{n})$, i.e., the space of piece-wise constant functions, we note that the group action does not alter the test functions.

Applying the group action (103) to the naive finite element scheme (104) we obtain

[TABLE]

which shows that the finite element formulation is not invariant. To obtain a symmetry-preserving formulation, we now implement the invariantisation process. Assuming ${U_{0}}\neq 0$, we choose the cross-section

[TABLE]

and obtain the moving frame

[TABLE]

Substituting the group normalisations (107) into (105) we obtain the invariant finite element scheme consisting of seeking ${U_{0}},{U_{1}}\in\mathbb{V}_{1}^{C}(I_{n})$ such that

[TABLE]

We observe that by preserving the Lie point symmetries of (100), we recover a consistent numerical scheme, even though consistency of the numerical scheme was destroyed by the initial naive discretisation. To illustrate the recovery of this long time structure we conduct a brief numerical experiment. Setting $\tau{}=0.01$ and using the initial conditions

[TABLE]

corresponding to the exact solution

[TABLE]

we display the point-wise errors for the naive discretisation (104) and the invariant discretisation (108) in Figure 2.

Similar to all previous examples in this section, we note that the invariant scheme (108) converges optimally, whereas the naive discretisation lacks consistency.

5. A comparison of symmetry-preserving finite element methodologies

The preservation of Lie point symmetries for finite element methods using the theory of moving frames has already been the subject of a preliminary investigation in [8]. In this final section, we highlight certain differences between the methodology presented in this paper and the one considered in [8].

First, the equation used to implement the finite element method is evidently different. In [8], the finite element method is applied to the original equation (1) while in the current work we consider the system of first order equations (7). Obviously, the two formulations represent the same equation, but the implementation of the finite element method differs. Indeed, when working with the original equation (1), it will be necessary to consider interpolating polynomials of degree greater than $m$ (or interpolating functions with smoother boundary conditions such as Hermite polynomials). Therefore, as the order of the equation increases, the degree of the interpolating polynomials will increase and the level of computational difficulty will follow the same trend. On the other hand, by recasting a differential equation as a system of first order equations, one can always work with low order interpolating polynomials. Though, it is worth mentioning that as the order of the original ODE increases, the number of auxiliary variables $\bm{U}$ will also increase, which introduces it own set of computational challenges.

The second, and most important, distinction between the two approaches is in the interpretation of the finite element functions approximating the solution to the differential equation, and how the symmetry group acts on the interpolating functions. In [8], finite element methods are viewed in terms of their underlying difference discretisations. As such, if

[TABLE]

is an approximation of the exact solution $y(t)$ to the differential equation (1), the induced action on $Y(t)$ was defined as the combination of the product action on the coefficients $Y_{n}$, i.e. $\widehat{Y_{n}}=g\cdot Y_{n}$, together with the action on the basis functions $\widehat{\phi_{n}}$ introduced in (24) to give

[TABLE]

In this approach, time derivatives were approximated using finite differences and moving frames were constructed by normalising these discrete approximations, resulting in what are known as discrete moving frames, [37]. For example, for the working example (8) considered in Examples 2.3, 2.6, 2.14 and 2.16, instead of considering the continuous cross-section (32), in [8] they introduce the discrete cross-section

[TABLE]

where the second term in the cross-section is the centred difference approximation of the first derivative $y_{t}$. We note that since $Y(t)$ is linear and the second term in the cross-section spans two elements, this cross-section is not equivalent to (32), although the two coincides in the continuous limit. One of the issues of working with this approach is that when considering higher order finite element schemes, one would need to specify cross-sections utilising more accurate difference quotients to obtain high order invariant schemes. On the other hand, in our proposed framework this issue does not arise.

In contrast, in the present work we utilise the continuous framework, and in particular, the continuously differentiable nature of the finite element approximation when restricted to the interior of a single element. As such, the action on $\bm{U}(t)$ is defined to be the same continuous action as the one acting on the exact solution $u(t)$. This significantly simplifies the procedure and allows us to generate invariant schemes for arbitrary degree approximations (assuming the action in $t$ is linear). However, there are also certain drawbacks to our methodology. In particular, as we do not evaluate the integrals of our approximation, we may not remove time dependent factors multiplying every term in the finite element approximation as in Examples 2.6 and 3.2. As such, we find ourselves multiplying the entirety of our approximation by some factor to ensure it is invariant under the integral. Such a problem would not occur if considering the underlying finite difference equation as in [8].

Another noticeable difference is in the group of transformations considered. In [8], to preserve the form of the approximation (111), only projectable group actions were considered. On the other hand, the current methodology applies to general Lie point symmetries.

Finally, we note that despite these differences, the two approaches have shown that in certain cases the symmetry-preserving schemes obtained can provide better long term numerical results than their non-invariant counterparts. Understanding when this happens and the varying benefits for different differential equations remains an open question.