Invariant Physics-Informed Neural Networks for Ordinary Differential
  Equations
Physics-informed neural networks have emerged as a prominent new method for
solving differential equations. While conceptually straightforward, they often
suffer training difficulties that lead to relatively large discretization
errors or the failure to obtain correct solutions. In this paper we introduce
invariant physics-informed neural networks for ordinary differential equations
that admit a finite-dimensional group of Lie point symmetries. Using the method
of equivariant moving frames, a differential equation is invariantized to
obtain a, generally, simpler equation in the space of differential invariants.
A solution to the invariantized equation is then mapped back to a solution of
the original differential equation by solving the reconstruction equations for
the left moving frame. The invariantized differential equation together with
the reconstruction equations are solved using a physics-informed neural
network, and form what we call an invariant physics-informed neural network. We
illustrate the method with several examples, all of which considerably
outperform standard non-invariant physics-informed neural networks.
Shivam Arora, Alex Bihlo, Francis Valiquette
# Invariant Physics-Informed Neural Networks for Ordinary Differential Equations

###### Abstract

Physics-informed neural networks have emerged as a prominent new method for solving differential equations. While conceptually straightforward, they often suffer training difficulties that lead to relatively large discretization errors or the failure to obtain correct solutions. In this paper we introduce _invariant physics-informed neural networks_ for ordinary differential equations that admit a finite-dimensional group of Lie point symmetries. Using the method of equivariant moving frames, a differential equation is invariantized to obtain a, generally, simpler equation in the space of differential invariants. A solution to the invariantized equation is then mapped back to a solution of the original differential equation by solving the reconstruction equations for the left moving frame. The invariantized differential equation together with the reconstruction equations are solved using a physcis-informed neural network, and form what we call an invariant physics-informed neural network. We illustrate the method with several examples, all of which considerably outperform standard non-invariant physics-informed neural networks.

## 1 Introduction

Physics-informed neural networks (PINN) are an emerging method for solving differential equations using deep learning, [23, 30]. The main idea behind this method is to train a neural network as an approximate solution interpolant for a system of differential equations. This is done by minimizing a loss function that incorporates both the differential equations and any associated initial and/or
boundary conditions. The method has a particular elegance as the derivatives in the differential equations can be computed using automatic differentiation rather than numerical discretization, which greatly simplifies the solution procedure, especially when solving differential equations on arbitrary surfaces, [35].

The ease of the discretization procedure in physics-informed neural networks, however, comes at the price of numerous training difficulties, and numerical solutions that are either not particularly accurate, or fail to converge at all to the true solution of the given differential equation. Since training a physics-informed neural network constitutes a non-convex optimization problem, an analysis of failure modes when physics-informed neural networks fail to train accurately is a non-trivial endeavour. This is why several modified training methodologies have been proposed, which include domain decomposition strategies, [18], modified loss functions, [38], and custom optimization, [3]. While all of these strategies, sometimes substantially, improve upon vanilla physics-informed neural networks, none of these modified approaches completely overcome all the inherent training difficulties.

Here we propose a new approach for training physics-informed neural networks, which relies on using Lie point symmetries of differential equations and the method of equivariant moving frames to simplify the form of the differential equations that have to be solved. This is accomplished by first projecting the differential equation onto the space of differential invariant to produce an _invariantized differential equation_. The solution to the invariantized equation is then mapped back to the solution of the original equation by solving a system of first order differential equations for the left moving frame, called _reconstruction equations_. The invariant physics-informed neural network architecture proposed in this paper consists of simultaneously solving the system of equations consisting of the invariantized differential equation and the reconstruction equations using a physics-informed neural network. The method proposed is entirely algorithmic, and can be implemented for any system of differential equations that is strongly invariant under the action of a group of Lie point symmetries. Since almost all equations of physical relevance admit a non-trivial group of Lie point symmetries, the proposed method is potentially a viable path for improving physics-informed neural networks for many real-world applications. The idea of projecting a differential equation into the space of invariants and then reconstructing its solution is reminiscent of the recent work [37], where the authors consider Hamiltonian systems with symmetries, although the tools used in our paper and in [37] to achieve the desired goals are very different. Moreover, in our approach we do not assume that our equations have an underlying symplectic structure.

To simplify the theoretical exposition, we focus on the case of ordinary differential equations in this paper. We show using several examples that the proposed approach substantially improves upon the numerical results achievable with vanilla physics-informed neural networks. Applications to partial differential equations will be considered elsewhere.

The paper is organized as follows. We first review relevant work on physics-informed neural networks and symmetry preserving numerical methods in Section 2. In Section 3 we introduce the methods of equivariant moving frames and review how it can be used to solve ordinary differential equations that admit a group of Lie point symmetries. Building on Section 3 we introduce a version of invariant physics-informed neural network in Section 4. We illustrate our method with several examples in Section 5. The examples show that our proposed invariant physics-informed neural
network formulation can yield better numerical results than its non-invariant version. A short summary and discussion about potential future research avenues concludes the paper in Section 6.

## 2 Previous work

Physics-informed neural networks were first proposed in [23], and later popularized through the work of Raissi et al., [30]. The main idea behind physics-informed neural networks is to train a deep neural network to directly approximate the solution to a system of differential equations. This is done by defining a loss function that incorporates the given system of equations, along with any relevant initial and/or boundary conditions. Crucially, this turns training physics-informed neural networks into a multi-task, non-convex optimization problem that can be challenging to minimize, [22]. There have been several solutions proposed to overcome the training difficulties and improve the generalization capabilities of physics-informed neural networks. These include modified loss functions, [26; 38], meta-learned optimization, [3], domain decomposition methods, [6; 18], and the use of operator-based methods, [11; 24].

The concepts of symmetries and transformation groups have also received considerable attention in the machine learning community. Notably, the equivariance of convolutional operations with respect to spatial translations has been identified as a crucial ingredient for the success of convolutional neural networks, [13]. The generalization of this observation for other types of layers of neural networks and other transformation groups has become a prolific subfield of deep learning since. For example, see [16] for some recent results.

Here we do not consider the problem of endowing a neural network with equivariance properties but rather investigate the question whether a better formulation of a given differential equation can help physics-informed neural networks better learn a solution. As we will be using symmetries of differential equations for this re-formulation, our approach falls within the framework of geometric numerical integration. The problem of symmetry-preserving numerical schemes, in other words the problem of designing discretization methods for differential equations that preserve the symmetries of a given differential equation, has been studied extensively over the past several decades, see [14; 34; 39] for some early work on the topic. Invariant discretization schemes have since been proposed for finite difference, finite volume, finite elements and meshless methods, [2; 4; 5; 7; 8; 12; 19; 28; 31; 32].

## 3 Method

In this section we introduce the theoretical foundations on which the invariant physics-informed neural network framework is based. In order to fix some notation, we begin by recalling certain well-known results pertaining to symmetries of differential equations, and refer the reader to [9; 10; 17; 27] for a more thorough exposition. Within this field, the use of moving frames to solve differential equations admitting symmetries is not as well-known. Therefore, the main purpose of this section is to introduce this solution procedure. In contrast to the approach proposed in [25], we avoid the introduction of computational variables. All computations are based on the differential invariants of the prolonged group action which results in less differential equations. Our approach is a simplified version of the algorithm presented in [36], which deals with partial differential equations admitting infinite-dimensional symmetry Lie pseudo-groups.
As mentioned in the introduction, in this paper we limit ourselves to the case of ordinary differential equations.

### Invariant differential equations

Let \(M\) be a \((q+1)\)-dimensional manifold with \(q\geq 1\). Given a one-dimensional curve \(C\subset M\), we introduce the local coordinates \(z=(t,u)=(t,u^{1},\ldots,u^{q})\) on \(M\) so that the curve \(C\) is locally specified by the graph of a function \(C=\{(t,f(t))\}\). Accordingly, the \(n\)-th order jet space \(\mathrm{J}^{(n)}\) is locally parametrized by \(z^{(n)}=(t,u^{(n)})\), where \(u^{(n)}\) denotes all the derivatives \(u^{\alpha}_{j}=u^{\alpha}_{t^{j}}\) of order \(0\leq j\leq n\), with \(\alpha=1,\ldots,q\).

Let \(G\) be an \(r\)-dimensional Lie group (locally) acting on \(M\):

\[(T,U)=Z=g\cdot z=g\cdot(t,u),\qquad\text{where}\qquad g\in G. \tag{1}\]

The group transformation (1) induces an action on curves \(C\subset M\), which prolongs to the jet space \(\mathrm{J}^{(n)}\):

\[Z^{(n)}=g\cdot z^{(n)}. \tag{2}\]

Coordinate expressions for the prolonged action (2) are obtained by applying the implicit total derivative operator

\[\mathrm{D}_{T}=\frac{1}{\mathrm{D}_{t}(T)}\,\mathrm{D}_{t},\qquad\text{where} \qquad\mathrm{D}_{t}=\frac{\partial}{\partial t}+\sum_{j=0}^{\infty}\sum_{ \alpha=1}^{q}\,u^{\alpha}_{j+1}\frac{\partial}{\partial u^{\alpha}_{j}}\]

denotes the standard total derivative operator, to the transformed dependent variables \(U^{\alpha}\):

\[U^{\alpha}_{j}=U^{\alpha}_{T^{j}}=\mathrm{D}^{j}_{T}(U^{\alpha}),\qquad\alpha= 1,\ldots,q,\qquad j\geq 0. \tag{3}\]

In the following we use the notation \(\Delta(z^{(n)})=\Delta(t,u^{(n)})=0\) to denote a system of differential equations, and use the index notation \(\Delta_{i}(z^{(n)})=0\), \(i=1,\ldots,l\), to label each equation in \(\Delta(z^{(n)})\). Also, a differential equation \(\Delta(z^{(n)})=0\), can either be a single equation or represent a system of differential equations.

**Definition 1**.: A nondegenerate1 ordinary differential equation \(\Delta(z^{(n)})=0\) is said to be _strongly invariant_ under the prolonged action of a connected local Lie group of transformations \(G\) if and only if

Footnote 1: A differential equation is nondegenerate if at every point in its solution space it is both locally solvable and of maximal rank, [27].

\[\Delta(g\cdot z^{(n)})=0\qquad\text{for all}\qquad g\in G\]

near the identity element.

**Remark 2**.: Strong invariance is more restrictive than the usual notion of symmetry, where invariance is only required to hold on the solution space. In the following, we require strong invariance to guarantee that our differential equation is an invariant function.

Invariance is usually stated in terms of the infinitesimal generators of the group action. To this end, let

\[\mathbf{v}_{\kappa}=\xi_{\kappa}(t,u)\frac{\partial}{\partial t}+\sum_{\alpha= 1}^{q}\phi^{\alpha}_{\kappa}(t,u)\frac{\partial}{\partial u^{\alpha}},\qquad \kappa=1,\ldots,r, \tag{4}\]
be a basis of infinitesimal generators for the group action \(G\). The prolongation of the vector fields (4), induced from the prolonged action (3), is given by

\[\mathbf{v}^{(n)}=\xi_{\kappa}(t,u)\frac{\partial}{\partial t}+\sum_{j=0}^{n}\sum _{\alpha=1}^{q}\phi_{\kappa}^{\alpha,j}(t,u^{(j)})\frac{\partial}{\partial u_ {j}^{\alpha}},\qquad\kappa=1,\ldots,r,\]

where the prolonged coefficients are computed using the standard prolongation formula

\[\phi_{\kappa}^{\alpha,j}=\mathrm{D}_{t}^{j}(\phi_{\kappa}^{\alpha}-\xi_{ \kappa}u_{1}^{\alpha})+\xi_{\kappa}u_{j+1}^{\alpha},\qquad\kappa=1,\ldots,r, \qquad\alpha=1,\ldots,q,\qquad 0\leq j\leq n.\]

**Proposition 3**.: A nondegenerate ordinary differential equation \(\Delta(z^{(n)})=0\) is strongly invariant under the prolonged action of a connected local Lie group of transformations \(G\) if and only if

\[\mathbf{v}_{\kappa}^{(n)}[\Delta_{i}(z^{(n)})]=0,\qquad\kappa=1,\ldots,r, \qquad i=1,\ldots,l,\]

where \(\mathbf{v}_{1},\ldots,\mathbf{v}_{r}\) is a basis of infinitesimal generators for the group of transformations \(G\).

**Remark 4**.: As one may observe, we do not include the initial conditions

\[u^{(n-1)}(t_{0})=u_{0}^{(n-1)} \tag{5}\]

when discussing the symmetry of the differential equation \(\Delta(z^{(n)})=\Delta(x,u^{(n)})=0\). This is customary when studying symmetries of differential equations. Of course, the initial conditions are necessary to select a particular solution and when implementing numerical simulations.

### Invariantization

Given a nondegenerate differential equation \(\Delta(z^{(n)})=0\) strongly invariant under the prolonged action of an \(r\)-dimensional Lie group \(G\) acting regularly on \(\mathrm{J}^{(n)}\), we now explain how to use the method of equivariant moving frames to "project" the differential equation onto the space of differential invariants. For the theoretical foundations of the method of equivariant moving frames, we refer the reader to the foundational papers [15, 21] and the textbook [25].

**Definition 5**.: A Lie group \(G\) acting smoothly on a \(\mathrm{J}^{(n)}\) is said to act _freely_ if the isotropy group

\[G_{z^{(n)}}=\{g\in G\,|\,g\cdot z^{(n)}=z^{(n)}\}=\{e\}\]

at the point \(z^{(n)}\) is trivial for all \(z^{(n)}\in\mathrm{J}^{(n)}\). The Lie group \(G\) is said to act _locally freely_ if \(G_{z^{(n)}}\) is a discrete subgroup of \(G\) for all \(z^{(n)}\in\mathrm{J}^{(n)}\).

**Remark 6**.: More generally we can restrict Definition 5, and the subsequent considerations, to a \(G\)-invariant submanifold \(\mathcal{V}^{(n)}\subset\mathrm{J}^{(n)}\). To simplify the discussion. we assume \(\mathcal{V}^{(n)}=\mathrm{J}^{(n)}\).

**Definition 7**.: A _right moving frame_ is a \(G\)-equivariant map \(\rho\colon\mathrm{J}^{(n)}\to G\) such that

\[\rho(g\cdot z^{(n)})=g\cdot\rho(z^{(n)})\]

for all \(g\in G\) where the prolonged action is defined. Taking the group inverse of a right moving frame yields the left moving frame

\[\overline{\rho}(z^{(n)})=\rho(z^{(n)})^{-1}\]

satisfying the equivariance condition

\[\overline{\rho}(g\cdot z^{(n)})=g\cdot\overline{\rho}(z^{(n)}).\]
**Theorem 8**.: A moving frame exists in the neighborhood of a point \(z^{(n)}\in\mathrm{J}^{(n)}\) provided the prolonged action of \(G\) on \(\mathrm{J}^{(n)}\) is (locally) free and regular.

A moving frame is obtained by selecting a cross-section \(\mathcal{K}\subset\mathrm{J}^{(n)}\) to the orbits of the prolonged action. Keeping with most applications, and to simplify the exposition, assume \(\mathcal{K}\) is a coordinate cross-section obtained by setting \(r\) coordinates of the jet \(z^{(n)}\) to constant values:

\[z^{a_{\kappa}}=c^{\kappa},\qquad\kappa=1,\dots,r. \tag{6}\]

Solving the normalization equations

\[Z^{a_{\kappa}}=c^{\kappa},\qquad\kappa=1,\dots,r,\]

for the group parameters, yields a right moving frame \(\rho\). Given a moving frame, there is a systemic procedure for constructing differential invariant functions.

**Definition 9**.: Let \(\rho\colon\mathrm{J}^{(n)}\to\mathbb{R}\) be a right moving frame. The invariantization of the differential function \(F\colon\mathrm{J}^{(n)}\to\mathbb{R}\) is the differential invariant function

\[\iota(F)(z^{(n)})=F(\rho(z^{(n)})\cdot z^{(n)}).\]

In particular, invariantization of the coordinate jet functions

\[\iota(z^{(n)})=\rho(z^{(n)})\cdot z_{n}\]

yields differential invariants that can be used as coordinates on the cross-section \(\mathcal{K}\). In particular, the invariantization of the coordinates used to define the cross-section in (6) are constant

\[\iota(z^{a_{\kappa}})=c^{\kappa},\qquad\kappa=1,\dots,r,\]

and are called _phantom invariants_. The remaining invariantized coordinates are called _normalized invariants_. In light of Theorem 5.32 in [29], assume there are \(q+1\) normalized invariants

\[H,I^{1},\,\dots,I^{q}, \tag{7}\]

such that locally the invariants

\[I^{\alpha}=I^{\alpha}(H),\qquad\alpha=1,\dots,q,\]

are independent functions of the invariant \(H\), and generate the algebra of differential invariants. This means that any differential invariant can be expressed in terms of (7) and their invariant derivatives with respect to \(\mathrm{D}_{H}\). In the following we let \(I^{(n)}\) denote the derivatives of \(I=(I^{1},\dots,I^{q})\) with respect to \(H\), up to order \(n\).

Assuming the differential equation \(\Delta(z^{(n)})=0\) is strongly invariant and its solutions are transverse to the prolonged action, this equation, once invariantized, will yield a differential equation in the space of invariants

\[\iota[\Delta(t,u^{(n)})]=\Delta_{\text{Inv}}(H,I^{(k)})=0,\qquad\text{where} \qquad k\leq n. \tag{8}\]

Initial conditions for (8) are obtained by invariantizing (5) to obtain

\[I^{(k-1)}(H_{0})=I_{0}^{(k-1)}. \tag{9}\]
**Example 10**.: To illustrate the concepts introduced thus far, we use the Schwarz equation

\[\frac{u_{ttt}}{u_{t}}-\frac{3}{2}\bigg{(}\frac{u_{tt}}{u_{t}}\bigg{)}^{2}=F(t) \tag{10}\]

as our running example. This equation admits a three-dimensional Lie group of point transformations given by

\[T=t,\qquad U=\frac{\alpha u+\beta}{\gamma u+\delta},\qquad\text{where}\qquad g =\begin{bmatrix}\alpha&\beta\\ \gamma&\delta\end{bmatrix}\in\text{SL}(2,\mathbb{R}), \tag{11}\]

so that \(\alpha\delta-\beta\gamma=1\). A cross-section to the prolonged action

\[U_{T} =\text{D}_{t}(U)=\frac{u_{t}}{(\gamma u+\delta)^{2}},\] \[U_{TT} =\text{D}_{t}(U_{T})=\frac{u_{tt}}{(\gamma u+\delta)^{2}}-\frac{2 \gamma u_{t}^{2}}{(\gamma u+\delta)^{3}}, \tag{12}\] \[U_{TTT} =\text{D}_{t}(U_{TT})=\frac{u_{ttt}}{(\gamma u+\delta)^{2}}- \frac{6\gamma u_{t}u_{tt}}{(\gamma u+\delta)^{3}}+\frac{6\gamma^{2}u_{t}^{3}}{ (\gamma u+\delta)^{4}},\]

is given by

\[\mathcal{K}=\{u=0,\,u_{t}=\sigma,\,u_{tt}=0\}\subset\mathcal{V}^{(n)}\subset \text{J}^{(n)}, \tag{13}\]

where \(\sigma=\text{sign}(u_{t})\), and \(\mathcal{V}^{(n)}=\{z\in\text{J}^{(n)}\,|\,u_{t}\neq 0\}\) with \(n\geq 2\). Solving the normalization equations

\[U=0,\qquad U_{T}=\sigma,\qquad U_{TT}=0, \tag{14}\]

together with the unitary constraint \(\alpha\delta-\beta\gamma=1\), we obtain the right moving frame

\[\alpha=\pm\frac{1}{\sqrt{|u_{t}|}},\qquad\beta=\mp\frac{u}{\sqrt{|u_{t}|}}, \qquad\gamma=\pm\frac{u_{tt}}{2|u_{t}|^{3/2}},\qquad\delta=\pm\frac{2u_{t}^{2 }-uu_{tt}}{2|u_{t}|^{3/2}}, \tag{15}\]

where the sign ambiguity comes from solving the normalization \(U_{T}=\sigma\), which involves the quadratic term \((\gamma u+\delta)^{2}\). Invariantizing the third order derivative \(u_{ttt}\) produces the differential invariant

\[\iota(u_{ttt})=\frac{u_{ttt}}{(\gamma u+\delta)^{2}}-\frac{6\gamma u_{t}u_{tt }}{(\gamma u+\delta)^{3}}+\frac{6\gamma^{2}u_{t}^{3}}{(\gamma u+\delta)^{4}} \bigg{|}_{(15)}=\sigma\bigg{(}\frac{u_{ttt}}{u_{t}}-\frac{3}{2}\bigg{(}\frac{ u_{tt}}{u_{t}}\bigg{)}^{2}\bigg{)}. \tag{16}\]

In terms of the general theory previously introduced, we have the invariants

\[H=t,\qquad I=\frac{u_{ttt}}{u_{t}}-\frac{3}{2}\bigg{(}\frac{u_{tt}}{u_{t}} \bigg{)}^{2}. \tag{17}\]

Since the independent variable \(t\) is an invariant, instead of using \(H\), we use \(t\) in the following computations. The invariantization of Schwarz equation (10) then yields the algebraic equation

\[I=F(t). \tag{18}\]

Since the prolonged action is transitive on the fibers of each component \(\{(t,u,u_{t},u_{tt})\,|\,u_{t}>0\}\cup\{(t,u,u_{t},u_{tt})\,|\,u_{t}<0\}= \mathcal{V}^{(2)}\), any initial conditions

\[u(t_{0})=u_{0},\qquad u_{t}(t_{0})=u_{t}^{0},\qquad u_{tt}(t_{0})=u_{tt}^{0},\]

is mapped, under invariantization, to the identities

\[0=0,\qquad\sigma=\sigma,\qquad 0=0.\]
### Recurrence relations

Using the recurrence relations we now explain how the invariantized equation (8) can be derived symbolically, without requiring the coordinate expressions for the moving frame \(\rho\) or the invariants \((H,I)\). The key observation is that the invariantization map \(\iota\) and the exterior differential do not, in general, commute

\[\iota\circ\mathrm{d}\neq\mathrm{d}\circ\iota.\]

The extend by which these two operations do not commute is encapsulated in the recurrence relations. To state these equations we need to introduce the (contact) invariant one-form

\[\varpi=\iota(\mathrm{d}t)=\rho^{*}(\mathrm{D}_{t}(T))\,\mathrm{d}t,\]

which comes from invariantizing the horizontal one-form \(\mathrm{d}t\), see [21] for more details.

Given a Lie group \(G\), let \(g\in G\) be represented by a faithful matrix. Then the _right Maurer-Cartan form_ is given by

\[\mu=\mathrm{d}g\cdot g^{-1}. \tag{19}\]

The pull-back of the Maurer-Cartan form (19) by a right moving frame \(\rho\) yields the invariant matrix

\[\nu=\mathrm{d}\rho\cdot\rho^{-1}=\begin{bmatrix}I_{ij}\end{bmatrix}\varpi, \tag{20}\]

where the invariants \(I_{ij}\) are called Maurer-Cartan invariants.

**Proposition 11**.: Let \(F\colon\mathrm{J}^{(n)}\to\mathbb{R}\) be a differential function. The recurrence relation for the invariantization map \(\iota\) is

\[\mathrm{d}[\iota(F)]=\iota[\mathrm{d}F]+\sum_{\kappa=1}^{r}\iota[\mathbf{v}_{ \kappa}^{(n)}(F)]\,\nu^{\kappa}, \tag{21}\]

where \(\nu^{1},\ldots,\nu^{r}\) is a basis of normalized Maurer-Cartan forms extracted from (20).

Substituting for \(F\) in (21) the jet coordinates (6) specifying the coordinate cross-section \(\mathcal{K}\) leads to \(r\) linear equations for the normalized Maurer-Cartan forms \(\nu^{1},\ldots,\nu^{r}\). Solving those equations and substituting the result back in (21) yields a symbolic expression for the differential of any invariantized differential function \(F\), without requiring the coordinate expressions for the moving frame \(\rho\).

**Example 12**.: Continuing Example 10, a basis of infinitesimal generators for the group action (11) is provided by

\[\mathbf{v}_{1}=\frac{\partial}{\partial u},\qquad\mathbf{v}_{2}=u\frac{ \partial}{\partial u},\qquad\mathbf{v}_{3}=u^{2}\frac{\partial}{\partial u}.\]

The prolongation of those vector fields, up to order 2, is given by

\[\mathbf{v}_{1}^{(3)} =\frac{\partial}{\partial u},\] \[\mathbf{v}_{2}^{(3)} =u\frac{\partial}{\partial u}+u_{t}\frac{\partial}{\partial u_{t} }+u_{tt}\frac{\partial}{\partial u_{tt}},\] \[\mathbf{v}_{3}^{(3)} =u^{2}\frac{\partial}{\partial u}+2uu_{t}\frac{\partial}{\partial u _{t}}+2(u_{t}^{2}+uu_{tt})\frac{\partial}{\partial u_{tt}}.\]
Applying the recurrence relation (21) to \(t\), \(u\), \(u_{t}\), and \(u_{tt}\) yields

\[\begin{split}\mathrm{d}[\iota(t)]&=\varpi,\\ \mathrm{d}[\iota(u)]&=\iota(u_{t})\varphi+\nu^{1},\\ \mathrm{d}[\iota(u_{t})]&=\iota(u_{tt})\varpi+ \iota(u_{t})\nu^{2}+2\iota(u)\iota(u_{t})\nu^{3},\\ \mathrm{d}[\iota(u_{tt})]&=\iota(u_{ttt})\varpi+ \iota(u_{tt})\nu^{2}+2[\iota(u_{t})^{2}+\iota(u)\iota(u_{tt})]\nu^{3}.\end{split} \tag{22}\]

Recalling the cross-section (13) and the invariants (16), (17), we make the substitutions \(H=\iota(t)=t\), \(\iota(u)=0\), \(\iota(u_{t})=\sigma\), \(\iota(u_{tt})=0\), \(\iota(u_{tt})=\sigma I\) into (22) and obtain

\[\mathrm{d}t=\varpi,\qquad 0=\varpi+\nu^{1},\qquad 0=\nu^{2},\qquad 0=I\sigma \varpi+2\nu^{3}.\]

Solving for the normalized Maurer-Cartan forms yields

\[\nu^{1}=-\sigma\,\varpi,\qquad\nu^{2}=0,\qquad\nu^{3}=-\frac{I\sigma}{2}\varpi.\]

In matrix form we have that

\[\nu=\begin{bmatrix}0&-\sigma\\ \frac{1}{2}\sigma I&0\end{bmatrix}\varpi=\begin{bmatrix}0&-\sigma\\ \frac{1}{2}\sigma F(t)&0\end{bmatrix}\varpi,\]

where we used the algebraic relationship (18) originating from the invariance of Schwarz' equation (10).

### Reconstruction

Let \(I(H)\) be a solution to the invariantized differential equation (8) with initial conditions (9). In this section we explain how to reconstruct the solution to the original equation \(\Delta(x,u^{(n)})=0\) with initial conditions (5). To do so, we introduce the reconstruction equations for the left moving frame \(\overline{\rho}=\rho^{-1}\):

\[\mathrm{d}\overline{\rho}=-\overline{\rho}\cdot\mathrm{d}\rho\cdot\overline{ \rho}=-\overline{\rho}\,\nu, \tag{23}\]

where \(\nu\) is the normalized Maurer-Cartan form introduced in (20). As we have seen in Section 3.3, the invariantized Maurer-Cartan matrix \(\nu\) can be obtained symbolically using the recurrence relations for the phantom invariants. Since \(\nu\) is invariant, it can be expressed in terms of \(H\), the solution \(I(H)\), and its derivatives. Thus, equation (23) yield a first order system of differential equations for the group parameters expressed in the independent variable \(H\). Integrating (23), we obtain the left moving frame that sends the invariant curve \((H,I(H))\) to the original solution

\[(t(H),u(H))=\overline{\rho}(H)\cdot\iota(t,u)(H). \tag{24}\]

Assuming \(t_{H}>0\), the initial conditions to the reconstruction equations (23) are given by

\[\overline{\rho}(H_{0})=\overline{\rho}_{0}\qquad\text{such that}\qquad \overline{\rho}_{0}\cdot\iota(t_{0},u_{0}^{(n-1)})=(t_{0},u_{0}^{(n-1)}). \tag{25}\]

If \(t_{H}<0\), one can always reparametrize the solution so that the derivative becomes positive.

The solution (24) is a parametric curve with the invariant \(H\) serving as the parameter. From a numerical perspective, this is sufficient to graph the solution. Though we note that by inverting \(t=t(H)\) to express the invariant \(H=H(t)\) in terms of \(t\), we can recover the solution as a function of \(t\):

\[u=u(H(t)).\]
**Example 13**.: The left moving frame

\[\overline{\rho}=\begin{bmatrix}\alpha&\beta\\ \gamma&\delta\end{bmatrix}\quad\in\quad\text{SL}(2,\mathbb{R})\]

that will send the solution (18) to the original function \(u(t)\) is a solution to the reconstruction equations

\[\begin{bmatrix}\alpha_{t}&\beta_{t}\\ \gamma_{t}&\delta_{t}\end{bmatrix}=-\begin{bmatrix}\alpha&\beta\\ \gamma&\delta\end{bmatrix}\begin{bmatrix}0&-\sigma\\ \frac{1}{2}\sigma F(t)&0\end{bmatrix}=\begin{bmatrix}\alpha&\beta\\ \gamma&\delta\end{bmatrix}\begin{bmatrix}0&\sigma\\ -\frac{1}{2}\sigma F(t)&0\end{bmatrix},\] (26a) with initial conditions \[\delta_{0}=\pm\frac{1}{\sqrt{|u_{t}^{0}|}},\qquad\beta_{0}=\pm\frac{u_{0}}{ \sqrt{|u_{t}^{0}|}},\qquad\gamma_{0}=\mp\frac{u_{tt}^{0}}{2(u_{t}^{0})^{3/2}}, \qquad\alpha_{0}=\pm\sqrt{|u_{t}^{0}|}\mp\frac{u_{0}u_{tt}^{0}}{2(|u_{t}^{0}|) ^{3/2}}. \tag{26b}\]

Then, the solution to Schwarz' equation (10) is

\[u(t)=\overline{\rho}\cdot 0=\frac{\beta}{\delta}. \tag{27}\]

### Summary

Let us summarize the algorithm for solving an ordinary differential equation \(\Delta(t,u^{(n)})=0\) admitting a group of Lie point transformations \(G\) using the method of moving frames.

1. Select a cross-section \(\mathcal{K}\) to the prolonged action.
2. Choose \(q+1\) invariants \(H\), \(I^{1},\ldots,I^{q}\) from \(\iota(t,u^{(n)})\), that generate the algebra of differential invariants, and assume \(I^{1}(H),\ldots,I^{q}(H)\) are functions of \(H\).
3. Invariantize the differential equation \(\Delta(t,u^{(n)})=0\) and use the recurrence relation (21) to write the result in terms of \(H\) and \(I^{(k)}\) to obtain the equation \(\Delta_{\text{Inv}}(H,I^{(k)})=0\).
4. Solve the equation \(\Delta_{\text{Inv}}(H,I^{(k)})=0\) subject to the initial conditions (9).
5. A parametric solution to the original equation \(\Delta(t,u^{(n)})=0\) is given by \(\overline{\rho}(H)\cdot\iota(t,u)(H)\), where the left moving frame \(\overline{\rho}(H)\) is a solution of the reconstruction equation (23) subject to the initial conditions (25).

## 4 Invariant physics-informed neural networks

Before introducing our invariant physics-informed neural network, we recall the definition of the standard physics-informed loss function that needs to be minimized when solving ordinary differential equations. To this end, assume we want to solve the ordinary differential equation \(\Delta(t,u^{(n)})=0\) subject to the initial conditions \(u^{(n-1)}(t_{0})=u_{0}^{(n-1)}\) on the interval \([t_{0},t_{f}]\).

First we introduce the collocation points \(\{t_{i}\}_{i=0}^{\ell}\) sampled randomly over the interval \([t_{0},t_{f}]\) with \(t_{0}=t_{0}<t_{1}<\cdots<t_{\ell}=t_{f}\). Then, a neural network of the form \(u_{\boldsymbol{\theta}}(t)=\mathcal{N}_{\boldsymbol{\theta}}(t)\), parameterized
by the parameter vector \(\mathbf{\theta}\), is trained to approximate the solution of the differential equation, i.e. \(u_{\mathbf{\theta}}(t)\approx u(t)\), by minimizing the physics-informed loss function

\[\mathcal{L}(\mathbf{\theta})=\mathcal{L}_{\Delta}(\mathbf{\theta})+\alpha\,\mathcal{L}_ {\text{I.C.}}(\mathbf{\theta})\] (28a) with respect to \[\mathbf{\theta}\], where \[\mathcal{L}_{\Delta}(\mathbf{\theta})=\sum_{i=0}^{\ell}\ \left[\Delta(t_{i},u_{\mathbf{\theta}}^{( n)}(t_{i}))\right]^{2}\] (28b) is the _differential equation loss_, \[\mathcal{L}_{\text{I.C.}}(\mathbf{\theta})=\left[u_{\mathbf{\theta}}^{(n-1)}(t_{0})-u_{ 0}^{(n-1)}\right]^{2} \tag{28c}\]

is the _initial condition loss_, and \(\alpha\) is a hyper-parameter to re-scale the importance of both loss functions. We note that the differential equation loss is the mean squared error of the differential equation evaluated at the collocation points \(\left\{t_{i}\right\}_{i=0}^{\ell}\subset[t_{0},t_{f}]\) over which the numerical solution is sought. The initial condition loss is likewise the mean squared error between the true initial conditions and the initial conditions approximated by the neural network. We note in passage that the initial conditions could alternatively be enforced as a hard constraint in the neural network, [11, 23], in which case the physics-informed loss function would reduce to \(\mathcal{L}_{\Delta}(\mathbf{\theta})\) only.

The physics-informed loss function (28) is minimized using gradient descent, usually using the Adam optimizer, [20], but also more elaborate optimizers can be employed, [3]. The particular elegance of the method of physics-informed neural networks lies in the fact that the derivatives \(u_{\mathbf{\theta}}^{(n)}\) of the neural network solution approximation are computed using _automatic differentiation_, [1], which is built into all modern deep learning frameworks such as JAX, TensorFlow, or PyTorch.

Similar to the above standard physics-informed neural network, an invariant physics-informed neural network is a feed-forward neural network approximating the solution of the invariantized

Figure 1: Solving a differential equation using moving frames.
differential equation and the reconstruction equations for the left moving frame. In other words, the loss function to be minimized is defined using symmetries of the given differential equation.

In light of the five step process given in Section 3.5, assume the invariantized equation \(\Delta_{\text{Inv}}(H,I^{(k)})=0\) and the reconstruction equation \(\mathrm{d}\overline{\rho}=-\overline{\rho}\nu\) have been derived. Introduce an interval of integration \([H_{0},H_{f}]\) over which the numerical solution is sought, and consider the collocation points \(\{H_{i}\}_{i=0}^{\ell}\subset[H_{0},H_{f}]\), such that \(H_{\ell}=H_{f}\). The neural network has to learn a mapping between \(H\) and the functions

\[I_{\boldsymbol{\theta}}(H)\qquad\text{and}\qquad\overline{\rho}_{\boldsymbol{ \theta}}(H),\]

where \(I_{\boldsymbol{\theta}}(H)\) denotes the neural network approximation of the differential invariants \(I(H)\) solving (8), and \(\overline{\rho}_{\boldsymbol{\theta}}(H)\) is the approximation of the left moving frame \(\overline{\rho}(H)\) solving the reconstruction equations (23). We note that the output size of the network depends on the numbers of invariants \(I(H)\) and the size of the symmetry group via \(\overline{\rho}(H)\).

The network is trained by minimizing the invariant physics-informed loss function consisting of the invariantized differential equation loss and the reconstruction equations loss defined as the mean squared error

\[\mathcal{L}_{\Delta_{\text{Inv}},\overline{\rho}}(\boldsymbol{ \theta})=\sum_{i=0}^{\ell}\big{(}\big{[}\Delta_{\text{Inv}}(H_{i},I_{ \boldsymbol{\theta}}^{(k)}(H_{i}))\big{]}^{2}+\big{[}\mathrm{d}\overline{ \rho}_{\boldsymbol{\theta}}(H_{i})+\overline{\rho}_{\boldsymbol{\theta}}(H_{ i})\,\nu(H_{i},I_{\boldsymbol{\theta}}^{(\kappa)}(H_{i}))\big{]}^{2}\big{)}. \tag{29}\]

We supplement the loss function (29) with the initial conditions (9) and (25) by considering the invariant initial conditions loss function

\[\mathcal{L}_{\text{I.C.}}(\boldsymbol{\theta})=\big{[}I_{\boldsymbol{\theta}} ^{(k-1)}(H_{0})-I_{0}^{(k-1)}\big{]}^{2}+\big{[}\overline{\rho}_{\boldsymbol{ \theta}}(H_{0})-\overline{\rho}_{0}\big{]}^{2}.\]

The final invariant physics-informed loss function is thus given by

\[\mathcal{L}_{\text{Inv}}(\boldsymbol{\theta})=\mathcal{L}_{\Delta_{\text{Inv} },\overline{\rho}}(\boldsymbol{\theta})+\alpha\,\mathcal{L}_{\text{I.C.}}( \boldsymbol{\theta}),\]

where \(\alpha\) is again a hyper-parameter rescaling the importance of the equation and initial condition losses.

## 5 Examples

We now implement the invariant neural network architecture introduced in Section 4 for several examples. We also train a standard physics-informed neural network to compare the solutions obtained. For both models we use feed-forward neural networks minimizing the invariant loss function and standard PINN loss function, respectively. For the sake of consistency, all networks used throughout this section have 5 layers, with 40 nodes per layer, and use the hyperbolic tangent as activation functions. For most examples, the loss stabilizes at fewer than 5,000 epochs, but for uniformity we trained all models for 5,000 epochs. The numerical errors of the two neural network solutions are obtained by comparing the numerical solutions to the exact solution, if available, or to the numerical solution obtained using odeint in scipy.integrate. We also compute the mean square error over the entire interval of integration for all examples together with the standard deviation averaged over \(5\) runs. These results are summarized in Table 1. Finally, the point-wise square error plots for each example are provided to show the error varying over the interval of integration.
**Example 14**.: As our first example, we consider the Schwarz equation (10), with \(F(t)=2\). For the numerical simulations, we used the initial conditions

\[u_{0}=u_{tt}^{0}=0,\qquad u_{t}^{0}=1. \tag{30}\]

According to (18) the invariantization of Schwarz' equation yields the algebraic constraint \(I=2\). Thus, the loss function (29) will only contain the reconstruction equations (26a). Namely,

\[\alpha_{t}+\beta=\beta_{t}-\alpha=\gamma_{t}+\delta=\delta_{t}-\gamma=0,\]

where we used the fact that \(\sigma=1\). Substituting (30) into (26b), yields the initial conditions

\[\delta_{0}=\alpha_{0}=\pm 1,\qquad\beta_{0}=\gamma_{0}=0\]

for the reconstruction equations. In our numerical simulations we worked with the positive sign. Once the reconstruction equations have been solved, the solution to the Schwarz equation is given by the ratio (27). The solution is integrated on the interval \(t\in[0,\pi]\). Error plots for the solutions obtained via the invariant PINN and the standard PINN implementations are given in Figure 2. These errors are obtained by comparing the numerical solutions to the exact solution \(u(t)=\tan(t)\). Clearly, the invariant implementation is substantially more precise near the vertical asymptote at \(t=\pi/2\).

**Example 15**.: As our second example, we consider the logistic equation

\[u_{t}=u(1-u) \tag{31}\]

occurring in population growth modeling. Equation (31) admits the one-parameter symmetry group

\[T=t,\qquad U=\frac{u}{1+\epsilon\,ue^{-t}},\qquad\text{where}\qquad\epsilon \in\mathbb{R}.\]

Implementing the algorithm outlined in Section 3.5, we choose the cross-section \(\mathcal{K}=\{u=1\}\). This yields the invariantized equation

\[I=\iota(u_{t})=0.\]

Figure 2: Time series of the squared error for the Schwarz equation (10).
The reconstruction equation is

\[\epsilon_{t}=I=0,\]

subject to the initial condition

\[\epsilon(t_{0})=\bigg{(}\frac{1-u_{0}}{u_{0}}\bigg{)}e^{t_{0}},\]

where \(u_{0}=0.5\) and our interval of integration is \([0,\pi]\). The solution to the logistic equation is then given by

\[u(t)=\frac{1}{1+\epsilon\,e^{-t}}.\]

As Figure 3 illustrates, the error incurred by the invariant PINN model is significantly smaller, by about a factor of more than 100, than the standard PINN implementation when compared to the exact solution \(u(t)=1/(1+e^{-t})\).

**Example 16**.: We now consider the driven harmonic oscillator

\[u_{tt}+u=\sin(t^{a}), \tag{32}\]

which appears in inductor-capacitor circuits, [33]. In the following we set \(a=0.99\), which yields bounded solutions close to resonance occurring when \(a=1\). The differential equation (32) admits the two-dimensional symmetry group of transformations

\[T=t,\qquad U=u+\alpha\sin(t)+\beta\cos(t),\qquad\text{where}\qquad\alpha, \beta\in\mathbb{R}.\]

A cross-section to the prolonged action is given by \(\mathcal{K}=\{u=u_{t}=0\}\). The invariantization of (32) yields

\[I=\iota(u_{tt})=\sin(t^{a}).\]

The reconstruction equations are

\[\alpha_{t}=\sin(t^{a})\cos(t),\qquad\beta_{t}=-\sin(t^{a})\sin(t), \tag{33}\]

Figure 3: Time series of the squared error for the logistic equation (31).
with initial conditions

\[\alpha(t_{0})=u_{0}\sin(t_{0})+u_{t}^{0}\cos(t_{0}),\qquad\beta(t_{0})=u_{0}\cos( t_{0})-u_{t}^{0}\sin(t_{0}),\]

where, in our numerical simulations, we set \(u_{0}=u_{t}^{0}=1\) and integrate over the interval \([0,10]\). Given a solution to the reconstruction equations (33), the solution to the driven harmonic oscillator (32) is

\[u(t)=\alpha(t)\,\sin(t)+\beta(t)\cos(t).\]

Figure 4 shows the error for the invariant PINN implementation and the standard PINN approach compared to the solution obtained using the odeint Runge-Kutta method in scipy.integrate. As in the previous two examples, the invariant version yields substantially better numerical results than the standard PINN method.

**Example 17**.: We now consider the second order ordinary differential equation

\[u_{tt}=\exp{[-u_{t}]} \tag{34}\]

with an exponential term. Equation (34) admits a three-dimensional symmetry group action given by

\[T=e^{\epsilon}t+a,\qquad U=e^{\epsilon}u+\epsilon\,e^{\epsilon}x+b,\]

where \(a,b,\epsilon\in\mathbb{R}\). In the following, we only consider the one-dimensional group

\[T=e^{\epsilon}t,\qquad U=e^{\epsilon}u+\epsilon\,e^{\epsilon}x.\]

We note that in this example the independent variable is not invariant as in the previous examples. A cross-section to the prolonged action is given by \(\mathcal{K}=\{u_{x}=0\}\). Introducing the invariants

\[H=\ln{\left[\frac{1}{1-\iota(t)}\right]}=\ln{\left[\frac{1}{1-tu_{t}}\right]},\qquad I=\iota(u)=\exp[-u_{t}](u-tu_{t}),\]

Figure 4: Time series of the squared error for the driven harmonic oscillator (32).
the invariantization of the differential equation (34) reduces to the first order linear equation

\[I_{H}+I=e^{-H}-1.\]

The reconstruction equation for the left moving frame is simply

\[\epsilon_{H}=1.\]

In terms of \(\epsilon\) and \(I\), the parametric solution to the original differential equation (34) is

\[t=e^{\epsilon}(1-e^{-H}),\qquad u=e^{\epsilon}(I+\epsilon(1-e^{-H})). \tag{35}\]

The solution to (34) is known and is given by

\[u(t)=(t+c_{1})\ln(t+c_{1})-t+c_{2}, \tag{36}\]

where \(c_{1}\), \(c_{2}\) are two integration constants. For the numerical simulations, we use the initial conditions

\[I_{0}=\exp[-u_{t}^{0}](u_{0}-t_{0}u_{t}^{0}),\qquad\epsilon_{0}=u_{t}^{0},\]

where \(u_{0}=u(t_{0})\), \(u_{t}^{0}=u_{t}(t_{0})\) with \(t_{0}=0\), and \(c_{1}=\exp(-5)\), \(c_{2}=0\) in (36). The interval of integration \([H_{0},H_{f}]\) is given by

\[H_{0}=\ln\bigg{[}\frac{1}{1-t_{0}u_{t}^{0}}\bigg{]},\qquad H_{f}=\ln\bigg{[} \frac{1}{1-t_{f}u_{t}^{f}}\bigg{]}, \tag{37}\]

where \(u_{t}^{f}=u_{t}(t_{f})\) and \(t_{f}=2\). We choose the interval of integration given by (37), so that when \(t\) is given by (35) it lies in the interval \([0,2]\).

Figure 5 shows the error obtained for the invariant PINN model when compared to the exact solution (34), and similarly for the non-invariant PINN model. As in all previous examples, the invariant version drastically outperforms the standard PINN approach.

**Example 18**.: As our final example, we consider a system of first order ODEs

\[u_{t}=-u+(t+1)v,\qquad v_{t}=u-tv. \tag{38}\]

This system admits a two-dimensional symmetry group of transformations given by

\[T=t,\qquad U=\alpha u+\beta t,\qquad V=\alpha v+\beta,\]

where \(\alpha>0\) and \(\beta\in\mathbb{R}\). Working with the cross-section \(\mathcal{K}=\{u=1,v=0\}\), the invariantization of (38) yields

\[I=\iota(u_{x})=-1,\qquad J=\iota(v_{x})=1.\]

The reconstruction equations are

\[\alpha_{t}=\alpha(1+t),\qquad\beta_{t}=\alpha\]

subject to the initial conditions \(\alpha_{0}=1\), \(\beta_{0}=1\), corresponding to the initial conditions \(u_{0}=v_{0}=1\), when \(t_{0}=0\). In our numerical simulations we integrated over the interval \([0,2]\). The solution to (38) is then given by

\[u(t)=\alpha(t)+t\,\beta(t),\qquad v(t)=\beta(t).\]
As in all previous examples, comparing the numerical solutions to the exact solution

\[u(t)=\sqrt{\frac{2}{\pi}}\,c\,e^{-(t+1)^{2}/2}+c\,t\,\text{erf}\!\left(\frac{t+1} {\sqrt{2}}\right)+kt,\qquad v(t)=c\,\text{erf}\!\left(\frac{t+1}{\sqrt{2}} \right)+k,\]

with \(c=\left(\sqrt{2/\pi}\,\exp(-1/2)\right)^{-1}\) and \(k=1-c\,\text{erf}(1/\sqrt{2})\), where \(\text{erf}(t)=2/\sqrt{\pi}\int_{0}^{t}e^{-x^{2}}\mathrm{d}x\) is the standard error function, we observe in Figure 6 that the invariant version of the PINN model considerably outperforms its non-invariant counterpart.

## 6 Summary and conclusions

In this paper we have introduced the notion of invariant physics-informed neural networks. These combine physics-informed neural networks with methods from the group analysis of differential

Figure 5: Time series of the squared error for the exponential equation (34).

Figure 6: Time series of the squared error for the system of equations (38).
equations to simplify the form of the differential equations that have to be solved. In turns, this simplifies the loss function that has to be minimized, and our numerical tests show that the solutions obtained with the invariant model outperformed their non-invariant counterparts, and typically considerably so. Table 1 summarizes the examples considered in the paper and shows that the invariant PINN outperforms vanilla PINN for the examples considered.

The proposed method is fully algorithmic and as such can be applied to any system of differential equations that is strongly invariant under the prolonged action of a group of Lie point symmetries. It is worth noting that the work proposed here parallels some of the work on invariant discretization schemes which, for ordinary differential equations, also routinely outperform their non-invariant counterparts. We have observed this to also be the case for physics-informed neural networks.

Lastly, while we have restricted ourselves here to the case of ordinary differential equations, our method extends to partial differential equations as well. Though, when considering partial differential equations, it is not sufficient to project the equations onto the space of differential invariants as done in this paper. As explained in [36], integrability conditions among the differential invariants must also be added to the invariantized differential equations. In the multivariate case, the reconstruction equations (23) will then form a system of first order partial derivatives for the left moving frame. Apart from these modifications, invariant physics-informed neural networks can also be constructed for partial differential equations, which will be investigated elsewhere.

#### Acknowledgments and Disclosure of Funding

This research was undertaken, in part, thanks to funding from the Canada Research Chairs program and the NSERC Discovery Grant program. The authors also acknowledge support from the Atlantic Association for Research in the Mathematical Sciences (AARMS) Collaborative Research Group on _Mathematical Foundations for Scientific Machine Learning_.
