Direct Poisson neural networks: learning non-symplectic mechanical systems

A Poisson bracket on a manifold $\mathcal{M}$ (physically corresponding to the state space, for instance position and momentum of the body) is a skew-symmetric bilinear algebra on the space $\mathcal{F}(\mathcal{M})$ of smooth functions on $\mathcal{M}$ given by

$$\begin{align} \left\{\bullet,\bullet\right\}:\mathcal{F}\left(\mathcal{M}\right)\times\mathcal{F}\left(\mathcal{M}\right)\rightarrow\mathcal{F}\left(\mathcal{M}\right). \end{align} \tag{ 2 }$$

Poisson brackets satisfy the Leibniz rule

$$\begin{align} \left\{F,H G\right\} = \left\{F,H\right\}G + H\left\{F,G\right\}, \end{align} \tag{ 3 }$$

and the Jacobi identity,

$$\begin{align} \left\{F,\left\{H,G\right\}\right\} + \left\{H,\left\{G,F\right\}\right\} + \left\{G,\left\{F,H\right\}\right\} = 0, \end{align} \tag{ 4 }$$

for arbitrary functions F, H and G, [12, 20, 21]. A manifold equipped with a Poisson bracket is called a Poisson manifold and is denoted by a two-tuple $(\mathcal{M},\{\bullet,\bullet\})$. A function C is called a Casimir function if it commutes with all other functions that is $\{F,C\} = 0$ for all F. For instance, the magnitude of the angular momentum of an RB is a Casimir function.

Hamiltonian dynamics. Hamiltonian dynamics can be seen as evolution on a Poisson manifold. For a Hamiltonian function (physically energy) H on $\mathcal{M}$, the Hamiltonian vector field and Hamilton's equation are

$$\begin{align} X_H\left(F\right): = \left\{F,H\right\}, \qquad \dot{\mathbf{x}} = X_H\left(\mathbf{x}\right) = \left\{\mathbf{x},H\right\}, \end{align} \tag{ 5 }$$

respectively, where $\mathbf{x} \in \mathcal{M}$ is a parametrization of manifold $\mathcal{M}$. The algebraic properties of the Poisson bracket have some physical consequences. Skew-symmetry implies energy conservation,

$$\begin{align} \dot{H} = \left\{H,H\right\} = 0, \end{align} \tag{ 6 }$$

while the Leibniz rule ensures that the dynamics does not depend on biasing the energy by a constant. Referring to a Poisson bracket, one may determine the Poisson bivector field according to

$$\begin{align} L\left(\mathrm{d}F,\mathrm{d}H\right): = \left\{F,H\right\}. \end{align} \tag{ 7 }$$

This identification makes it possible to define a Poisson manifold by a tuple $(\mathcal{M},L)$ consisting of a manifold and a Poisson bivector. In this notation, the Jacobi identity can be rewritten as $\mathcal{L}_{\mathbf{X}_H} L = 0$, that is the Lie derivative of the Poisson bivector with respect to the Hamiltonian vector field is zero [25]. In other words, the Jacobi identity expresses the self-consistency of the Hamiltonian dynamics in the sense that both the building blocks (Hamiltonian function and the bivector field) are constant along the evolution.

Assuming a local coordinate system $\mathbf{x} = (x^i)$ on $\mathcal{M}$, Poisson bivector determines Poisson matrix $L = [L^{kl}]$ which enables us to write [26]

$$\begin{align} L = L^{kl}\frac{\partial}{\partial x^k}\wedge \frac{\partial}{\partial x^l}. \end{align} \tag{ 8 }$$

In this realization, the Poisson bracket and Hamilton's equations are written as

$$\begin{align} \left\{F,H\right\} = L^{kl}\frac{\partial F}{\partial x^k}\frac{\partial H}{\partial x^l}, \qquad \dot{x}^k = L^{kl}\frac{\partial H}{\partial x^l}, \end{align} \tag{ 9 }$$

respectively. Here, we have assumed summation over the repeated indices. Further, the Jacobi identity (4) turns out to be the following system of partial differential equations

$$\begin{align} J^{ijl} = L^{kl}\frac{\partial L^{ij}}{\partial x^k} +L^{ki}\frac{\partial L^{jl}}{\partial x^k} +L^{kj}\frac{\partial L^{li}}{\partial x^k} = 0. \end{align} \tag{ 10 }$$

The left-hand side of this equation is called Jacobiator, and in the case of Hamiltonian systems, it is equal to zero. Jacobi identity (10) is a system of differential equations consisting of PDEs. In 3D systems (for instance the RB, that has three degrees of freedom), Jacobi identity (10) is a single PDE whose general solution is known [27]. In four-dimensional (4D) systems (for instance a particle moving in a plane, or the Shivamoggi equations), Jacobi identity (10) consists of four PDEs, and there are some partial results, but for an arbitrary n, according to our knowledge, there is no general solution yet. We shall focus on 3D, 4D, and six-dimensional (6D) (HT or particle in a volume) cases in the upcoming subsections.

Symplectic manifolds. If there is no non-constant Casimir function for a Poisson manifold, then it is also a symplectic manifold. Although we can see symplectic manifolds as examples of Poisson manifolds, it is possible to define a symplectic manifold in a direct way without referring to a Poisson manifold. A manifold $\mathcal{M}$ is called symplectic if it is equipped with a closed non-degenerate two-form (called a symplectic two-form) Ω. A two-form is called non-degenerate when

$$\begin{align} \Omega \left( X,Y \right) = 0, \qquad \forall X\in \mathfrak{X}\left(\mathcal{M}\right) \end{align} \tag{ 11 }$$

implies Y = 0. A two-form is closed when being in the kernel of de Rham exterior derivative, $\mathrm{d}\Omega = 0$. A Hamiltonian vector field X_H on a symplectic manifold $(\mathcal{M},\Omega)$ is defined as

$$\begin{align} \iota_{X_{H}}\Omega =\mathrm{d}H, \end{align} \tag{ 12 }$$

where ι is the contraction operator (more precisely the interior derivative) [25]. Referring to a symplectic manifold one can define a Poisson bracket

$$\begin{align} \left\{F,H\right\}: = \Omega\left(X_F,X_H\right), \end{align} \tag{ 13 }$$

where the Hamiltonian vector fields are defined through equation (12). The closedness of the symplectic two-form Ω guarantees the Jacobi identity (4). The non-degeneracy condition of Ω puts an extra condition to the bracket (13) that the Casimir functions are only the constant functions, in contrast with Poisson manifolds, which may have also non-constant Casimirs. The Darboux–Weinstein coordinates show more explicitly the relationship between Poisson and symplectic manifolds in a local picture.

Darboux–Weinstein coordinates. We start with $n = (2m+k)$-dimensional Poisson manifold $\mathcal{M}$ equipped with Poisson bivector L. Near every point of the Poisson manifold, the Darboux–Weinstein coordinates $(x^i) = (q^a,p_b,u^\alpha)$ (here a runs from 1 to m, and α runs from 1 to k) give a local form of the Poisson bivector

$$\begin{align} L = \frac{\partial }{\partial q^a}\wedge \frac{\partial }{\partial p_a} + \frac{1}{2}\lambda ^{\alpha\beta} \frac{\partial }{\partial u^\alpha}\wedge \frac{\partial }{\partial u^\beta} \end{align} \tag{ 14 }$$

with the coefficient functions $\lambda^{\alpha\beta}$ equal zero at the origin. If k = 0 in the local formula (14), then there remains only the first term on the right-hand side and the Poisson manifold turns out to be a 2m-dimensional symplectic manifold. Newtonian mechanics, for instance, fits this kind of realization. On the other hand, if m = 0 in (14), then there remains only the second term which is a full-degenerate Poisson bivector. A large class of Poisson manifolds is of this form, namely Lie–Poisson structure on the dual of a Lie algebra including RB dynamics, Vlasov dynamics, etc. In general, Poisson bivectors have both the symplectic part as well as the fully degenerate part, for instance, the HT dynamics in section 2.4.

When the Poisson bivector is non-degenerate, it generates a symplectic Poisson bracket, and it commutes with the canonical Poisson bivector

$$\begin{align} \mathbf{L}_\mathrm{can} = \begin{pmatrix} \mathbf{0} & \mathbf{I}\\ -\mathbf{I} & \mathbf{0}\end{pmatrix} \end{align} \tag{ 15 }$$

in the sense that

$$\begin{align} \mathbf{L}\cdot\mathbf{L}_\mathrm{can} - \mathbf{L}_\mathrm{can}\cdot\mathbf{L} = \mathbf{0}. \end{align} \tag{ 16 }$$

This compatibility condition is employed later in section 3 to measure the error of DPNNs when learning symplectic Poisson bivectors.

In this subsection, we focus on 3D Poisson manifolds, following [28–30]. One of the important observations in 3D is the isomorphism between the space of vectors and the space of skew-symmetric matrices given by

$$\begin{equation} \mathbf{L} = \begin{pmatrix} 0 & -J_z & J_y\\ J_z & 0 & -J_x \\ -J_y & J_x & 0\end{pmatrix} \leftrightarrow \mathbf{J} = \left(J_x, J_y, J_z\right). \end{equation} \tag{ 17 }$$

This isomorphism lets us write Jacobi identity (10) as a single scalar equation

$$\begin{align} \mathbf{J}\cdot \left(\nabla \times \mathbf{J}\right) = 0, \end{align} \tag{ 18 }$$

see, for example, [29, 31–33]. The general solution of Jacobi identity (18) is

$$\begin{align} \mathbf{J} = \frac{1}{\phi}\nabla C \end{align} \tag{ 19 }$$

for arbitrary functions φ and C, where C is a Casimir function. Hamilton's equation then takes the particular form

$$\begin{align} \dot{\mathbf{x}} = \mathbf{J}\times \nabla H = \frac{1}{\phi}\nabla C \times \nabla H . \end{align} \tag{ 20 }$$

Note that by changing the roles of the Hamiltonian function H and the Casimir C one can arrive at another Hamiltonian structure for the same system. In this case, the Poisson vector is defined as $\mathbf{J} = -(1/\phi)\nabla H$ and the Hamiltonian function is C. This is an example of a bi-Hamiltonian system, that manifests integrability [27, 34–36]. A bi-Hamiltonian system admits two different but compatible Hamilton formulations. In 3D, two Poisson vectors, say J₁ and J₂ are compatible if

$$\begin{align} \mathbf{J}_1 \cdot \left(\nabla\times \mathbf{J}_2\right) = \mathbf{J}_2 \cdot \left(\nabla\times \mathbf{J}_1\right). \end{align} \tag{ 21 }$$

Indeed, if a 3D Hamiltonian system admits two Poisson bivectors (represented by two vector J₁ and J₂) and two Hamiltonians H₁ and H₂,

$$\begin{align} \dot{\mathbf{x}} = \mathbf{J}_1 \times \nabla H_1 = \mathbf{J}_2 \times \nabla H_2, \end{align} \tag{ 22 }$$

then both the Hamiltonians are conserved and thus H₂ serves as a Casimir function in the formulation with J₁ while H₁ being a Casimir for J₂. Compatibility condition (21) then follows from the solution of Jacobi identity (19). This compatibility condition will be used later in section 3 to measure the error of learning Poisson bivectors in 3D by DPNNs.

RB dynamics. Let us consider an example of a 3D Hamiltonian system, a freely rotating RB. The state variable $\mathbf{M}\in\mathcal{M}$ is the angular momentum in the frame of reference co-rotating with the RB. The Poisson structure is

$$\begin{align} \left\{F,H\right\}^\mathrm{(RB)}\left(\mathbf{M}\right) = - \mathbf{M} \cdot \frac{\partial F}{\partial \mathbf{M}}\times\frac{\partial H}{\partial \mathbf{M}}, \end{align} \tag{ 23 }$$

see [37]. Poisson bracket (23) is degenerate because it preserves any function of the magnitude of M. The Hamiltonian function is the energy

$$\begin{align} H = \frac{1}{2}\left(\frac{M_x^2}{I_x}+\frac{M_y^2}{I_y}+\frac{M_z^2}{I_z}\right), \end{align} \tag{ 24 }$$

where I_x , I_y and I_z are moments of inertia of the body. In this case, Hamilton's equation

$$\begin{align} \dot{\mathbf{M}} = \mathbf{M} \times \frac{\partial H}{\partial \mathbf{M}} \end{align} \tag{ 25 }$$

gives Euler's RB equation [38].

In 4D, we consider the following local coordinates $(u,x,y,z) = \left( u,\mathbf{x} \right)$. A skew-symmetric matrix L can be identified with a couple of vectors $\mathbf{U} = ( U^{1},U^{2},U^{3})$ and $\mathbf{V} = ( V^{1},V^{2},V^{3} )$ as

$$\begin{align} L = \begin{pmatrix} 0 & -U^{1} & -U^{2} & -U^{3} \\ U^{1} & 0 & -V^{3} & V^{2} \\ U^{2} & V^{3} & 0 & -V^{1} \\ U^{3} & -V^{2} & V^{1} & 0 \end{pmatrix}. \end{align} \tag{ 26 }$$

After this identification, Jacobi identity (10) turns out to be a system of PDEs consisting of four equations [39]

$$\begin{align} \partial_{u}\left(\mathbf{U}\cdot\mathbf{V}\right) & = \mathbf{V}\cdot\left( \partial_{u} \mathbf{U}-\nabla\times\mathbf{V}\right) , \end{align} \tag{ 27a }$$

$$\begin{align} \nabla\left(\mathbf{U}\cdot\mathbf{V}\right) & = \mathbf{V}\left(\nabla\cdot\mathbf{U}\right)- \mathbf{U}\times\left( \partial_{u}\mathbf{U}-\nabla\times\mathbf{V}\right) . \end{align} \tag{ 27b }$$

Note that L is degenerate (its determinant being zero) if and only if $\mathbf{U}\cdot\mathbf{V} = 0$. So, for degenerate Poisson matrices, the Jacobi identity is satisfied if

$$\begin{align} \nabla\cdot\mathbf{U} = 0, \qquad \partial_u\mathbf{U} -\nabla\times\mathbf{V} = \mathbf{0} . \end{align} \tag{ 28 }$$

Superintegrability. Assume that a given dynamics admits two time-independent first integrals, say H₁ and H₂. Then, when vectors U and V have the form

$$\begin{align} \mathbf{U} = \nabla H_{1}\times\nabla H_{2}, \qquad \mathbf{V} = \partial_{u}H_{1}\nabla H_{2}-\partial_{u}H_{2}\nabla H_{1}, \end{align} \tag{ 29 }$$

they constitute a Poisson bivector, and in particular Jacobi identity (28) is satisfied. Functions H₁ and H₂ are Casimir functions. For a Hamiltonian function H₃, the Hamiltonian dynamics is

$$\begin{align} \dot{u} & = -\left( \nabla H_{1}\times\nabla H_{2}\right) \cdot\nabla H_{3}, \end{align} \tag{ 30a }$$

$$\begin{align} \mathbf{\dot{x}} & = \left(\nabla H_{1}\times\nabla H_{2}\right)\partial _{u}H_{3}+\left(\nabla H_{2}\times\nabla H_{3}\right)\partial_{u}H_{1}+\left(\nabla H_{3} \times\nabla H_{1}\right)\partial_{u}H_{2}. \end{align} \tag{ 30b }$$

By permuting the roles of the functions H₁, H₂, and H₃ (two of them are Casimirs and one of them is Hamiltonian), one arrives at two additional Poisson realizations of the dynamics. This is the case of a superintegrable (tri-Hamiltonian) formulation, [40].

Two Poisson bivectors of a superintegrable 4D Hamiltonian system must satisfy the compatibility condition

$$\begin{align} \mathbf{U}_1\cdot\mathbf{V}_2 + \mathbf{V}_1\cdot\mathbf{U}_2 = 0 \end{align} \tag{ 31 }$$

where $\mathbf{U}_{1,2}$ and $\mathbf{V}_{1,2}$ are the vectors identified from formula (26). We shall use this compatibility condition to measure the error of DPNNs when learning Poisson bivectors in 4D.

Shivamoggi equations. An example of a 4D Poisson (and non-symplectic) system is the Shivamoggi equations, which arise in the context of magnetohydrodynamics,

$$\begin{align} \dot{u} = -uy,\qquad \dot{x} = zy,\qquad \dot{y} = zx-u^{2},\qquad \dot{z} = xy, \end{align} \tag{ 32 }$$

see [41, 42]. The first integrals of this system of equations are

$$\begin{align} H_{1} = x^{2}-z^{2},\qquad H_{2} = z^{2}+u^{2}-y^{2},\qquad H_{3} = u\left(z+x\right). \end{align} \tag{ 33 }$$

Vectors U_i and V_i of Poisson matrices $N^{\left( i\right) }$ ($i = 1,2,3$) for the Hamiltonian functions H₁, H₂, and H₃ are

$$\begin{align} \mathbf{U}_{1} & = 2u\left( -y,z,y\right) ,\qquad \mathbf{V} _{1} = 2\left( u^{2},y\left(x+z\right),u^{2}-z\left(x+z\right)\right) , \nonumber\\ \mathbf{U}_{2} & = 2\left( x+z\right) \left( 0,u,0\right) ,\qquad \mathbf{V}_{2} = 2\left( x+z\right)\left( x,0,-z\right) , \nonumber\\ \mathbf{U}_{3} & = -4\left( yz,zx,xy\right) ,\qquad \mathbf{V} _{3} = 4u\left( -x,0,z\right) , \end{align} \tag{ 34 }$$

respectively. Note that all these three Poisson matrices are degenerate, since $\mathbf{U}_{i}\cdot\mathbf{V}_{i} = 0$ holds for all $i = 1,2,3$. The equations of motion can be written as

$$\begin{align*} X = \vartheta N^{\left( 1\right) }\bar{\nabla}H_{1} = \vartheta N^{\left( 2\right) }\bar{\nabla}H_{2} = \vartheta N^{\left( 3\right) }\bar{\nabla}H_{3}, \ \ \ \vartheta = -\frac{1}{4\left(x+z\right)} \end{align*}$$

up to multiplication with a conformal factor ϑ in all three cases. Note that the 4D gradient is denoted by $\bar{\nabla}H = (\partial_u H,\partial_x H,\partial_y H,\partial_z H)$.

6D Hamiltonian systems can again be symplectic or non-symplectic (degenerate). The former case is represented by a particle in 3D (P3D) while the latter is for instance the HT dynamics. Since the evolution of a P3D is canonical and thus analogical to the 2D dynamics, we shall recall only the HT dynamics.

A supported rotating RB in a uniform gravitational field is called HT [43]. The mechanical state of the body is described by the position of the center of mass r and angular momentum M. In this case, the Poisson bracket is

$$\begin{align} \left\{F,G\right\}^{\left(\mathrm{HT}\right)} \left(\mathbf{r},\mathbf{M}\right) = -\mathbf{M} \cdot \left(\frac{\partial F}{\partial \mathbf{M}}\times \frac{\partial G}{\partial \mathbf{M}}\right) - \mathbf{r} \cdot \left(\frac{\partial F}{\partial \mathbf{M}}\times \frac{\partial G}{\partial \mathbf{r}}-\frac{\partial G}{\partial \mathbf{M}}\times \frac{\partial F}{\partial \mathbf{r}}\right). \end{align} \tag{ 35 }$$

Even though the model is even dimensional, it is not symplectic. Two non-constant Casimir functions are r² and $\mathbf{M}\cdot\mathbf{r}$. In this case, we assume the Hamiltonian function as

$$\begin{align} H = \frac{1}{2}\left(\frac{M_x^2}{I_x}+\frac{M_y^2}{I_y}+\frac{M_z^2}{I_z}\right) + Mgl \mathbf{r}\cdot \boldsymbol{\chi}, \end{align} \tag{ 36 }$$

where $-g\boldsymbol{\chi}$ is the vector of gravitational acceleration. Hamilton's equation is then

$$\begin{align} \dot{\mathbf{M}} & = \mathbf{M} \times \frac{\partial H}{\partial \mathbf{M}} + \mathbf{r} \times \frac{\partial H}{\partial \mathbf{r}} \end{align} \tag{ 37a }$$

$$\begin{align} \dot{\mathbf{r}} & = - \mathbf{r} \times \frac{\partial H}{\partial \mathbf{M}}. \end{align} \tag{ 37b }$$

In the following sections, we apply DPNNs to the models here recalled and show that DPNNs are capable to extract the Poisson bivector and Hamiltonian from simulated trajectories of the models.

Neural networks serve as robust interpolators to high-dimensional functions. As they consist of several interconnected layers of neurons, they can be seen as compositions of high-dimensional mappings, which provides the robustness of approximation. Actually, the universal approximation theorem states roughly speaking that every continuous function from $\mathbb{R}^n$ to $\mathbb{R}^m$ can be approximated by a neural network with one hidden layer and a non-linear activation function [44–46]. These multi-layered architectures have shown a great capability for tasks like image recognition, natural language processing, and other applications of deep learning [47].

A prototype of such neural networks is the fully connected (or dense) network, see figure 1, where every neuron in a given layer is connected to every neuron in the subsequent layer. The mathematical operation each neuron typically carries out is an affine transformation of the inputs followed by the application of a non-linear activation function [48]. Parameters of the affine transformation (weights and biases) are then subject to the learning procedure, where they are sought so that a prescribed loss function is minimized.

Zoom In Zoom Out Reset image size

The affine transformation between two consecutive layers is associated with weights (the strengths of particular neuron–neuron connections) and biases, which are added to the weighted sum of the inputs, turning a zero input into a non-zero value. The activation function plays a crucial role because it provides non-linearity to the neural network, enabling it to approximate more complicated functions. In our approach, we use the softplus function,

$$\begin{align} \sigma_{\mathrm{softplus}}\left(x\right) = \log\left(1 + e^x\right), \end{align} \tag{ 38 }$$

which ensures a differentiability of the outputs [49, 50].

Training a neural network then involves tuning its weights and biases to minimize the difference between the actual and the expected outputs. This difference is quantified by a so called loss function. The weights and biases are then adapted in such a way that the loss function is minimized by a backpropagation algorithm, which computes the gradient of the loss function with respect to the parameters by applying the chain rule of calculus [51] and then shifts the parameters along the gradient. The simplest implementation called stochastic gradient descent (SGD) works by splitting a large dataset into smaller chunks—minibatches and then performing incremental updates along the loss gradient of the corresponding batch. Many variations of this algorithm exist such as the Adam algorithm (short for adaptive moment estimation) [52] that works like SGD but utilizes momentum term to make the optimization smoother. We shall use Adam optimizer in this work. Although much time could be spent on the various techniques of deep learning, we refer an interested reader to book [48], and proceed to the actual application in Hamiltonian systems.

When we have a collection of snapshots of a trajectory of a Hamiltonian system, how to identify the underlying mechanics? In other words, how to recognize the model (Poisson bivector and energy) from the snapshots? Machine learning provides a robust method for such task. Machine learning has previously been shown to reconstruct GENERIC models [13, 14, 53], but the Poisson bivector is typically known and symplectic. PNNs [24] provide a method for learning also non-symplectic mechanics, which however relies on the identification of dimension of the symplectic subdynamics in the Darboux–Weinstein coordinates and on a transformation to the coordinates. Here, we show a robust method that does not need to know the dimension of the symplectic subsystem and that satisfies Jacobi identity also in the coordinates in which the data are prescribed. Therefore, we refer to the method as DPNNs.

DPNNs learn Hamiltonian mechanics directly by training a model for the $\mathbf{L}(\mathbf{x})$ matrix and a model for the Hamiltonian $H(\mathbf{x})$ simultaneously. The neural network that encodes $\mathbf{L}(\mathbf{x})$ only learns the upper triangular part of L and skew-symmetry is then automatically satisfied. The network has one hidden fully connected layer (with 64 neurons) equipped with the softplus activation. The network that learns $H(\mathbf{x})$ has the same structure. The input layer has the dimension of the physical system (for instance three nodes in the case of RB), and dimension of the output layer is the number of independent components of the trained quantity (number of independent components of the L bivector, or one when learning H). The actual learning was implemented within the standard framework PyTorch [54], using the Adam optimizer [52], and is publically available on Github [55].

The loss function contains squares of deviation of the training data and predicted trajectories that are obtained by the implicit midpoint rule (IMR) numerically solving the exact equations (for the training data) or Hamilton's equation with the trained models for $\mathbf{L}(\mathbf{x})$ and $H(\mathbf{x})$ (for the predicted data). Although such a model leads to a good match between the grand truth trajectories and predicted trajectories, it does not need to satisfy Jacobi identity. Therefore, we use also an alternative model where squares of the Jacobiator (10) are added to the loss function, which enforces Jacobi identity in a soft way, see figure 2. Finally, in 3D we know the form of the Poisson bivector since we have the general solution of Jacobi identity (20). In such a case, the neural network encoding L can be simplified to a network learning $C(\mathbf{x})$ and Jacobi identity is automatically satisfied, see figure 3. In summary, we use three methods:

• (WJ) Training $\mathbf{L}(\mathbf{x})$ and $H(\mathbf{x})$ without the Jacobi identity. The loss function consists only of the L₂ error measuring the discrepancy between the training steps and steps obtained by IMR iterations with bivector L and Hamiltonian $H(\mathbf{x})$ encoded by the networks:
  $$\begin{align} \mathcal{L}_{\mathrm{WJ}} = c_{\mathrm{mov}} \sum_{n} |\left(\left(\mathbf{x}_{n+1}\right)_{\mathrm{exact}}-\mathbf{x}_n\right) - \left(\left(\mathbf{x}_{n+1}\right)_{\mathrm{NN}} - \mathbf{x}_n\right)|^2 \end{align} \tag{ 39a }$$
  where n goes over all training steps. Each point $(\mathbf{x}_{n+1})_{\mathrm{NN}}$ is calculated from x_n by the IMR method (using the L and H that are being trained and encoded by the neural networks). Prefactor $c_{\mathrm{mov}}$ is used to scale the loss function to numerically advantageous values (typically we use $c_{\mathrm{mov}} = 10$).
• (SJ) Training $\mathbf{L}(\mathbf{x})$ and $H(\mathbf{x})$ with soft Jacobi identity, where the L₂-norm of the Jacobiator (10) is a part of the loss function. The loss function is thus
  $$\begin{align} \mathcal{L}_{\mathrm{SJ}} = \mathcal{L}_{\mathrm{WJ}} + c_{\mathrm{Jac}}\sum_n \sum_{ijk} |J^{ijk}\left(\mathbf{x}_n\right)|^2, \end{align} \tag{ 39b }$$
  where i, j, and k go over the dimension of the system (here 3, 4, or 6). Prefactor $c_{\mathrm{Jac}}$ is used to scale the loss function into better numericalvalues (typically we use $c_{\mathrm{Jac}} = 10$).
• (IJ) Training $C(\mathbf{x})$ and $H(\mathbf{x})$ with implicitly valid Jacobi identity, based on the general solution of Jacobi identity in 3D (20). The loss function is the same as in the case of the WJ method, but this time Jacobi identity is valid automatically,
  $$\begin{align} \mathcal{L}_{\mathrm{IJ}} = \mathcal{L}_{\mathrm{WJ}}. \end{align} \tag{ 39c }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The training itself then proceeds in the following steps:

• 1.  
  Simulation of the training and validation data. For a randomly generated set of initial conditions, we simulate a set of trajectories by IMR. These trajectories are then split into steps (couples of two consecutive states) and the collection of steps is split into a training set and a validation set (40% used for the validation set). The collection of steps is shuffled and split to batches (20 steps per a batch).
• 2.  
  Parameters of the neural networks WJ, SJ, and IJ are trained by back-propagation on the training data (using the Adam optimizer), minimizing the loss function. The loss function differs in the three flavors of DPNNs, as in equation (39). Then, the loss function is evaluated on the validation data to report the errors. The learning rate was set to 0.001 while the number of epochs is chosen appropriately for each physical system.
• 3.  
  A new set of initial conditions is randomly chosen and new trajectories are generated using IMR, which gives the ground truth (GT).
• 4.  
  Trajectories with the GT initial conditions are simulated using the trained models for L and H and compared with GT.

In the following sections, we illustrate this procedure for learning RB mechanics, a particle in two dimensions (P2D), Shivamoggi equations, a P3D, and HT dynamics.

Figure 4(a) shows a sample trajectory of RB dynamics (25) from the GT set, as well as trajectories with the same initial conditions, generated the DPNNs. The training was carried out on 200 trajectories while GT consisted of 400 trajectories. Number of training epochs was 50. The initial conditions for the training and validation datasets were sampled uniformly from a ball of radius 11.18 (in the m variable), each trajectory consisted of 200 steps, and the moments of inertia were chosen as $I_x = 10.0$, $I_y = 20.0$, and $I_z = 40.0$.

Zoom In Zoom Out Reset image size

Errors of the three learning methods (WJ, SJ, and IJ) are shown in table 1. All three methods were capable to learn the dynamics well. Figure 4(b) shows the norm of the Jacobiator evaluated on the validation set. Jacobiator is zero for IJ and small in SJ, while it does not go to zero in WJ. Therefore, IJ and SJ satisfy Jacobi identity while WJ does not.

Table 1. Summary of the learning errors for a rigid body (RB), particle in two dimensions (P2D), Shivamoggi equations (Sh), particle in three dimensions (P3D), and heavy top (HT).

Mod.	Met.	$\Delta\mathbf{M}$	$\Delta\mathbf{r}$	$\Delta\mathbf{L}$	$\Delta\mathbf{M}^2$	$\Delta(\mathbf{r}\times\mathbf{M})$	$\Delta\mathbf{r}^2$	$\Delta(\mathbf{M}\cdot\mathbf{r})$	$\det \mathbf{L}$
RB	WJ	$1.5 \times 10^{-03}$	N/A	−5.0	$5.2 \times 10^{-04}$	N/A	N/A	N/A	0
	SJ	$4.6 \times 10^{-04}$	N/A	−6.5	$3.9 \times 10^{-04}$	N/A	N/A	N/A	0
	IJ	$1.7 \times 10^{-03}$	N/A	−7.4	$4.8 \times 10^{-04}$	N/A	N/A	N/A	0
P2D	WJ	$3.2 \times 10^{-03}$	$3.2 \times 10^{-03}$	−1.8	N/A	$9.7 \times 10^{-03}$	N/A	N/A	0.83
	SJ	$2.9 \times 10^{-03}$	$2.9 \times 10^{-03}$	−1.8	N/A	$8.6 \times 10^{-03}$	N/A	N/A	0.93
Sh	WJ	N/A	$3.4 \times 10^{-04}$	−2.3	N/A	N/A	N/A	N/A	−2.1
	SJ	N/A	$3.0 \times 10^{-04}$	−3.3	N/A	N/A	N/A	N/A	−3.4
P3D	WJ	$1.9 \times 10^{-02}$	$1.9 \times 10^{-02}$	−1.6	N/A	$3.7 \times 10^{-02}$	N/A	N/A	0.91
	SJ	$1.0 \times 10^{-02}$	$9.1 \times 10^{-03}$	−1.9	N/A	$3.5 \times 10^{-02}$	N/A	N/A	0.70
HT	WJ	$3.1 \times 10^{-02}$	$9.9 \times 10^{-03}$	N/A	N/A	N/A	$1.9 \times 10^{-03}$	$2.8 \times 10^{-03}$	−2.5
	SJ	$1.6 \times 10^{-02}$	$4.6 \times 10^{-03}$	N/A	N/A	N/A	$1.3 \times 10^{-03}$	$1.9 \times 10^{-03}$	−2.6

Figure 5 shows the compatibility error of learning the Poisson bivector (21). All three methods learn the Poisson bivector well, but IJ is the most precise, followed by SJ and WJ. Finally, figure 6 shows errors in learning the trajectories $\mathbf{M}(t)$. All three methods learn the trajectories well, but in this case, the SJ method works slightly better.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1 shows the medians of the errors for all the models and applied learning methods. N/A indicates quantities that are not to be conserved. Error $\Delta\mathbf{M}$ is calculated as the median of square deviation of M² over all time steps. Error $\Delta\mathbf{r}$, $\Delta \mathbf{M}\cdot\mathbf{r}$, $\Delta \mathbf{M}^2$, and $\Delta \mathbf{r}^2$ are calculated analogically. Error $\Delta\mathbf{L}$ in the RB case is calculated as $\log_{10}$ of the L² norm of the compatibility condition (21), calculated for the learned J divided by its trace and multiplied by factor 1000 (using the exact J when generating the GT). In the P2D and P3D cases, where the Poisson bivector is symplectic, the error is calculated as log₁₀ of squares of the symplecticity condition (16). In the case of Shivamoggi equations, the $\Delta \mathbf{L}$ error is the log₁₀ of the squared superintegrable compatibility condition (31). $\Delta \det\mathbf{L}$ errors are medians of squares of learned $\det\mathbf{L}$, and in the Shivamoggi and the HT cases, the values are logarithmic, since determinants are supposed to be zero in those cases.

A particle moving in a 2D potential field represents a 4D symplectic system. Initial conditions were sampled uniformly from balls of radius 11.18 in the momentum and 10.48 in the positions (200 points), and each trajectory had 200 steps. The Hamiltonian was chosen as $H(\mathbf{r},\mathbf{m}) = \mathbf{m}^2 + \mathbf{r}^2$. GT was generated in the same way, but with 400 sampling points.

The simulated trajectories were learned by WJ and SJ methods (each with 50 epochs). No implicit IJ method was used because no general solution of Jacobi identity in 4D is available that would work for both the degenerate and symplectic Poisson bivectors. Results of the learning are in table 1, and both WJ and SJ learn the dynamics comparably well. Figure 7 shows a sample trajectory, and figure 8 shows the distribution of learned $\det(\mathbf{L})$. The median determinant (after a normalization such that the determinant is equal to 1.0 in GT), was close to this value for both SJ and WJ, indicating a symplectic system.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Shivamoggi equations (32) represent a 4D Hamiltonian system that is not symplectic, and thus has a degenerate Poisson bivector. The equations were solved within the range of parameters $u\in [-0.5, 0.5]$, $x\in [-0.5,0,5]$, $y\in [-0.1, 0.1]$, $z\in [-0.5,0.5]$. It was necessary to constraint the range of $\mathbf{r} = (u,x,y,z)$ because for instance when u = 0, the solutions explode [41]. Figure 9 shows the u-component over a sample trajectory. The initial conditions were sampled uniformly (500 points) as well as the GT (also 500 points), and each trajectory consisted of 150 steps.

Zoom In Zoom Out Reset image size

We used 80 epochs in the learning, and results of the WJ and SJ methods are in figures 9 and 10. Figure 10 shows the distribution of $\log_{10}(\det(\mathbf{L}))$, indicating that the system is indeed degenerate. In comparison with determinants of L learned in the symplectic case of a 2D particle (P2D), see table 1, the learned determinants are quite low in the Shivamoggi case (after the same normalization as in the P2D case). Therefore, DPNNs are able to distinguish between symplectic and non-symplectic Hamiltonian systems.

Zoom In Zoom Out Reset image size

Motion of a P3D was simulated and learned with the same parameters as the P2D (section 3.4). Figure 11 shows momentum during a sample trajectory of a P3D space taken from the GT set, as well as trajectories with the same initial conditions obtained by DPNNs (WJ and SJ flavors). Training (and validation) dataset and GT dataset consisted of 200 and 400 trajectories, respectively. Table 1 contains numerical values of the learning errors. Both WJ and SJ learn the Poisson bivector as well as the trajectories (and thus also the energy) well. The median determinant is close to unity, which indicates a symplectic system.

Zoom In Zoom Out Reset image size

The HT dynamics is 6D, as well as the P3D. Parameters of the HT were the same as those of RB (section 3.3), except for taking 300 sampling points for the training and validation data, both $c_{\mathrm{mov}}$ and $c_{\mathrm{Jac}}$ were set to 100, number of epochs was 100, and constant Mgl was set to 0.981.

Figures 12(a) and (b) show a sample trajectory of a HT from the GT set and trajectories with the same initial conditions obtained by DPNNs. The training and validation were done in two sets of trajectories (with 300 and 400 trajectories, respectively). Numerical values of the learning errors can be found in table 1. For instance, the L matrix is close to being singular, indicating a non-symplectic system, but SJ learns slightly better than WJ. Similarly, as in the four-dimensional case, DPNNs distinguish between symplectic (P3D) and non-symplectic cases (HT).

Zoom In Zoom Out Reset image size