For a given differential equation, a numerical method is called a geometric numerical integrator if it can accurately preserve some of the geometric characteristics of its solution [HLW06]. The structure-preserving algorithm of differential equations has made a lot of important development, and the basic idea of maintaining the important structure of original differential equation for numerical methods is widely accepted. For the Hamiltonian system, the symplectic integrator for preserving the symplectic geometry property of the solution has been proved to have very good orbital tracking ability for a long time and various excellent properties [FEN86, FQ03, HLW06, MW01]. Many effective methods for constructing symplectic integrators, such as the methods based on variational integrators [MW01], generating functions [FQ03] and Rung-Kutta (RK) methods [SUN93, SUN95], have been developed and investigated. With the establishment of backward error analysis [3, 27] and discrete KAM theory 
, the symplectic integrators for Hamiltonian system is becoming more and more perfect. Many new related numerical methods such as multi-symplectic methods for Hamiltonian partial differential equations[QW11] and stochastic symplectic methods for random Hamiltonian systems [CHJ20] have been well developed. In addition to the symplectic property, another extremely important feature of the Hamiltonian system is the conservation of energy. Therefore, it is very important to keep the system energy by numerical method. The important energy preserving methods include discrete gradient method [GON96, 20]
, average vector field method, HBVMs [BIT09] and spectral methods [JWZ19]
When multiple important structures or physical quantities exist for a system, such as the symplectic structure and the energy for Hamiltonian system, a natural question arises: is it possible to construct a numerical method that preserves several of them? Regarding the symplectic structure and the energy for Hamiltonian system, one has unfortunately a negative answer in general for constant step size. In fact, it is proved  that for non-ingregrable systems if the method is symplectic and can conserve the Hamiltonian energy exactly, then it is the time advance map for the exact Hamiltonian systems up to a reparametrization of time. A similar negative results is proved in  for general Hamiltonian system by B-series method. But the above negative results do not prevent people from constructing numerical methods to maintain both the energy and symplectic structure of the Hamilton system in some weaker sense.
A breakthrough work in this regard is the parameterized Gauss collocation method firstly developed by Brugnano in . Consider the canonical Hamiltonian systems in the form
is the identity matrix andis the Hamiltonian energy. Brugnano et.al. introduced a nice idea to develop a new family of Gauss type methods which share both the symplecticity-like and energy conservation features under suitable conditions. More precisely, they define a family of RK methods , where is the step size of integration, is a real parameter. This method satisfies the following three conditions simultaneously: (i) for one gets the Gauss collocation method of order , is the number of stages of the RK method; (ii) for any fixed choice of , the corresponding method is of order and satisfies the conditions , thus being a quadratic invariant-preserving RK method; (iii) for any choice of and in a neighborhood of , there exists a value of parameter such that (energy conservation). The resulting method has order , preserves the energy and quadratic invariants 111For stand RK method with constant step size , the sufficient condition of symplecticity is given by for all . For irreducible RK method, which can be intuitively understood to mean that this RK method is not equivalent to a RK method with a lower series, see more details in [HW06, pp.187], this condition is also necessary, and also sufficient and necessary for quadratic invariants-preserving. However, for the parametric -RK methods that preserve the energy, the parameter depends on the initial values, i.e., . Then, -RK methods with condition are in general not symplectic from the definition, but can preserve all the quadratic invariants..
This method has been rewritten in the framework of discrete linear integral methods , which leads to a more refined theoretical analysis and a nice practical implementation strategy for seeking the parameter . Another extension form -RK to -PRK to make use of the additive structure of Hamiltonian system is given in . The key of constructing parameterized RK or PRK methods in [6, 28] lies in the so-called -transform technique. However, this technique is in general not suitable for constrained Hamiltonian systems since it will broke the constrained manifolds. See more details in Section 4.
In this paper, we focus on the geometric integrators for Hamiltonian systems subject to holonomic constraints. Consider the Hamiltonian with holonomic constraints (constraints that depend on only) [HLW06, Chap. VII]
where we introduced to be the Jacobian matrix of so that Note that we use the convention which is commonly used in the community of fluid mechanics. Throughout this paper, we assume that the matrix has full rank and is invertible, which will allow us to express in terms of . The dynamics given by (1.3) is on the manifold
which is the cotangent bundle of the manifold given by See [HLW06, VII.1.2] for more details. The constraint is also called the hidden constraint.
It is easily to verify that the ODE (1.3) has many geometric structures. First of all, the Hamiltonian (or energy) is preserved:
Secondly, the symplecticity is also preserved. Consider a map and its tangent map . One says the map is symplectic if
For being a submanifold of , so that
for any smooth curve in satisfying . It can be shown that the flow map of (1.3) is symplectic. Plainly speaking, the symplecticity means that the flow preserves the quadratic forms, which further implies that the volume form is conserved (see the introduction of [JAY96]).
From the pointview of structure-preserving methods, if one solves the Hamiltonian system (1.3) numerically, besides the basic requirement that the solution falls onto the manifold , one also tends to ask that the energy preservation or the symplecticity is satisfied. This additional requirement for staying on manifolds for constrained Hamiltonian systems presents new difficulties in the construction of structure preserving numerical schemes.
There are many symplectic methods for the holonomic constrained Hamiltonian system. One typical class is the PRK methods as detailed explained in [JAY96, MW01, JAY02]. The first order method in [HLW06, sec VII.1.3] and the Shake-Rattle method in [AND83, HLW06] are special cases of the PRK methods. Symplectic variational integrators for constrained Hamiltonian systems are developed in [MW01, WOL17]. Those methods are symplectic and preserve the manifolds exactly, but they are not Hamiltonian energy converved. The second order energy-conserving method in the line integral framework developed in  is not symplectic in general and can not conserve the the hidden constraint . The projected RK methods are proposed in , but it is not symplectic.
The goal of this work is to introduce a parameter in the PRK methods so that both the energy and the quadratic invariants can be preserved for the constrained Hamiltonian system following the spirit of [6, 28]. Note that if we do this method for every trajectory, then the energy can be preserved for all trajectories and the method is in general not symplectic. However, if we fix down the parameters, the method is symplectic and can be preserve energy for one trajectory.
The rest of the paper is arranged as follows. In section 2, we recall some of the basic concepts and properties of symplectic PRK methods for constrained Hamiltonian systems. In order to show the symplectic PRK methods preserve the quadratic invariants for Hamiltonian systems subject to holonomic constraints, we provide a new variational formulation for symplectic PRK methods in section 3. In section 4, we construct parametrization -PRK methods based on Shake-Rattle method and Lobatto IIIA-IIIB pairs, and check various properties of the new schemes. Numerical examples are included in section 5 to illustrate the energy and quadratic invariants conservations of the -PRK methods.
2 Symplectic Partitioned Runge-Kutta methods
In this subsection, we give a brief review of the symplectic PRK methods. See [JAY96, HLW06] for more details. The -stage PRK methods denoted by for the Hamiltonian ODE (1.3) is given by the following
[JAY96] For a given a PRK method, if the conditions for symplecticity
is satisfied, then the method for Hamiltonian systems with holonomic constraints is also symplectic.
For -stage PRK methods, if are given, then there are equations for the unknowns (). However, for general symplectic PRK methods, may not fall onto . To ensure that satisfies the constraint, it is suggested in [JAY96] to impose the conditions
With this, one has correspondingly
from the symplectic condition given in Lemma 2.1. This requirement removes the unknowns and the equations , . Hence, the number of equations is now equal to the number of unknowns . In fact, it is shown in [JAY96] that these variables can be uniquely solved when is small enough under some reasonable assumptions on and . We remark that even though one imposes , the Lagrange multiplier is still in the system of equations. Lastly, and can be determined by the hidden constraint condition
For the PRK methods suggested in [JAY96], we hence solve for the next two nonlinear systems in turn at each step. The first one
with Then, we can solve by
As commented in [HLW06, sec VI.6.2], the symplectic schemes can be written as certain variational integrators. In [HLW06, sec VI.6.3], it was explained in detail how the PRK schemes for unconstrained problems can be formulated into a variational integrator. The same is true for Hamiltonian systems with holonomic constraints. In fact, Marsden and West addressed this issue in [MW01, section 3.5.6]. In particular, given , one can implicitly define quantities for and for through the following system of equations
Note that the definitions of and above are different from other (in fact, for while there should be some corrections for and to define ). Moreover, and are not in (2.1). They are different by some terms involving the corresponding Lagrange multipliers. See the proof of [MW01, Theorem 3.5.1] for the details. With the quantities defined in (2.8), one can define the discrete Lagrangian
where is the Lagrangian for the time-continuous dynamics. The Euler-Lagrange equation for this discrete Lagrangian under the holonomic constraints reads
See [MW01, eqs. (3.5.2a-d)]. Note that though the Euler-Lagrange equations are derived by fixing , one should regard as knowns and solve dynamically. The good thing of this variational formulation is that many conservation properties can be derived via the discrete Noether theorem. See [MW01] and also section 3 below for more details on conservation properties and discrete Noether theorem. However, for quadratic invariants, the form (2.8)-(2.9) is not convenient to use and we will propose another variational formulation below in section 3.
2.1 PRK schemes: Shake-Rattle and Lobatto IIIA-IIIB pairs
We first recall some famous PRK methods, including Shake-Rattle algorithm and Lobatto IIIA-IIIB pairs, which will be used to construct the parameterization -PRK later. For typical separable Hamiltonian
the traditional Shake method reads
Then, the momentum is The corresponding Hamiltonian formulation for the Shake method is
Unfortunately, so-defined may not satisfy the constraint. Anderson [AND83] proposed the Rattle algorithm. That is to determine by using another multiplier which can be different from :
Then, one uses the condition
to determine and thus . Clearly, this Shake-Rattle algorithm is one of the symplectic PRK method mentioned above. The more general form of Shake-Rattle for the general Hamiltonian system with holonomic constraints reads
It is proved that the Shake-Rattle algorithm is symmetric, symplectic and convergent of order two [HLW06, sec VII.I].
The Shake-Rattle algorithm is second order. A nice extension to higher order is the Lobatto IIIA-IIIB pairs developed in [JAY96]. This method combined with the projection step by the hidden constraint condition defined in the manifold has the following merits. (a) preserving the numerical solutions on the manifold exactly; (b) it is symmetric, symplectic and super convergent of order , the is the stage of PRK method. In particular, for , the coefficients for Lobatto IIIA-IIIB pairs are given by
3 An alternative variational formulation for symplectic PRK schemes with holonomic constraints
In this section, we aim to study the quadratic invariants of the PRK schemes via the discrete Noether’s theorem. For this purpose, we need to find some equivalent variational forms for the PRK scheme. However, the form in the work of Marsden and West (i.e. (2.8)-(2.9)) is not convenient, as it involves Lagrangian multipliers and there is nonlinearity. Instead, we will propose an alternative formulation so that the symmetry can be studied better.
It is known that symplectic PRK schemes for Hamiltonian systems without constraints can conserve quadratic invariants of the form
where is some fixed matrix ([HLW06, sec VI.2.2]). Naturally, one is curious whether there is an analogue for Hamiltonian systems with holonomic constraints. Clearly, due to the constraint , the quadratic invaraints must be of some particular forms. Instead, motivated by Noether’s theorem [NOE71, ARN13], we will consider that
Clearly, by Noether’s theorem, we have the following.
For any , the quantity is a first integral of the Hamiltonian system with holonomic constraints.
We now aim to show that the symplectic PRK schemes will conserve all the quadratic invaraints of the form for .
The symplectic PRK schemes for Hamiltonian system with holonomic constraints defined in (2.1) conserve all quadratic invariants of the form
To do this, we first recall the discrete Noether’s theorem for variational integrators without constraints.
Lemma 3.2 (Discrete Noether’s theorem).
Assume the discrete Lagrangian is invariant under a one-parameter group of transformation : for all and . Then the corresponding variational integrators for Lagrangian systems have the first integral in the form of
where is the generator for the group. Furthermore, when there are holonomic constraints, besides the conditions above, if moreover leaves the constraints manifold invariant, then the claim still holds.
Proof of the second part of Lemma 3.2.
Taking derivative on and setting in , one has
Since leaves the constraints manifold invaraint,
Using (2.10), one easily finds . ∎
Below, we propose an alternative variational formulation for the PRK scheme that appears different from (2.8)-(2.9). This formulation is suitable to verify the conditions for the discrete Noether’s theorem.
Note that we are not defining as , so that , . In this definition, do not have to be in .
Proof of Proposition 3.1.
Step–1: First of all, we note that the discrete Lagrangian is well-defined. The conditions for the extremizers are given by
Note that we have built in the constraint by using the multiplier . Since , this will not change anything. These equations together with ( equations) and the constraint ( equations) can yield the solutions for sufficiently small.
Here we are solving in terms of for any in , not just for . With this observation, we can take differentiation on for the constraint (3.6) to find
and for the constraints (3.5) to find
Note that here we used the convention
The reason to do this is that we have made a column vector.
Defining that one has
Step–3: We verify that the Euler-Lagrange equations (2.10) given by so-defined discrete Lagrange is equivalent to the symplectic PRK methods with constraints. Recall that the (continuous) Hamiltonian is related to Lagrangian by where is determined by the second equation. It can be computed that Hence, we have
where The first equation in (3.7) can then be rewritten as
Using (3.12), we then have
Using the condition , this is further simplified to
Defining one has the desired form for .
Lastly, we verify Theorem 3.1.
Proof of Theorem 3.1.
The constraint condition is obvious. We verify that