Symmetry Preserving Discretization of the Hamiltonian Systems with Holonomic Constraints

: Symmetry preserving difference schemes approximating equations of Hamiltonian systems are presented in this paper. For holonomic systems in the Hamiltonian framework, the symmetrical operators are obtained by solving the determining equations of Lie symmetry with the Maple procedure. The difference type of symmetry preserving invariants are constructed based on the three points of the lattice and the characteristic equations. The difference scheme is constructed by using these discrete invariants. An example is presented to illustrate the applications of the results. The solutions of the invariant numerical schemes are compared to the noninvariant ones, the standard and the exact solutions.


Introduction
All variables are measured in a certain time interval, and all simulations can be implemented in meshes. Difference equations serve as primary and fundamental mathematical models in physics and mechanics. Symmetries concern the basic principles of physics and provide a most valuable heuristic guide in the search for dynamical laws. Continuous symmetries should be incorporated in the discretization of continuous systems, but they are usually lost in discrete descriptions [1,2].
Lie group theory is an efficient tool for solutions of differential equations [3]. It is also more convenient for the discretization of both ordinary differential equations (ODEs) [4][5][6] and partial differential equations (PDEs) [7][8][9]. For a given finite-difference equation (system), one can have the operators of the group admitted by the equation. Different types of difference equations which admit the different symmetry transformation operators can be obtained [10,11], and this theory is applied to search for discrete symmetries for dynamical systems [12,13]. Generally speaking, it needs to solve some linear partial finite-difference equations. But the problem of exact integration of partial finite-difference equations is always intractable, even when the equations are linear. Therefore, in order to construct the symmetry preserving difference schemes, the inverse problem needs to be solved, i.e., the difference equations and meshes are be constructed when the symmetries (transformation groups) are known. It is also a simpler way of having the symmetry preserving discretization of the differential equations.
It is better to keep the geometry structure during the discrete process. Payen, Matignon and Haine [14] addressed the discretization of the boundary controlled 3D Maxwell's equations as a port-Hamiltonian system, and the proposed scheme preserves the Dirac structure at the discrete level. The symmetry structure is the geometry structure of the systems. It is meaningful to maintain the symmetry of the system in the discrete process. Verstappen and Veldman [15] studied the symmetry-preserving discretization of turbulent flow. Bourlioux, Rebelo and Winternitz [16] discretized the nonlinear ODEs under the group SL(2, R) using the symmetry preserving discretization method. The results showed that solution methods incorporating the Lie point symmetries provide better results than standard methods. The symmetry preserving schemes not only serve as fundamental mathematical models but also preserve main geometric properties of the dynamical systems [17,18]. Dorodnitsyn et al. proposed approaches on the applications of Lie group theory to difference equations and discrete mechanical systems [19]. This approach is to start with a differential equation and to introduce a symmetry adapted mesh and difference equation in such a way that all the symmetries of the original differential equation are preserved [20]. Levi's [21], Winternitz's [22], and Rebelo's [23] contributions are the typical representatives for this approach.
There have been some papers on constructing the invariant discretizations of dynamical systems based on the inverse problem method above, such as the one-dimensional gas systems [24], the shallow water equations in Lagrangian [25], the holonomic dynamical systems [23], etc. The symmetry preserving discretization of ODE and PDE are studied deeply [9,10,[19][20][21][22][23]. The theory has been applied to dynamical systems recently. The references [7,12,18,[24][25][26] studied symmetry preserving schemes in a Lagrangian framework. In the Lagragian framework, the symmetry-preserving structure is studied in the space M = Q × R or M = TQ × R. The transformation group ψ : M → M . The Lie transformation generators are τ α , ξ µ α , η µ α : M → R . In this paper, the symmetry-preserving structure is in the dual space M = T * Q × R of the tangent space. The conditions of preserving the symmetries of the Hamiltonian systems are studied, and the corresponding numerical simulation are proposed. The goal of this article is to extend the method of invariant discretization of differential equations to the holonomic Hamiltonian systems. This new method can preserve fundamental symmetric properties and improve the features of numerical algorithms in a Hamiltonian framework. The main goal of this paper is to find a way to discretize the constrained systems in the Hamiltonian framework. Generally speaking, the constrained dynamical systems are universal in engineering but relatively difficult in terms of getting a solution. Suppose the systems are subject to holonomic constraint. Although it is simpler, the results could make some preliminary analysis in the numerical simulation of complex systems.
The paper is structured as follows. In Section 2, the Lie symmetry of Hamiltonian systems with holonomic constraints are recalled, based on the infinitesimal generators and the prolongations. In Section 3, the Lie point symmetry group in Section 2 is prolonged to the three points of lattice, and the invariant difference scheme is obtained by solving the discrete characteristic equations. The symmetry preserving discretization of the Hamiltonian systems with holonomic constraints are constructed by the difference invariants. Section 4 is devoted to an example to illustrate the theoretical results, and the Hamiltonian system's behaviors are simulated using the invariant discretization. Some conclusions are drawn in Section 5. Throughout the paper Einstein's summation is used.

The Equations of the Constrained Hamiltonian System
Let X = R m be the space representing the independent variables, and U = R n represent the dependent variables. Consider a mechanical system whose configuration is determined by n generalised coordinates q s (s = 1, . . . , n). Let the motion of the system be subjected to the g holonomic constraints Q s = Q s t, q s , . q s , the differential equations of motion of the system can be written in the form where L = L(t, q s , . q s ) is the Lagrangian, Q s are the non-potential generalised forces. Introduce the generalized momentum and Hamilton function and the canonical form of equations of the Hamiltonian systems with holonomic constraints are involving one independent variables t and 2n dependent variables p s = (p 1 , . . . , p n ) and q s = (q 1 , . . . , q n ) where ∂ q s = ∂/∂q s , ∂ p s = ∂/∂p s , H is the Hamiltonian, Q s are the nonpotential generalized forces as functions of t, q s and p s . Expanding Equation (3), we can solve . q s and . p s as . q s = g s (t, p s , q s ),

Infnitesimal Generators and the Prolongations
A symmetry group of the system will be a local group of transformation G r in Z acting on the open subset M ⊂ X × U. For the holonomic Hamiltonian systems with X = R and U = R 2n , we denote a local one parameter continuous transformation group G 1 concerning t, q s and p s to a continuous parameter ε ∈ R 2n can be taken from Accordingly, the transformations in the group G 1 acting on the variables can be written as since ε→0 forms the identity of the group. The infinitesimals are defined by the new functions: Then, Equation (5) becomes The first prolongation of the infinitesimal generators should also be known for the Hamiltonian systems with holonomic constraints. By obtaining . q * s and . p * s : where D t is the total derivative operator for time. Ignoring the high-order terms, we have One can determine the tangent vector field γ = {τ, ξ s , η s } T , of the group at the point (t, q s , p s ), where τ, ξ s and η s are unknown and demand to be determined to achieve the symmetry transformation. The Lie algebra L of the symmetry group G 1 is realized by vector fields of the form X (0) and the first prolongation is where ∂ t = ∂/∂t. In the same way, the vector fields and the prolongation of systems (1) in Lagrangian framework are X (0)

Lie Symmetries for Hamiltonian Systems
The algorithm for finding the symmetry algebra and the symmetry group for a given Equation (3) goes back to S. Lie and is given in many books on the subject [3]. For the equations of Hamiltonian systems, the symmetry algebra L and the symmetry group G 1 are given as follows.

Definition 1.
If G 1 is a local group of transformations acting on M, and for every infinitesimal generator of G 1 , then G 1 is a Lie symmetry group of the systems.
Considering Equation (4), the determining equations of the Lie symmetries for Equation (3) under the transformation (8) have the form where g s = ∂ p s H and h s = −∂ q s H + Q s . This produces 2n partial differential equations to determine the functions τ, ξ s and η s. These equations provide an effective computational procedure for finding the most general symmetry group of the systems. In this procedure, the coefficients τ, ξ s and η s of the infinitesimal transformation generator X (0) s of oneparameter symmetry group of the systems be unknown functions of t, q s and p s . The determining Equations (14) will thus involve t, q s and p s , as well as τ, ξ s , η s and their partial derivatives with respect to t, q s and p s . We can equate the coefficients of the remaining derivative of q s and p s to zero after eliminating any dependencies among the derivative of the q s 's and p s 's caused by the system itself. This will result in solving a large number for the coefficient functions τ, ξ s and η s for the infinitesimal generator. In most instances, these determining equations can be solved by elementary method, and the general solution will determine the most general infinitesimal symmetry of the system.
We can also use an existing software program designed to automate the major steps of this algorithm in order to calculate the symmetry group of the system more efficiently, such as Mathematica [27], Maple [28], etc. In this paper, the symmetry operators are obtained in Maple. In general, the Low-order differential equations can be solved by the determine() provided by liesymm, an embedded software package in Maple system [29]. The solution of symmetric equations of higher order nonlinear partial differential equations can be obtained by means of some standard software packages like PDEtools [30], SADE [31], GeM [32], etc.

Invariance of Difference Equations of Hamiltonian Systems
Consider the space Z of sequences (t, q s , p s ), and denote by D the first order linear discrete operator Fixing arbitrary parameter value h > 0, we form a pair of the operators of discrete transformation to the right and left where D is a derivation in Z. The operators S +h and S −h commute with each other and satisfy Consider Hamiltonian difference equations at some lattice points (t, q s and p s ). Generally, the lattice is irregular.
Using S +h and S −h we define a pair of right and left discrete (finite-difference) differentiation operators by setting where D can be written as where This equation is writ ten on finitely many points of the difference mesh ω h . We assume that the mesh is determined by the equation where Ω ∈ A h . The function Ω is uniquely determined by the discretization of the space of independent variables. In the continuous limit one equation, the Equation (16) reduces to an identity (like 0 = 0). The vector fields are the same as in the continuous case, however they must be prolonged to all points of the lattice, involved in the system (15) and (16). We have G 1 be a one-parameter group in Z with operator then the symmetry is the Lie symmetry of discrete Hamiltonian systems with holonomic constraint.
The mesh spacing transformations are h + and h − . By using the uniform mesh h + = h − and applying the infinitesimal transformation (7), we obtain prX The meshes satisfying (20) are said to be invariantly uniform for the systems.

The Lie Symmetry-Preserving Difference Scheme for Hamiltonian Systems
Suppose the Hamiltonian Equation (3) admit a known transformation group G 1 without loss of generality, after the complete set of r functionally independent invariants (I C 1 (z), . . . , I C r (z)) of order k is constructed, system(3) can be rewritten in the invariant representation Ψ a (I C 1 (z), . . . , I C r (z)) = 0. For the systems of the difference Equation (15) and the difference mesh (16) which admits the same group G 1 and provide a first order approximation to the original systems (3). We wish to construct the difference scheme that not only approximates Equation (3), but has the same Lie point symmetry group G 1 . This is achieved by constructing the scheme out of difference invariants of the group G 1 , or out of invariant manifolds. These are found using the vector fields (8) and (9) corresponding to the invariance algebra of Equation (3). We can proceed as follows: (i) The Lie symmetry groups of the original Equation (3) of Hamiltonian systems are obtained. Then the corresponding vectors fields are prolonged to three point on a line and the corresponding values qi = qi (t), q _ i = qi (t _ ), q + I = qi (t+), p _ i = pi (t _ ), p _ i = pi (t+).
(ii) The set of difference invariants can be found by solving the following standard linear problem of group analysis: the local invariant I of G 1 is the solution of Equation (21). The method of characteristics gives us a set of independent elementary invariants of the corresponding Lie group action. The invariance conditions (18) and (19) yield r + nm functionally independent difference invariants (I D 1 (z), . . . , I D r+nm (z)) because there are two sets (right and left) of first difference derivatives in the difference space.
(iii) The invariant mesh from the general equation for invariant mesh generation is where the ω β are arbitrary smooth functions. The choice can be made from different standpoints. Any other invariant, i.e., any other solution of(16), will be necessarily be a function of I D 1 (z), . . . , I D r+mn (z). (IV) Form a set of r invariants with the desired approximation property; i.e., for each invariant I D α represen table in Z by the Taylor group, the relation is which holds on the mesh chosen above. In practice, as a rule, it suffices to have less τ of such invariants. The Equation (15) can be written in the invariant representation From the above statement, we have the following proposition: (22) and (24) satisfy: (i) Equation (22) reduce to zero, and Equation (24) reduce to Equation (3) when h → 0 ;

Proposition 2. For the Hamiltonian system (3), if the vector field of the Equations
(ii) prX ω β (I D 1 (z), . . . , I D r+mn (z)) Ψ D γ (I D 1 (z), . . . , I D r (z)) = 0. Then the Equations (22) and (24) are symmetry preserving discretization of the Hamiltonian system. (22) and (24), if they can reduce to the continuous Equation (3) when h → 0 , they are the approximation of the continuous systems. Equations (22) and (24) are Lie symmetry if they satisfy Equation (40). They have the same Lie group as the continuous system. Equations (22) and (24) can preserve the symmetry of the Hamiltonian system.

Proof. For Equations
Since the functions Ψ D γ = 0 are assumed to be locally analytic in their arguments, it follows from(21) that the difference equations Ψ D γ = 0 model the corresponding differential equations. The system of difference Equations (22) and (24) constructed admits the complete group G 1 , and they are the symmetry-preserving discretization of the Hamiltonian systems with holonomic constraints.

Examples
Example 1. Consider a nonlinear holonomic system with Lagrangian L = q 2 /2 with the nonpotential generalized forces Q = exp q. The Euler-Lagrangian function is Applying the formula Equation (12) in Equation (25), we obtain By separation of the coefficient of the above system (26), the determining equations of the holonomic system are Then the general solutions of the Equations (27)-(30) are The symmetrical generators are given as The second prolongation of (32) is The equations of holonomic system (25) in Hamiltonian framework is which admits the symmetry group associated with the Lie algebra spanned by the operators The symmetries of the holonomic system in the Lagrangian framework can be calculated in Maple.
We have at least a three-point difference stencil, on which we should construct (to be define) a first-order approximation to the original Equation (34). Let us represent the operators (35) and (36) by extending to the space (t, q, p, q t , p t , h + , h -) as follows (the corresponding coordinates in X 1 are zero): It was shown that a transformation defined by (7) conserves uniformity of a mesh in time if and only if D +h D −h (τ) = 0, so the mesh is uniform. In the following, the independent invariants which can be obtained by solving the standard linear problem of group analysis (38): We should write out two invariant difference equations as the form The simplest version I D are taken. This mesh has the obvious integral h + = h -= ε, ε = const; 0 < ε 1, where the constant ε characterize the mesh spacing smallness. Using Taylor series expansions, we can obtain the relation h + = h -+ O(h 2 ) for the spacing of the invariant mesh (41). On such a uniform mesh, the equation ω(I D α ) = 0 for invariant mesh generation provides infinitely many possibilities. We choose Thus, for the difference equation approximating the differential equation of the systems we can take The constructed invariant model is not unique. In the continuous limit, operators (37) and (38) reduces to (35) and (36) respectively, so the discrete invariants will reduce to the continuous differential equations of the Hamiltonian system when h → 0 . Equations (45) and (46) can reduce to (34), but (47) can not although it can be written in invariant form. So Equation (47) is not the discretization of the systems. Equations (45) and (46) give two approximations to the differential invariant . p s /exp q model O(h 2 ). But, even on approximation to (34), some difference ones are still not the symmetry-preserved discretization. For example, the equation has approximation to Equation (34) as well but can not be written in invariant form. Equation (48) can not be constructed in the form of I D 1 (z), . . . , I D 5 (z), it can not be written as Equation (24). It is not the symmetry preserving scheme from Proposition 2. Equation (48) is not the symmetrypreserving discretization of the systems. Of course, the equation for invariant mesh generation provides many possibilities. In this case, we choose two symmetry preserving schemes according to Equations (45) and (46) p = q t , p t = exp q. (49) The constructed invariant models (49) and (50) admit X 1 and X 2 . It is obvious that Equation (48) admits X 1 but does not admit X 2 , it approximates(34) but does not admit the same group as the original equation. So Equations (49) and (50) are both the symmetry-preserving discretization of the systems, but Equation (48) is not.
The solution obtained with the three schemes (48), (49) and (50) is displayed in Figure 1 for the coarse resolution h = 0.01. The initial conditions are taken from the exact solution: q 0 = 0.1; p 0 = 0.105; h = 0.01. The solution is also computed by the Matlab standard and adaptive Runge-Kutta scheme ODE45 in Figure 1. They are compared with the exact solution of the Hamiltonian system. From Figure 1, the symmetry preserving schemes (49) and (50) provide approximations of the differential Equation (34). The accuracy and stability of symmetry preserving schemes is better than the standard scheme. In Figure 2, the solution is shown with the invariant scheme for three resolutions h = 0.1, h = 0.01, h = 0.001, the graphics show that smaller steps provide a better approximation. Table 1 compares the values at several points using the standard scheme, the symmetry preserving schemes (sym. pres. I and II) and the analytical solutions (exact). The total elapsed time for standard schemes and the symmetry preserving ones is listed on in Table 2. These two tables have shown that the accuracy of the symmetry preserving schemes is better than that of standard schemes at no significant additional cost.   Example 2. Consider a particle in the one-dimensional Coulomb field. The Lagrangian is L = q 2 /2. It is subjected to repulsive forces Q = 1/q 2 which is the holonomic constraint. The differential equation of motion of the system is In Hamiltonian framework, the Equation (51) can be written as Using Maple procedure, the Lie symmetry vector field for Equation (52) are The operators (53) and (54) can be extended to the space (t, q, p, q t , p t , h + , h -) as follows From Proposition 1 and the generators in vector fields (55) and (56), we have D +h D −h (τ) = 0, the mesh defined by (7) is uniform. Applying the method of characteristics to Equation (56), the independent invariants can be obtained We write out two invariant difference equations as the form The simplest version (58) of difference equation is just the map scheme. The Equation (59) can be expressed as which satisfy Proposition 1, it is symmetry-preserving discretization of the systems. There are different types of discretization for the system (52), such as Although it has approximation to (52), it can not satisfy (18) expressed asprX 2 p t − q (−2) − h 2 q = 0. So expression (61) is not the symmetry-preserving scheme.
The system (51) is a conservative mechanical system. The energy is constant. The energy with the schemes (60) and (61) is displayed in Figure 3 for the step h = 0.01. The initial conditions are q 0 = 10 and p 0 = 10. The results show that the symmetry preserving schemes can keep the energy constant. But the noninvariant scheme can not simulate the energy correctly.

Conclusions
The basic motivation for this research program is applying the method of symmetry preserving discretization to the holonomic Hamiltonian systems. The essential feature (symmetry) of the Hamiltonian systems is incorporated in the symmetry preserving difference mathematical model. The main theoretical results are the following: In the Hamiltonian framework, the Lie symmetry generators can be obtained from the determining Equation (14). The generators generate the Lie group and form the vector field. The discrete invariants can be obtained using the prolonged Lie symmetry vector field. The difference schemes (equations and meshes) are constructed out of the invariants of the corresponding Lie groups. The numerical examples show that the invariant schemes can simulate the dynamical behaviors accurately. Although the noninvariant schemes can give the original equations under continuous limit, they can not reflect the solution correctly.
It is significant to preserve the geometrical structure when the continuous systems are discretized. The invariant schemes can preserve the symmetry structure of the holonomic Hamiltonian systems. It proposes a preliminary analysis of the structure-preserving discretization of the constrained systems in the Hamiltonian framework.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.