# Energy and quadratic invariants preserving methods for Hamiltonian systems
with holonomic constraints

Lei Li E-mail<EMAIL_ADDRESS>School of Mathematical Sciences,
Institute of Natural Sciences, MOE-LSC, Shanghai Jiao Tong University,
Shanghai, 200240, P. R. China.  Dongling Wang E-mail<EMAIL_ADDRESS>Department of Mathematics and Center for Nonlinear Studies, Northwest
University, Xi’an, Shaanxi, 710127, P. R. China.

###### Abstract

We introduce a new class of parametricization structure-preserving partitioned
Runge-Kutta ($\alpha$-PRK) methods for Hamiltonian systems with holonomic
constraints. When the scalar parameter $\alpha=0$, the methods are reduced to
the usual symplectic PRK methods like Shake-Rattle method or PRK schemes based
on Lobatto IIIA-IIIB pairs, which can preserve all the quadratic invariants
and the constraints. When $\alpha\neq 0$, the methods are also shown to
preserve all the quadratic invariants and the constraints manifold exactly. At
the same time, for any given consistent initial values $(p_{0},q_{0})$ and
small step size $h>0$, it is proved that there exists
$\alpha^{*}=\alpha(h,p_{0},q_{0})$ such that the Hamiltonian energy can also
be exactly preserved at each step. We provide a new variational formulation
for symplectic PRK schemes and use it to prove that the parametrized PRK
methods can preserve the quadratic invariants for Hamiltonian systems subject
to holonomic constraints. The parametric $\alpha$-PRK methods are shown to
have the same convergence rate as the usual PRK methods and perform very well
in various numerical experiments.

##  1 Introduction

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 [BG, Tan] and discrete KAM theory
[Sha], 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, MQR],
average vector field method [QM], 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 [GM] 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 [CFM]
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 [BIT]. Consider the canonical
Hamiltonian systems in the form

$\begin{split}\begin{cases}&\dot{y}=J\nabla H(y),\\\
&y(t_{0})=y_{0}\in\mathbb{R}^{2d},\\\ \end{cases}\quad J=\begin{pmatrix}0&I\\\
-I&0\end{pmatrix}\in\mathbb{R}^{2d\times 2d},\end{split}$ (1.1)

where $y=(p^{T},q^{T})^{T}$ and $I$ is the identity matrix and $H$ is 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 $y_{1}(\alpha)=\Phi_{h}(y_{0},\alpha)$, where
$h$ is the step size of integration, $\alpha$ is a real parameter. This method
satisfies the following three conditions simultaneously: (i) for $\alpha=0$
one gets the Gauss collocation method of order $2s$, $s$ is the number of
stages of the RK method; (ii) for any fixed choice of $\alpha\neq 0$, the
corresponding method is of order $2s-2$ and satisfies the conditions
$b_{i}a_{i,j}+b_{j}a_{j,i}=b_{i}b_{j}$, thus being a quadratic invariant-
preserving RK method; (iii) for any choice of $y_{0}$ and in a neighborhood of
$h$, there exists a value of parameter $\alpha^{\ast}=\alpha^{\ast}(y_{0},h)$
such that $H(y_{1})=H(y_{0})$(energy conservation). The resulting method
$y_{1}=\Phi_{h}(y_{0},\alpha^{*})$ has order $2s$, preserves the energy and
quadratic invariants 111For stand RK method $(c,A,b)$ with constant step size
$h>0$, the sufficient condition of symplecticity is given by
$b_{i}a_{i,j}+b_{j}a_{j,i}=b_{i}b_{j}$ for all $i,j=1,2,...,s$. 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 $\alpha$-RK methods $y_{1}(\alpha)=\Phi_{h}(y_{0},\alpha^{\ast})$
that preserve the energy, the parameter depends on the initial values, i.e.,
$\alpha^{\ast}=\alpha^{\ast}(y_{0},h)$. Then, $\alpha$-RK methods with
condition $b_{i}a_{i,j}+b_{j}a_{j,i}=b_{i}b_{j}$ 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 [BGI], which leads to a more refined theoretical analysis and a nice
practical implementation strategy for seeking the parameter $\alpha^{*}$.
Another extension form $\alpha$-RK to $\alpha$-PRK to make use of the additive
structure of Hamiltonian system is given in [WXL]. The key of constructing
parameterized RK or PRK methods in [BIT, WXL] lies in the so-called
$W$-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 $q$ only) [HLW06, Chap. VII]

$\displaystyle\widetilde{H}(p,q,\lambda)=H(p,q)+g(q)^{T}\lambda,$ (1.2)

where $H:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$ is the Hamiltonian
function without constraint and $g(\cdot):\mathbb{R}^{d}\to\mathbb{R}^{m}$ is
the constraint function. Here, $\lambda\in\mathbb{R}^{m}$ is the Lagrangian
multiplier. The ODE corresponding to (1.2) is given by (see [HLW06, VII.1.2])

$\begin{split}&\dot{q}=\nabla_{p}H(p,q)=:H_{p}(p,q),\\\
&\dot{p}=-\nabla_{q}H(p,q)-G^{T}(q)\lambda=:-H_{q}(p,q)-G^{T}(q)\lambda,\\\
&0=g(q),\end{split}$ (1.3)

where we introduced $G(q)=\left(\frac{\partial g_{i}}{\partial
q_{j}}\right)_{m\times d}$ to be the Jacobian matrix of $g(q)$ so that
$G^{T}(q)=\nabla_{q}g\in\mathbb{R}^{d\times m}.$ Note that we use the
convention $(\nabla_{q}g)_{ij}=\frac{\partial g_{j}}{\partial q_{i}}$ which is
commonly used in the community of fluid mechanics. Throughout this paper, we
assume that the matrix $G(q)$ has full rank and $G(q)H_{pp}(p,q)G(q)^{T}$ is
invertible, which will allow us to express $\lambda$ in terms of $(p,q)$. The
dynamics given by (1.3) is on the manifold

$\begin{split}\mathcal{M}=\left\\{(p,q);\quad
g(q)=0,~{}G(q)H_{p}(p,q)=0\right\\},\end{split}$ (1.4)

which is the cotangent bundle of the manifold given by
$\mathcal{Q}:=\\{q:g(q)=0\\}.$ See [HLW06, VII.1.2] for more details. The
constraint $G(q)H_{p}(p,q)=0$ 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) $H$ is preserved:

$\displaystyle\frac{d}{dt}H(p(t),q(t))=0.$ (1.5)

Secondly, the symplecticity is also preserved. Consider a map
$\varphi:\mathcal{M}\to\mathcal{M}$ and its tangent map
$\varphi^{\prime}:T_{x}\mathcal{M}\to T_{x}\mathcal{M}$. One says the map
$\varphi$ is symplectic if

$\displaystyle(\varphi^{\prime}(x)\xi_{1})^{T}J(\varphi^{\prime}(x)\xi_{2})=\xi_{1}^{T}J\xi_{2},~{}\forall\xi_{1},\xi_{2}\in
T_{x}\mathcal{M}.$ (1.6)

For $\mathcal{M}$ being a submanifold of $\mathbb{R}^{2d}$,
$\varphi^{\prime}=(\nabla_{x}\varphi)^{T}$ so that

$\varphi^{\prime}\xi=\xi\cdot\nabla_{x}\varphi=\left.\frac{d}{d\tau}\right|_{\tau=0}\varphi(\gamma(\tau))$

for any smooth curve $\gamma$ in $\mathcal{M}$ satisfying
$\gamma(0)=x,\gamma^{\prime}(0)=\xi$. 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 $\mathcal{M}$, one also tends to ask that the
energy preservation or the symplecticity is satisfied. This additional
requirement for staying on manifolds $\mathcal{M}$ 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 $\mathcal{M}$ exactly, but they are not Hamiltonian energy
converved. The second order energy-conserving method in the line integral
framework developed in [BGIW] is not symplectic in general and can not
conserve the the hidden constraint $G(q)H_{p}(p,q)=0$. The projected RK
methods are proposed in [WDQW], 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 [BIT, WXL]. 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 $\alpha$-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 $\alpha$-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 $s$-stage PRK methods denoted by
$\left(c,A,b;\widehat{A},\widehat{b}\right)$ for the Hamiltonian ODE (1.3) is
given by the following

$\displaystyle\begin{split}&p_{1}=p_{0}+h\sum_{i=1}^{s}\hat{b}_{i}\ell_{i},~{}~{}~{}q_{1}=q_{0}+h\sum_{i=1}^{s}b_{i}k_{i}.\\\
&\ell_{i}=-H_{q}(P_{i},Q_{i})-G^{T}(Q_{i})\Lambda_{i},~{}k_{i}=H_{p}(P_{i},Q_{i}),~{}g(Q_{i})=0,\\\
&P_{i}=p_{0}+h\sum_{j=1}^{s}\hat{a}_{ij}\ell_{j},~{}~{}~{}Q_{i}=q_{0}+h\sum_{j=1}^{s}a_{ij}k_{j}.\end{split}$
(2.1)

The equations considered here are autonomous, so the parameters $c$ for the
Runge-Kutta methods are absent in (2.1). See section 4.1.2 for more details.

######  Lemma 2.1.

[Jay96] For a given a PRK method, if the conditions for symplecticity

$\displaystyle\begin{split}&b_{i}=\hat{b}_{i},\\\
&b_{i}\hat{a}_{ij}+\hat{b}_{j}a_{ji}=b_{i}\hat{b}_{j},\end{split}$ (2.2)

is satisfied, then the method for Hamiltonian systems with holonomic
constraints is also symplectic.

For $s$-stage PRK methods, if $a_{ij},\hat{a}_{ij},b_{i},\hat{b}_{i}$ are
given, then there are $3s$ equations for the $3s$ unknowns
($P_{i},Q_{i},\Lambda_{i}$). However, for general symplectic PRK methods,
$(p_{1},q_{1})$ may not fall onto $\mathcal{M}$. To ensure that $q_{1}$
satisfies the constraint, it is suggested in [Jay96] to impose the conditions

$\displaystyle a_{1j}=0,~{}~{}~{}a_{sj}=b_{j},~{}~{}1\leq j\leq s,$ (2.3)

so that

$\displaystyle Q_{1}\equiv q_{0},~{}~{}~{}Q_{s}\equiv q_{1}.$ (2.4)

With this, one has correspondingly

$\displaystyle\hat{a}_{js}=0,~{}\hat{a}_{j1}=b_{1},~{}~{}j=1,2,...,s,$ (2.5)

from the symplectic condition given in Lemma 2.1. This requirement removes the
unknowns $Q_{1},\Lambda_{s}$ and the equations $g(Q_{1})=0$,
$Q_{1}=q_{0}+h\sum_{j=1}^{s}a_{ij}k_{j}$. Hence, the number of equations is
now equal to the number of unknowns
$\left(Q_{2},\ldots,Q_{s},P_{1},\ldots,P_{s},\Lambda_{1},\ldots,\Lambda_{s-1}\right)$.
In fact, it is shown in [Jay96] that these variables can be uniquely solved
when $h$ is small enough under some reasonable assumptions on $H$ and
$\left(c,A,b;\widehat{A},\widehat{b}\right)$. We remark that even though one
imposes $Q_{1}=q_{0}$, the Lagrange multiplier $\Lambda_{1}$ is still in the
system of equations. Lastly, $p_{1}$ and $\Lambda_{s}$ can be determined by
the hidden constraint condition $0=G(q_{1})\cdot H_{p}(p_{1},q_{1}).$

For the PRK methods suggested in [Jay96], we hence solve for the next two
nonlinear systems in turn at each step. The first one

$\begin{split}\begin{cases}&Q_{i}=q_{0}+h\sum\limits_{j=1}^{s}a_{ij}H_{p}(P_{j},Q_{j}),\\\
&P_{i}=p_{0}-h\sum\limits_{j=1}^{s-1}\hat{a}_{ij}\left(H_{q}(P_{j},Q_{j})-G^{T}(Q_{j})\Lambda_{j}\right),\\\
&0=g(Q_{i}).\\\ \end{cases}\end{split}$ (2.6)

with $Q_{1}\equiv q_{0},Q_{s}\equiv q_{1}.$ Then, we can solve
$(p_{1},\Lambda_{s})$ by

$\begin{cases}p_{1}&=p_{0}-h\sum\limits_{i=1}^{s-1}\hat{b}_{i}\left(H_{q}(P_{i},Q_{i})-G^{T}(Q_{i})\Lambda_{i}\right)-h\hat{b}_{s}\left(H_{q}(P_{s},Q_{s})-G^{T}(Q_{s})\Lambda_{s}\right),\\\
0&=G(q_{1})H_{p}(p_{1},q_{1}).\end{cases}$ (2.7)

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
$(q_{0},q_{1})\in\mathbb{R}^{d}\times\mathbb{R}^{d}$, one can implicitly
define quantities
$\left(\bar{p}_{0},\bar{p}_{1},\bar{Q}_{i},\bar{P}_{i},\dot{\bar{P}}_{i},\dot{\bar{Q}}_{i}\right)$
for $i=1,\ldots,s$ and $\bar{\Lambda}_{i}$ for $i=2,\ldots,s-1$ through the
following system of equations

$\displaystyle\begin{split}q_{1}&=q_{0}+h\sum_{j=1}^{s}b_{j}\dot{\bar{Q}}_{j},\quad\bar{p}_{1}=\bar{p}_{0}+h\sum_{j=1}^{s}\hat{b}_{j}\dot{\bar{P}}_{j},\\\
\bar{Q}_{i}&=q_{0}+h\sum_{j=1}^{s}a_{ij}\dot{\bar{Q}}_{j},\quad\bar{P}_{i}=\bar{p}_{0}+h\sum_{j=1}^{s}\hat{a}_{ij}\dot{\bar{P}}_{j},~{}~{}i=1,\ldots,s,\\\
\dot{\bar{Q}}_{i}&=H_{p}(\bar{Q}_{i},\bar{P}_{i}),~{}~{}i=1,\ldots,s;\dot{\bar{P}}_{1}=-H_{q}(\bar{Q}_{1},\bar{P}_{1}),~{}~{}\dot{\bar{P}}_{s}=-H_{q}(\bar{Q}_{s},\bar{P}_{s})\\\
\dot{\bar{P}}_{i}&=-H_{q}(P_{i},Q_{i})-G^{T}(Q_{i})\Lambda_{i},~{}~{}~{}0=g(Q_{i}),~{}~{}i=2,\ldots,s-1.\end{split}$
(2.8)

Note that the definitions of $\dot{\bar{P}}_{1}$ and $\dot{\bar{P}}_{s}$ above
are different from other $\dot{\bar{P}}_{i}$ (in fact,
$\dot{\bar{P}}_{i}=\dot{P}_{i}$ for $i\neq 1,s$ while there should be some
corrections for $\dot{\bar{P}}_{1}$ and $\dot{\bar{P}}_{s}$ to define
$\dot{P}_{1},\dot{P}_{s}$). Moreover, $\bar{p}_{1}$ and $\bar{p}_{0}$ are not
$p_{1},p_{0}$ 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

$\displaystyle
L_{h}(q_{0},q_{1})=h\sum_{i=1}^{s}b_{i}L(\bar{Q}_{i},\dot{\bar{Q}}_{i}),$
(2.9)

where $L(q,\dot{q})$ is the Lagrangian for the time-continuous dynamics. The
Euler-Lagrange equation for this discrete Lagrangian under the holonomic
constraints $(q_{0},q_{1})\in\mathcal{Q}\times\mathcal{Q}$ reads

$\displaystyle\begin{split}&p_{0}=-\frac{\partial L_{h}}{\partial
x}(q_{0},q_{1})+hG^{T}(q_{0})\cdot\lambda_{0},\\\ &p_{1}=\frac{\partial
L_{h}}{\partial y}(q_{0},q_{1})-hG^{T}(q_{1})\cdot\lambda_{1},\\\
&g(q_{1})=0,\\\ &G(q_{1})H_{p}(p_{1},q_{1})=0.\end{split}$ (2.10)

See [MW01, eqs. (3.5.2a-d)]. Note that though the Euler-Lagrange equations are
derived by fixing $q_{0},q_{1}$, one should regard $p_{0},q_{0}$ as knowns and
solve $(p_{1},q_{1})$ 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
$\alpha$-PRK later. For typical separable Hamiltonian

$H(p,q)=\frac{1}{2}p^{T}M^{-1}p+U(q),$

the traditional Shake method reads

$\begin{split}&q_{n+1}-2q_{n}+q_{n-1}=-h^{2}M^{-1}(U_{q}(q_{n})+G(q_{n})^{T}\Lambda_{n}),\\\
&0=g(q_{n+1})\end{split}$

Then, the momentum is $p_{n}=M(q_{n+1}-q_{n-1})/(2h).$ The corresponding
Hamiltonian formulation for the Shake method is

$\begin{split}&p_{n+1/2}=p_{n}-\frac{h}{2}(U_{q}(q_{n})+G(q_{n})^{T}\Lambda_{n}),\\\
&q_{n+1}=q_{n}+hM^{-1}p_{n+1/2},~{}~{}g(q_{n+1})=0,\\\
&p_{n+1}=p_{n+1/2}-\frac{h}{2}(U_{q}(q_{n+1})+G(q_{n+1})^{T}\Lambda_{n+1})\end{split}$

Unfortunately, so-defined $p_{n}$ may not satisfy the constraint. Anderson
[And83] proposed the Rattle algorithm. That is to determine $p_{n+1}$ by using
another multiplier $\mu_{n}$ which can be different from $\Lambda_{n+1}$:

$p_{n+1}=p_{n+1/2}-\frac{h}{2}(U_{q}(q_{n+1})+G(q_{n+1})^{T}\mu_{n}).$

Then, one uses the condition

$G(q_{n+1})\cdot M^{-1}p_{n+1}=0$

to determine $\mu_{n}$ and thus $p_{n+1}$. 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

$\displaystyle\begin{split}&p_{n+1/2}=p_{n}-\frac{h}{2}(\nabla_{q}H(p_{n+1/2},q_{n})+G^{T}(q_{n})\Lambda_{n}),\\\
&q_{n+1}=q_{n}+\frac{h}{2}(\nabla_{p}H(p_{n+1/2},q_{n})+\nabla_{p}H(p_{n+1/2},q_{n+1})),~{}~{}~{}~{}g(q_{n+1})=0,\\\
&p_{n+1}=p_{n+1/2}-\frac{h}{2}(\nabla_{q}H(p_{n+1/2},q_{n+1})+G^{T}(q_{n+1})\Lambda_{n+1})\\\
&G(q_{n+1})\cdot p_{n+1}=0.\end{split}$ (2.11)

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
$\mathcal{M}$ has the following merits. (a) preserving the numerical solutions
on the manifold $\mathcal{M}$ exactly; (b) it is symmetric, symplectic and
super convergent of order $2s-2$, the $s$ is the stage of PRK method. In
particular, for $s=3$, the coefficients for Lobatto IIIA-IIIB pairs are given
by

$\begin{split}\begin{array}[]{c|ccc}0&0&0&0\\\
\frac{1}{2}&\frac{5}{24}&\frac{1}{3}&-\frac{1}{24}\\\
1&\frac{1}{6}&\frac{2}{3}&\frac{1}{6}\\\
\hline\cr&\frac{1}{6}&\frac{2}{3}&\frac{1}{6}\\\
\end{array},\quad\quad~{}\begin{array}[]{c|ccc}0&\frac{1}{6}&-\frac{1}{6}&0\\\
\frac{1}{2}&\frac{1}{6}&\frac{1}{3}&0\\\ 1&\frac{1}{6}&\frac{5}{6}&0\\\
\hline\cr&\frac{1}{6}&\frac{2}{3}&\frac{1}{6}\\\ \end{array}.\end{split}$
(2.12)

##  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

$\displaystyle I(p,q)=q^{T}Dp,$ (3.1)

where $D$ 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 $g(q)=0$, the quadratic invaraints
must be of some particular forms. Instead, motivated by Noether’s theorem
[Noe71, Arn13], we will consider that

$\displaystyle\mathscr{D}:=\Big{\\{}D:g\left(e^{sD}q\right)=g(q),~{}~{}L\left(e^{sD}q,e^{sD}\dot{q}\right)=L(q,\dot{q}),~{}\text{for
all}~{}s\in\mathbb{R}\Big{\\}},$ (3.2)

Clearly, by Noether’s theorem, we have the following.

######  Lemma 3.1.

For any $D\in\mathscr{D}$, the quantity $q^{T}Dp$ 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 $q^{T}Dp$ for $D\in\mathscr{D}$.

######  Theorem 3.1.

The symplectic PRK schemes for Hamiltonian system with holonomic constraints
defined in (2.1) conserve all quadratic invariants of the form

$\displaystyle I(p,q)=q^{T}Dp,~{}~{}D\in\mathscr{D}.$ (3.3)

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 $L_{h}(q_{0},q_{1})$ is invariant under a one-
parameter group of transformation $\\{\mathfrak{g}_{s}:s\in\mathbb{R}\\}$:
$L_{h}\left(\mathfrak{g}_{s}(q_{0}),\mathfrak{g}_{s}(q_{1})\right)=L_{h}(q_{0},q_{1})$
for all $s\in\mathbb{R}$ and $(q_{0},q_{1})$. Then the corresponding
variational integrators for Lagrangian systems have the first integral in the
form of

$p_{n+1}^{T}a(q_{n+1})=p_{n}^{T}a(q_{n})$

where $a(q)=\frac{d}{ds}\mathfrak{g}_{s}(q)|_{s=0}$ is the generator for the
group. Furthermore, when there are holonomic constraints, besides the
conditions above, if moreover $\\{\mathfrak{g}_{s}:s\in\mathbb{R}\\}$ leaves
the constraints manifold $\mathcal{Q}$ invariant, then the claim still holds.

The first part can be found in [HLW06, sec VI.6]. The second part has been
remarked in [MW01, sec. 3.4.2]. For the convenience, we provide a direct
verification here.

###### Proof of the second part of Lemma 3.2.

Taking derivative on $s$ and setting $s=0$ in
$L_{h}\left(\mathfrak{g}_{s}(q_{0}),\mathfrak{g}_{s}(q_{1})\right)=L_{h}(q_{0},q_{1})$,
one has

$a(q_{0})\cdot\frac{\partial L_{h}}{\partial
x}(q_{0},q_{1})+a(q_{1})\cdot\frac{\partial L_{h}}{\partial
y}(q_{0},q_{1})=0.$

Since $\\{\mathfrak{g}_{s}:s\in\mathbb{R}\\}$ leaves the constraints manifold
invaraint,

$G(q)a(q)=0\Rightarrow a(q)\cdot G^{T}(q)\lambda=0,~{}\text{for
all~{}}\lambda.$

Using (2.10), one easily finds $p_{1}^{T}a(q_{1})=p_{0}^{T}a(q_{0})$. ∎

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.

######  Proposition 3.1.

The discrete Lagragian for
$(q_{0},q_{1})\in\mathbb{R}^{d}\times\mathbb{R}^{d}$ for the symplectic PRK
methods ((2.1) with conditions (2.2)-(2.3)) can be defined as

$\displaystyle
L_{h}(q_{0},q_{1})=\mathrm{ext}_{\left\\{\dot{Q}_{i}\right\\}}h\sum_{i=1}^{s}b_{i}L(Q_{i},\dot{Q}_{i})$
(3.4)

with constraints

$\displaystyle g(Q_{i})=0,~{}~{}i=2,\ldots,s-1,$ (3.5) $\displaystyle
q_{1}=q_{0}+h\sum_{j=1}^{s}a_{ij}\dot{Q}_{j}.$ (3.6)

Above, “$\mathrm{ext}$” means extremizing and
$Q_{i}=q_{0}+h\sum_{j=1}^{s}a_{ij}\dot{Q}_{j}$.

######  Remark 3.1.

Note that we are not defining $Q_{1},Q_{s}$ as
$a_{si}=b_{i},~{}~{}a_{1i}=0,b_{i}\neq 0$, so that $Q_{1}=q_{0}$,
$Q_{s}=q_{1}$. In this definition, $q_{0},q_{1}$ do not have to be in
$\mathcal{Q}$.

###### 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

$\displaystyle\begin{split}&\left(\sum_{i=1}^{s}b_{i}\frac{\partial
L}{\partial q}(Q_{i},\dot{Q}_{i})ha_{ij}\right)+b_{j}\frac{\partial
L}{\partial\dot{q}}(Q_{j},\dot{Q}_{j})=\sum_{i=2}^{s-1}b_{i}ha_{ij}G^{T}(Q_{i})\Lambda_{i}+b_{j}\lambda,\\\
&g(Q_{i})=0,~{}i=2,\ldots,s-1,\\\
&q_{1}=q_{0}+h\sum_{j=1}^{s}a_{ij}\dot{Q}_{j}.\end{split}$ (3.7)

Note that we have built in the constraint $g(Q_{i})=0$ by using the multiplier
$b_{i}\Lambda_{i}$. Since $b_{i}\neq 0$, this will not change anything. These
$sd$ equations together with $g(Q_{i})=0$ ($s-2$ equations) and the constraint
($d$ equations) can yield the solutions
$(\\{\dot{Q}_{j}\\}_{j=1}^{s},\\{\Lambda_{i}\\}_{i=2}^{s-1},\lambda)$ for $h$
sufficiently small.

Here we are solving $\dot{Q}_{i}$ in terms of $q_{0},q_{1}$ for any
$(q_{0},q_{1})$ in $\mathbb{R}^{d}\times\mathbb{R}^{d}$, not just for
$(q_{0},q_{1})\in\mathcal{Q}\times\mathcal{Q}$. With this observation, we can
take differentiation on $q_{\ell},\ell=0,1$ for the constraint (3.6) to find

$\displaystyle h\sum_{j=1}^{s}b_{j}\frac{\partial\dot{Q}_{j}}{\partial
q_{0}}=-I,~{}~{}~{}h\sum_{j=1}^{s}b_{j}\frac{\partial\dot{Q}_{j}}{\partial
q_{1}}=I,$ (3.8)

and for the constraints (3.5) to find

$\displaystyle h\left(\sum_{j=1}^{s}a_{ij}\frac{\partial\dot{Q}_{j}}{\partial
q_{0}}\right)G^{T}(Q_{i})=-G^{T}(Q_{i}),~{}~{}~{}h\left(\sum_{j=1}^{s}a_{ij}\frac{\partial\dot{Q}_{j}}{\partial
q_{1}}\right)G^{T}(Q_{i})=0.$ (3.9)

Note that here we used the convention

$\frac{\partial\dot{Q}_{j}}{\partial
q_{0}}:=\nabla_{q_{0}}\dot{Q}_{j},~{}~{}~{}\left(\frac{\partial\dot{Q}_{j}}{\partial
q_{0}}\right)_{\ell,m}=\frac{\partial\dot{Q}_{j}^{(m)}}{\partial
q_{0}^{(\ell)}}.$

The reason to do this is that we have made $\frac{\partial L}{\partial q}$ a
column vector.

Step–2: We derive the equation for $p_{i}$ defined in (2.10). By the first
equation in (2.10), one has

$\displaystyle\begin{split}p_{0}&=-\left(h\sum_{i=1}^{s}b_{i}\left[I+h\sum_{j=1}^{s}a_{ij}\frac{\partial\dot{Q}_{j}}{\partial
q_{0}}\right]\frac{\partial L}{\partial
q}(Q_{i},\dot{Q}_{i})+h\sum_{i=1}^{s}b_{i}\frac{\partial\dot{Q}_{i}}{\partial
q_{0}}\frac{\partial
L}{\partial\dot{q}}(Q_{i},\dot{Q}_{i})\right)+hG^{T}(q_{0})\cdot\lambda_{0}\\\
&=-h\sum_{i=1}^{s}b_{i}\frac{\partial L}{\partial
q}(Q_{i},\dot{Q}_{i})-h\sum_{j=1}^{s}\frac{\partial\dot{Q}_{j}}{\partial
q_{0}}\left[\sum_{i=2}^{s-1}hb_{i}a_{ij}G^{T}(Q_{i})\Lambda_{i}+b_{j}\lambda\right]+hG^{T}(q_{0})\cdot\lambda_{0}\\\
&=-h\sum_{i=1}^{s}b_{i}\frac{\partial L}{\partial
q}(Q_{i},\dot{Q}_{i})+h\sum_{i=2}^{s-1}b_{i}G^{T}(Q_{i})\Lambda_{i}+\lambda+hG^{T}(q_{0})\cdot\lambda_{0}.\end{split}$

Above, the second inequality follows from (3.7) while last inequality follows
from (3.8) and (3.9). Similarly, we can compute by the second equation in
(2.10) that

$\displaystyle\begin{split}p_{1}&=h\sum_{i=1}^{s}b_{i}h\sum_{j=1}^{s}a_{ij}\frac{\partial\dot{Q}_{j}}{\partial
q_{1}}\frac{\partial L}{\partial
q}(Q_{i},\dot{Q}_{i})+h\sum_{i=1}^{s}b_{i}\frac{\partial\dot{Q}_{i}}{\partial
q_{1}}\frac{\partial
L}{\partial\dot{q}}(Q_{i},\dot{Q}_{i})-hG^{T}(q_{1})\cdot\lambda_{1}\\\
&=h\sum_{j=1}^{s}\frac{\partial\dot{Q}_{j}}{\partial
q_{1}}\left[\sum_{i=2}^{s-1}hb_{i}a_{ij}G^{T}(Q_{i})\Lambda_{i}+b_{j}\lambda\right]-hG^{T}(q_{1})\cdot\lambda_{1}\\\
&=\lambda-hG^{T}(q_{1})\cdot\lambda_{1}.\end{split}$

Defining that $\Lambda_{1}:=\lambda_{0}/b_{1},\Lambda_{s}:=\lambda_{1}/b_{s},$
one has

$\displaystyle p_{1}=p_{0}+h\sum_{i=1}^{s}b_{i}\frac{\partial L}{\partial
q}(Q_{i},\dot{Q}_{i})-h\sum_{i=1}^{s}b_{i}G^{T}(Q_{i})\Lambda_{i}.$ (3.10)

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 $H(p,q)=p\cdot\dot{q}-L(q,\dot{q}),~{}~{}p=\frac{\partial
L}{\partial\dot{q}},$ where $\dot{q}=\dot{q}(p,q)$ is determined by the second
equation. It can be computed that $\dot{q}=\frac{\partial H(p,q)}{\partial
p},~{}~{}~{}\frac{\partial H(p,q)}{\partial q}=-\frac{\partial
L(q,\dot{q})}{\partial q}.$ Hence, we have

$\displaystyle
Q_{i}=q_{0}+h\sum_{j=1}^{s}a_{ij}H_{p}(P_{i},Q_{i}),~{}~{}~{}q_{1}=q_{0}+h\sum_{j=1}^{s}a_{ij}H_{p}(P_{i},Q_{i}),$
(3.11) $\displaystyle
p_{1}=p_{0}-h\sum_{i=1}^{s}b_{i}H_{q}(P_{i},Q_{i})-h\sum_{i=1}^{s}b_{i}G^{T}(Q_{i})\Lambda_{i},$
(3.12)

where $P_{i}=\frac{\partial L(Q_{i},\dot{Q}_{i})}{\partial\dot{Q}_{i}}.$ The
first equation in (3.7) can then be rewritten as

$\sum_{i=1}^{s}hb_{i}a_{ij}(-H_{q}(P_{i},Q_{i}))+b_{j}P_{j}=\sum_{i=2}^{s-1}b_{i}ha_{ij}G^{T}(Q_{i})\Lambda_{i}+b_{j}[p_{1}+hb_{s}G^{T}(q_{1})\cdot\Lambda_{s}]$

Using (3.12), we then have

$\displaystyle
P_{j}=p_{0}-h\sum_{i=1}^{s}[b_{i}-b_{i}a_{ij}/b_{j}]H_{q}(P_{i},Q_{i})-\sum_{i=2}^{s-1}[b_{i}-b_{i}a_{ij}/b_{j}]G^{T}(Q_{i})\Lambda_{i}-hb_{1}G^{T}(Q_{1})\cdot\Lambda_{1}.$
(3.13)

Using the condition $a_{1j}=0,a_{sj}=b_{j}$, this is further simplified to

$\displaystyle
P_{j}=p_{0}-h\sum_{i=1}^{s}[b_{i}-b_{i}a_{ij}/b_{j}]H_{q}(P_{i},Q_{i})-\sum_{i=1}^{s}[b_{i}-b_{i}a_{ij}/b_{j}]G^{T}(Q_{i})\Lambda_{i}.$
(3.14)

Defining $\hat{a}_{ji}=b_{i}-b_{i}a_{ij}/b_{j},$ one has the desired form for
$P_{j}$.

Eventually, we find (3.11), (3.12), (3.14) and (3.5)–(3.6) then form a
complete system that is the same as the symplectic PRK method. ∎

Lastly, we verify Theorem 3.1.

###### Proof of Theorem 3.1.

The constraint condition is obvious. We verify that

$\displaystyle L_{h}(e^{sD}q_{0},e^{sD}q_{1})=L_{h}(q_{0},q_{1}).$ (3.15)

This is in fact a direct corollary of Proposition 3.1. Since $e^{sD}$ is a
linear group, then for any $(q_{0},q_{1})$ and any
$\\{\dot{Q}_{j}\\}_{j=1}^{s}$, we have a corresponding sequence of data
$\\{\dot{\bar{Q}}_{j}\\}_{j=1}^{s}:=\\{e^{sD}\dot{Q}_{j}\\}_{j=1}^{s}$ for
$(e^{sD}q_{0},e^{sD}q_{1})$. Correspondingly the new $Q_{i}$ (denoted by
$\bar{Q}_{i}$) is given by

$\bar{Q}_{i}=e^{sD}q_{0}+h\sum_{j=1}^{s}a_{ij}e^{sD}\dot{Q}_{j}=e^{sD}Q_{i},$

by the linearity of $e^{sD}:\mathbb{R}^{d}\to\mathbb{R}^{d}$. Since for the
continuous Lagrangian one has
$L(e^{sD}Q_{i},e^{sD}\dot{Q}_{i})=L(Q_{i},\dot{Q}_{i})$, we find the sums to
extremize in (3.4) have the same value at the corresponding data. Hence, their
extrema must be the same and thus (3.15) follows. Applying the discrete
Noether’s theorem (Lemma 3.2), we obtain the desired result. ∎

##  4 Energy and quadratic invariants preserving methods

In this section, we propose a method that preserve the energy and quadratic
invariants for Hamiltonian system with constraints based on the symplectic PRK
schemes, following the spirit of [BIT, WXL]. Recall that the following
conditions guarantee that the $s$-stage PRK method is symplectic and preserves
the quadratic invariants:

$\displaystyle\begin{split}&\text{The symplecticity and quadratic
invariants:~{}}b_{i}=\hat{b}_{i},b_{i}\hat{a}_{ij}+\hat{b}_{j}a_{ji}=b_{i}\hat{b}_{j},i,j=1,...,s.\\\
&\text{For the constraints:~{}}a_{1j}=0,~{}a_{sj}=b_{j},\text{or
equivalently,~{}}\hat{a}_{is}=0,~{}~{}\hat{a}_{i1}=b_{1},~{}i=1,\ldots,s.\end{split}$
(4.1)

These schemes do not in general preserve the Hamiltonian energy. In fact, such
methods conserve the symplectic properties and Hamiltonian energy for constant
stepsize $h>0$ has been proven not to exist in general [CFM, GM]. The key
observation now is that the conditions imposed above for preserving the
symplecticity and quadratic invariants and constraints manifolds cannot
determine $\left(c,A,b;\hat{A},\hat{b}\right)$ totally, so there is actually
some freedom to choose these parameters.

Inspired by [BIT, WXL], we are motivated to construct $\alpha$-PRK methods so
that the energy is also preserved. However, the specific construction
technique developed in [BIT, WXL] is not feasible here in general. In fact,
the key in the construction of $\alpha$-PRK methods in the two paper is the so
called $W$-transformation, which allows us to express the symplectic condition
in a very compact form with the help of matrix $X_{G}$, where the three
diagonal anti-symmetric matrix $X_{G}$ is the $W$-transformation matrix of
Gauss collocation method. Thus, we can parameterize the transformed matrix so
that it remains to satisfy the same conditions as usual symplectic RK or PRK
methods, and then obtain the parameterized $\alpha$-PRK method by inverse
$W$-transformation. For example, the $2$-stage and $3$-stage $\alpha$-PRK
methods based on Lobatto IIIA-IIIB pairs are given by [WXL]

$\begin{split}\begin{array}[]{c|cc}0&\frac{1}{6}\alpha&-\frac{1}{6}\alpha\\\
1&\frac{1}{2}-\frac{1}{6}\alpha&\frac{1}{2}+\frac{1}{6}\alpha\\\
\hline\cr&\frac{1}{2}&\frac{1}{2}\\\
\end{array},\quad\quad~{}\begin{array}[]{c|cc}0&\frac{1}{2}-\frac{1}{6}\alpha&\frac{1}{6}\alpha\\\
1&\frac{1}{2}+\frac{1}{6}\alpha&-\frac{1}{6}\alpha\\\
\hline\cr&\frac{1}{2}&\frac{1}{2}\\\ \end{array}.\end{split}$ (4.2)

$\begin{split}\begin{array}[]{c|ccc}0&\frac{1}{5}\alpha&-\frac{2}{5}\alpha&\frac{1}{5}\alpha\\\
\frac{1}{2}&\frac{5}{24}-\frac{1}{15}\alpha&\frac{1}{3}+\frac{2}{15}\alpha&-\frac{1}{24}-\frac{1}{15}\alpha\\\
1&\frac{1}{6}+\frac{1}{5}\alpha&\frac{2}{3}-\frac{2}{5}\alpha&\frac{1}{6}-\frac{1}{5}\alpha\\\
\hline\cr&\frac{1}{6}&\frac{2}{3}&\frac{1}{6}\\\
\end{array},\quad\quad~{}\begin{array}[]{c|ccc}0&\frac{1}{6}-\frac{1}{5}\alpha&-\frac{1}{6}+\frac{2}{5}\alpha&-\frac{1}{5}\alpha\\\
\frac{1}{2}&\frac{1}{6}+\frac{1}{15}\alpha&\frac{1}{3}-\frac{2}{15}\alpha&\frac{1}{15}\alpha\\\
1&\frac{1}{6}-\frac{1}{5}\alpha&\frac{5}{6}+\frac{2}{5}\alpha&-\frac{1}{5}\alpha\\\
\hline\cr&\frac{1}{6}&\frac{2}{3}&\frac{1}{6}\\\ \end{array}.\end{split}$
(4.3)

We can see that although the method keeps symplecticity and quadratic
invariants, it generally does not satisfy the constraint conditions
$a_{1j}=0,a_{sj}=b_{j}$, so its numerical solution cannot be guaranteed to
fall on the manifold $\mathcal{M}$.

In this article, we introduce another new parameterization $\alpha$-PRK
methods. According to (4.1), we focus only on $(A,b)$. Note that the first and
last rows of $A$ are fixed and depended on $b$ while the other rows of $A$ can
be determined in a great freedom. Since the energy $H(p,q)$ is a scalar, we
can let these coefficients change depending on a scalar $\alpha\in\mathbb{R}$
in a particular chosen way. Of course, due to the degrees of freedom, such a
choice may not be unique. Then, we hope to choose $\alpha$ such that the
conditions (4.1) are satisfied. Let’s first take a closer look at classic
examples and then introduce the general higher-order methods.

###  4.1 New modified algorithms based on existing methods

####  4.1.1 $\alpha$-Rattle

In this section, we aim to propose a method based on the Shake-Rattle
algorithm in section 2.1. As known, this method is two stage method with order
$2$. If we perturb $b$, the order can be reduced. Since $s=2$, with the
conditions (4.1), there is only one degree of freedom left for choosing the
parameters. In particular, we let

$\displaystyle
b_{1}=1/2+\alpha,~{}~{}b_{2}=1/2-\alpha,~{}~{}\alpha\in\mathbb{R}.$ (4.4)

Then, the matrix $A(\alpha)$-$\widehat{A}(\alpha)$ is fully determined, and
the $\alpha$-PRK, which we call $\alpha$-Rattle, in the form Butcher table is
given:

$\begin{split}\begin{array}[]{c|cc}0&0&0\\\
1&\frac{1}{2}+\alpha&\frac{1}{2}-\alpha\\\
\hline\cr&\frac{1}{2}+\alpha&\frac{1}{2}-\alpha\\\
\end{array},\quad\quad~{}\begin{array}[]{c|cc}\frac{1}{2}+\alpha&\frac{1}{2}+\alpha&0\\\
\frac{1}{2}+\alpha&\frac{1}{2}+\alpha&0\\\
\hline\cr&\frac{1}{2}+\alpha&\frac{1}{2}-\alpha\\\ \end{array}.\end{split}$
(4.5)

We can see that this new $\alpha$-PRK method is different from the previous
method given in (4.2). For fixed $\alpha$, it not only can preserve the
symplecticity and quadratic invariants, but also can guarantee the numerical
solutions to fall on the manifold $\mathcal{M}$. When we adjust $\alpha$ each
step, the method can preserve energy and quadratic invariants, and also to
guarantee the numerical solutions to fall on the manifold $\mathcal{M}$.

We write out the details of the $\alpha$-Rattle method below for convenience:

$\begin{split}\begin{cases}&P_{1}=p_{0}+h(1/2+\alpha)\left(-H_{q}(P_{1},Q_{1})-G^{T}(Q_{1})\Lambda_{1}\right),\\\
&P_{2}=P_{1},\\\ &Q_{1}=q_{0},\\\
&Q_{2}=q_{0}+h\left((\frac{1}{2}+\alpha)H_{p}(P_{1},Q_{1})+(\frac{1}{2}-\alpha)H_{p}(P_{2},Q_{2})\right),\\\
&g(Q_{1})=g(Q_{2})=0,\\\ \end{cases}\end{split}$ (4.6)

The numerical solution $(p_{1},q_{1})$ is then solved by

$\begin{cases}q_{1}&=Q_{2},\\\
\ell_{1}&=-H_{q}(P_{1},Q_{1})-G^{T}(Q_{1})\Lambda_{1},\\\
\ell_{2}&=-H_{q}(P_{2},Q_{2})-G^{T}(Q_{2})\Lambda_{2},\\\
p_{1}&=p_{0}+h\left((\frac{1}{2}+\alpha)\ell_{1}+(\frac{1}{2}-\alpha)\ell_{2}\right),\\\
0&=G(q_{1})H_{p}(p_{1},q_{1}).\end{cases}$ (4.7)

The system (4.6) can be reduced to that

$\begin{split}\begin{cases}&P=p_{0}-h(1/2+\alpha)H_{q}(P,q_{0})-hG^{T}(q_{0})\lambda_{1},\\\
&q_{1}=q_{0}+h\left((\frac{1}{2}+\alpha)H_{p}(P,q_{0})+(\frac{1}{2}-\alpha)H_{p}(P,q_{1})\right),\\\
&0=g(q_{1}).\\\ \end{cases}\end{split}$ (4.8)

where $\lambda_{1}:=(1/2+\alpha)\Lambda_{1}$. There are three unknowns in this
system. Then, one can solve $(p_{1},\lambda_{2})$ by

$\begin{cases}p_{1}&=p_{0}-h(\frac{1}{2}+\alpha)H_{q}(P,q_{0})-h(\frac{1}{2}-\alpha)H_{q}(P,q_{1})-hG^{T}(q_{0})\lambda_{1}-hG^{T}(q_{1})\lambda_{2},\\\
0&=G(q_{1})H_{p}(p_{1},q_{1}),\end{cases}$ (4.9)

where $\lambda_{2}:=(1/2-\alpha)\Lambda_{2}$.

####  4.1.2 $\alpha$-PRK III methods

We now construct the $\alpha$-PRK III methods from Lobatto IIIA or IIIB
methods. For unconstrained Hamiltonian systems, this method has been used by
Sun [Sun95] to construct symplectic PRK method. However, as we pointed out
earlier, this method generally does not conserve the Hamiltonian energy.

The general Runge-Kutta method with parameters $(c,A,b)$ for non-autonomous
equations $\dot{y}=f(t,y)$ is given by

$y_{1}=y_{0}+h\sum_{j=1}^{s}b_{j}f(t_{0}+c_{j}h,k_{j}),~{}~{}~{}~{}k_{i}=y_{0}+h\sum_{j=1}^{s}a_{ij}f(t_{0}+c_{j}h,k_{j}),~{}i=1,\ldots,s.$

The classical $B(p),C(\eta)$ and $D(\zeta)$ conditions for the RK method with
parameters $(c,A,b)$ are given by

$\displaystyle\begin{split}B(p):&\quad\sum_{i=1}^{s}b_{i}c_{i}^{q-1}=\frac{1}{q}\quad
q=1,\ldots,p;\\\
C(\eta):&\quad\sum_{j=1}^{s}a_{ij}c_{j}^{q-1}=\frac{c_{i}^{q}}{q}\quad
i=1,\ldots,s,~{}q=1,\ldots,\eta;\\\
D(\zeta):&\quad\sum_{i=1}^{s}b_{i}c_{i}^{q-1}a_{ij}=\frac{b_{j}}{q}(1-c_{j}^{q})\quad
i=1,\ldots,s,~{}q=1,\ldots,\zeta.\end{split}$ (4.10)

See [HW06, IV.5, pp71] for more details on the $B(p),C(\eta)$ and $D(\zeta)$
assumptions for RK methods.

######  Theorem 4.1.

[Sun95] Suppose that an $s$-stage RK method $(c,A,b)$ with distinct abscissae
$c_{i}$ and $b_{i}\neq 0$ satisfies $B(p),C(\eta)$ and $D(\zeta)$. Then the
PRK method generated by
$\left(c_{i},a_{ij},\hat{a}_{ij}=b_{j}(1-\frac{a_{ji}}{b_{i}}),b_{i}\right)$
is symplectic and of order $q=\mathrm{min}(p,2\eta+2,2\zeta+2,\eta+\zeta+1)$
for Hamiltonian systems.

As a remark, the equations considered in this work are autonomous and the
parameters $c$ do not appear in the scheme. However, the conditions
$B(p),C(\eta),D(\zeta)$ still make sense: one may find the parameters $c$ from
$C(\eta)$ condition and then require $B(p),D(\zeta)$ to hold.

Here we construct several specific parameterized $\alpha$-PRK III methods for
Hamiltonian systems with holonomic constraints following the idea described in
Algorithm 1.

One method to parameterize the Lobatto IIIA or IIIB methods is to keep $b$ to
be constant at first but to change $a_{i2}$ in some suitable ways. This leads
to that $3$-stage parameterized $\alpha$-Lobatto IIIA-III$\widehat{A}$

$\begin{split}\begin{array}[]{c|ccc}0&0&0&0\\\
\frac{1}{2}&\frac{5}{24}-\alpha&\frac{1}{3}-\alpha&2\alpha-\frac{1}{24}\\\
1&\frac{1}{6}&\frac{2}{3}&\frac{1}{6}\\\
\hline\cr&\frac{1}{6}&\frac{2}{3}&\frac{1}{6}\\\
\end{array},\quad\quad~{}\begin{array}[]{c|ccc}0&\frac{1}{6}&4\alpha-\frac{1}{6}&0\\\
\frac{1}{2}&\frac{1}{6}&\frac{1}{3}+\frac{2}{3}\alpha&0\\\
1&\frac{1}{6}&\frac{5}{6}-8\alpha&0\\\
\hline\cr&\frac{1}{6}&\frac{2}{3}&\frac{1}{6}\\\ \end{array}.\end{split}$
(4.11)

Or the $3$-stage parameterized $\alpha$-Lobatto III$\widehat{B}$-IIIB

$\begin{split}\begin{array}[]{c|ccc}0&0&0&0\\\
\frac{1}{2}&\frac{5}{24}-\frac{\alpha}{2}&\frac{1}{3}+\alpha&-\frac{1}{24}-\frac{\alpha}{2}\\\
1&\frac{1}{6}&\frac{2}{3}&\frac{1}{6}\\\
\hline\cr&\frac{1}{6}&\frac{2}{3}&\frac{1}{6}\\\
\end{array},\quad\quad~{}\begin{array}[]{c|ccc}0&\frac{1}{6}&2\alpha-\frac{1}{6}&0\\\
\frac{1}{2}&\frac{1}{6}&\frac{1}{3}-\alpha&0\\\
1&\frac{1}{6}&\frac{5}{6}+2\alpha&0\\\
\hline\cr&\frac{1}{6}&\frac{2}{3}&\frac{1}{6}\\\ \end{array}.\end{split}$
(4.12)

One can check that all the $\alpha$-PRK III methods used in (4.11)-(4.12)
satisfy that $B(2)$ and $C(1)$. Hence, those four $\alpha$-PRK methods have
second order accuracy according to Theorem 4.1. We can construct general
$s$-stage $\alpha$-PRK methods based on IIIA-IIIB methods in the same way.
Note that those methods are different from the previous $\alpha$-PRK based on
$W$-transformation method given in (4.3). When $\alpha=0$, all the
$\alpha$-PRK methods presented in (4.3), (4.11)-(4.12) are reduced to the
classical Lobatto IIIA-IIIB pairs.

###  4.2 The general method and some theoretical results

As we can see from the above example, we can first parameterize the
coefficients $b(\alpha),A(\alpha)$, and then determine the entire $\alpha$-PRK
method according to the conditions (4.1).

####  4.2.1 The general method and properites

With the discussion above, we propose the method as follows:

Algorithm 1 (Energy and quadratic invariant preserving method)

1:Choose a particular form of functions $b=b(\alpha)$, and $A=A(\alpha)$ such
that there exists an interval $I\subset\mathbb{R}$ satisfying

  * •

For all $1\leq i\leq s$, $b_{i}(\alpha)\neq 0$ for all most all $\alpha\in I$;
$\sum_{i}b_{i}(\alpha)=1$.

  * •

$a_{1j}(\alpha)=0$, $a_{sj}(\alpha)=b_{j}(\alpha)$ for all $1\leq j\leq s$ and
all $\alpha\in I$.

  * •

Determine $\widehat{b}$ and $\widehat{A}$ by relation (4.1).

2:Pick computational time $T>0$; time step $h>0$. Let $N=\lceil T/h\rceil$.
Pick initial data $(p_{0},q_{0})$.

3:for $n\text{ in }1:N$ do

4: Let $(p_{n}(\alpha),q_{n}(\alpha))$ be given by (2.1) with coefficients
above and initial data $(p_{n-1},q_{n-1})$.

5: Choose $\alpha^{*}\in I$ such that
$H(p_{1}(\alpha^{*}),q_{1}(\alpha^{*}))=H(p_{0},q_{0})$.

6: Let $(p_{n},q_{n})=(p_{n}(\alpha^{*}),q_{n}(\alpha^{*}))$.

7:end for

From the algorithm, we can see that in order to preserve the constrained
manifold $\mathcal{M}$ we can impose the conditions
$a_{1j}=0,~{}a_{sj}=b_{j},$ or equivalently,
$\hat{a}_{is}=0,~{}~{}\hat{a}_{i1}=b_{1},~{}i=1,\ldots,s,$ and then to
determine the symplectic PRK methods form the relation (4.1).

For the general $s$-stage ($s\geq 3$) RK methods, due to the increase of
degrees of freedom, the parameterization methods is not unique, and there may
be many other methods. We now explore the properties of the new methods.

######  Theorem 4.2.

If the method in Algorithm 1 has a solution, then the energy $H(p_{n},q_{n})$
is a constant. Moreover, for all $D\in\mathscr{D}$, the quadratic invariant
$q_{n}^{T}Dp_{n}$ is also a constant.

###### Proof.

The energy constant is obvious by the method. For the quadratic invariant, we
note that Theorem 3.1 tells us that for any $\alpha$,
$q_{n}^{T}(\alpha)Dp_{n}(\alpha)=q_{0}Dp_{0}$. Hence, the claims follow. We
skip the details. ∎

Here, as that for symplectic RK methods, let us point out the difference
between the conservation of quadratic invariants and symplecticity. As we have
seen, for PRK methods with fixed coefficients, (4.1) guarantees the
symplecticity and conservation of quadratic invariants. In fact, for
Hamiltonian system, conservation of quadratic invariants and symplecticity are
equivalent. Now, if we apply the method for all initial data $(p_{0},q_{0})$
(all trajectories), the method given by Theorem 4.2 may not be symplectic
although (4.1) is satisfied. The reason is that the coefficients in the method
depend on the initial data now, and thus $dp_{n}\wedge dq_{n}$ may not be
equal to $dp_{0}\wedge dq_{0}$. Compared to the fact that symplectic involves
the differential of numerical solutions, the quadratic invariants depend only
on the values of the numerical solutions. The conservation of quadratic
invariants still holds by relations (4.1). Based on the discussion above, our
method can be applied in two possible cases.

  1. 1.

We apply the method for all trajectories so that the energy and quadratic
invariants are preserved for all trajectories. However, the method may not be
symplectic though.

  2. 2.

We fix down the coefficients obtained by a particular trajectory. Then, the
method is symplectic, preserves the quadratic invariants, and preserves the
energy for the chosen trajectory.

####  4.2.2 Existence of the numerical solutions and order of the methods

Clearly, the method in Algorithm 1 depends crucially on the solvability of

$\displaystyle\mu(\alpha,h):=H(p_{1}(\alpha,h),q_{1}(\alpha,h))-H(p_{0},q_{0})=0,$
(4.13)

for small enough $h>0$ and some given consistent initial values
$(p_{0},q_{0})\in\mathcal{Q}\times\mathcal{Q}$. Clearly, $\mu(\alpha,h)$ is
the error measuring the error for energy. As $h\to 0^{+}$, $\mu(\alpha,h)\to
0$ for all $\alpha$. We assume the order of the energy satisfies the follows.

######  Assumption 4.1.

We assume the method for $\alpha=0$ is a special base method that satisfies
the following.

  1. 1.

Suppose $\mu(\alpha,h)$ is analytic in a neighborhood of $(0,0)$.

  2. 2.

The order of energy approximation is higher for $\alpha=0$ than other $\alpha$
values near $0$. In other words, there exist $r,r_{1}$ with $r\geq r_{1}>1$
such that

$\displaystyle\lim_{h\to 0^{+}}\frac{1}{h^{r}}\mu(0,h)\neq
0,\quad~{}\lim_{h\to 0^{+}}\frac{1}{h^{r_{1}}}\mu(\alpha,h)\neq 0.$ (4.14)

######  Remark 4.1.

To understand the second term in Assumption 4.1, we note that

$\mu(\alpha,h)=H(p_{1}(\alpha,h),q_{1}(\alpha,h))-H(p(h),q(h))$

where $(p(t),q(t))$ is the exact solution of the Hamiltonian system with
initial data $(p_{0},q_{0})$. Hence, $\mu$ is in fact the numerical error of
the energy for the scheme. The second iterm in Assumption 4.1 is the
convergence order of the energy, which is often $O(h^{\beta+1})$ if the method
has a convergence order $\beta$.

By the analyticity, one has in a neighborhood of $(0,0)$ that

$\mu(\alpha,h)=\sum_{\ell\geq 0,m>1}c_{\ell,m}\alpha^{\ell}h^{m}$

Since $\mu$ cannot be identically zero, there must some $\ell,m$ such that
$c_{\ell,m}\neq 0$. Hence, there always exists some $r(\alpha)>0$ such that
$\lim_{h\to 0^{+}}\mu(\alpha,h)/h^{r(\alpha)}$ exists and is nonzero. Hence,
the second item in the assumption roughly says we have

$\displaystyle\mu(\alpha,h)=\left[c_{0,r}h^{r}+\sum_{m>r}c_{0,m}h^{m}\right]+c_{kr_{1}}\alpha^{k}h^{r_{1}}+\cdots,~{}~{}\quad
k>0,r\geq r_{1}>1.$ (4.15)

We further assume that the dependence in $\alpha$ is not degenerate, and
without los of generalization in (4.15) for $k=1$. If $k\neq 1$, we may set
$\tilde{\alpha}=\alpha^{k}$. Of course, in terms of $(\tilde{\alpha},h)$, the
analyticity of $\mu$ may break down. However, as soon as $\mu$ is $C^{1}$ in
$\tilde{\alpha}$ and $h$, the claim below still holds. Provided $k=1$, we can
prove that he method can give a sequence of numerical solutions, as below.

######  Theorem 4.3.

Suppose Assumption 4.1 holds and $r,r_{1}$ are the numbers as in Assumption
4.1. If it holds that

$\displaystyle\lim_{h\to
0^{+}}\frac{1}{h^{r_{1}}}\frac{\partial\mu}{\partial\alpha}(0,0)\neq 0,$
(4.16)

then the equation (4.13) has a solution $\alpha^{*}$ for $h$ sufficiently
small, and this solution can be written as

$\displaystyle\alpha^{*}(h)=\eta(h)h^{r-r_{1}}$ (4.17)

for some smooth function $\eta(\cdot)$, with $\lim_{h\to 0^{+}}\eta(h)\neq 0$.

###### Proof.

The proof can be done by the implicit function theorem, following similar
approach as in [BIT]. Define a new variable $\eta=h^{-(r-r_{1})}\alpha.$ Then,

$\tilde{\mu}(\eta,h)=h^{-r}\mu(h^{(r-r_{1})}\eta,h).$

By the assumptions, $\tilde{\mu}$ is analytic in a neighborhood of $(0,0)$.
Moreover,

$\frac{\partial\tilde{\mu}}{\partial\eta}(0,0)\neq 0,$

by (4.16). Clearly, $\tilde{\mu}=0$ implies $\mu=0$. By implicit function
theorem, $\tilde{\mu}=0$ implicitly defines a smooth function $\eta=\eta(h)$
for $h\in J$, where $J$ is a small interval containing $0$. By Assumption 4.1,
it is not hard to see (depending on $c_{0,r}$ and $c_{1,r_{1}}$ in (4.15))
$\eta(0)\neq 0.$ The claims then follow if we define
$\alpha^{*}=h^{r-r_{1}}\eta$. ∎

With the relation $\alpha^{*}=\eta(h)h^{r-r_{1}}$, the convergence order of
the method given in Algorithm 1 is in fact can be the same as that for
$\alpha=0$. This means that the order of convergence can be preserved. We in
fact have the following claim.

######  Corollary 4.1.

Suppose that for all $\alpha\neq 0$ near $0$, the order of PRK is $r_{1}-1$
and that the order of PRK scheme for $\alpha=0$ is $r-1$. If the conditions in
Theorem 4.3 hold, then the order of convergence for the method given in
Algorithm 1 is equal to $r-1$.

###### Proof.

By the assumptions,

$\|p_{1}(\alpha,h)-p(h)\|=\eta_{1}(\alpha,h)h^{r_{1}},~{}~{}\|q_{1}(\alpha,h)-q(h)\|=\eta_{2}(\alpha,h)h^{r_{1}},$

where $(p(h),q(h))$ are the exact solutions for the Hamiltonian system with
initial data $(p_{0},q_{0})$. Here, $\eta_{1}(\alpha,h)$ and
$\eta_{2}(\alpha,h)$ are bounded smooth functions. Then, it is easy to see
$|\eta_{i}(\alpha,h)-\eta_{i}(0,h)|\leq C|\alpha|$ for $i=1,2$. Hence, with
$\alpha=\alpha^{*}$

$\|p_{1}(\alpha,h)-p_{1}(0,h)\|\leq Ch^{r_{1}}h^{r-r_{1}}=Ch^{r}.$

The same is true for $q_{1}$. By the assumption on the accuracy for
$\alpha=0$, the claim follows. ∎

As a direct application of the Theorem 4.3 and Corrollary 4.1, we in the next
try to find more explicit conditions corresponding to the condition in Theorem
4.3 to hold (namely $\lim_{h\to 0^{+}}h^{-r_{1}}\partial_{\alpha}\mu\neq 0$)
and determine $r_{1}$ for the modified Shake and Rattle algorithm described in
section 4.1.1. We first of all assume the following.

######  Assumption 4.2.

The matrix $G(q)$ is of rank $m$ for all $q$, and $H_{pp}(p,q)$ is invertible
for all $(p,q)$.

Clearly, the first is reasonable as we assume the $m$ holonomic constraints
are independent. Below, we perform the discussion with Assumption 4.2.
Motivated by the condition in Theorem 4.3, we define the operator

$\displaystyle\mathcal{D}_{s}f:=\lim_{h\to
0^{+}}\Big{[}h^{-s}\partial_{\alpha}f(\alpha,h)|_{\alpha=0}\Big{]}.$ (4.18)

Since $\lim_{h\to 0^{+}}p_{1}(0,h)=p_{0}$ and $\lim_{h\to
0^{+}}q_{1}(0,h)=q_{0}$, one has

$\displaystyle\mathcal{D}_{s}\mu=(\mathcal{D}_{s}p_{1})\cdot
H_{p}(p_{0},q_{0})+(\mathcal{D}_{s}q_{1})\cdot H_{q}(p_{0},q_{0}).$ (4.19)

Now, we aim to find the minimum $r$ suh that $\mathcal{D}_{r}\mu\neq 0$. First
consider $s=1$. By (4.8) and $P(0,h)\to p_{0}$, one has

$\mathcal{D}_{1}P=-H_{q}(p_{0},q_{0})-\frac{1}{2}(\mathcal{D}_{0}P)\cdot
H_{pq}(p_{0},q_{0})-G^{T}(q_{0})\cdot\mathcal{D}_{0}\lambda_{1}.$

Clearly, $\mathcal{D}_{0}P=0$. Hence,
$\mathcal{D}_{1}P=-H_{q}(p_{0},q_{0})-G^{T}(q_{0})\cdot\mathcal{D}_{0}\lambda_{1}.$
Consequently,

$\begin{split}\mathcal{D}_{1}q_{1}&=H_{p}(p_{0},q_{0})-H_{p}(p_{0},q_{0})+\frac{1}{2}[\mathcal{D}_{0}P\cdot
H_{pp}(p_{0},q_{0})+\mathcal{D}_{0}PH_{pp}(p_{0},q_{0})+\mathcal{D}_{0}q_{1}H_{qp}(p_{0},q_{0})]\\\
&=0.\end{split}$

Similarly,

$\begin{split}\mathcal{D}_{1}p_{1}&=[-H_{q}(p_{0},q_{0}+H_{q}(p_{0},q_{0}))]-G^{T}(q_{0})\mathcal{D}_{0}(\lambda_{1}+\lambda_{2})-\mathcal{D}_{0}q_{1}\cdot\nabla_{qq}g(q_{0})\cdot\lim_{h\to
0^{+}}\lambda_{2}(0,h)\\\
&=-G^{T}(q_{0})\mathcal{D}_{0}(\lambda_{1}+\lambda_{2}).\end{split}$

By the constraint of $p_{1}$ in (4.9), $0=\mathcal{D}_{0}q_{1}\otimes
H_{p}(p_{0},q_{0}):\nabla^{2}g(q_{0})+G(q_{0})\cdot[\mathcal{D}_{1}p_{1}\cdot
H_{pp}(p_{0},q_{0})]=G(q_{0})\cdot[\mathcal{D}_{1}p_{1}\cdot
H_{pp}(p_{0},q_{0})].$ Since the $m$ constraints are independent and $H_{pp}$
is invertible, $GH_{pp}G^{T}$ is invertible and thus
$\mathcal{D}_{0}(\lambda_{1}+\lambda_{2})=0.$ Hence, $\mathcal{D}_{1}p_{1}=0.$

Now, we move to $s=2$. By the $q_{1}$ equation in (4.8),

$\begin{split}\mathcal{D}_{2}q_{1}&=\lim_{h\to
0^{+}}h^{-1}[H_{p}(P,q_{0})-H_{p}(P,q_{1})]+\mathcal{D}_{1}P\cdot
H_{pp}(p_{0},q_{0})+\frac{1}{2}\mathcal{D}_{1}q_{1}\cdot
H_{qp}(p_{0},q_{0})\\\ &=-H_{p}(p_{0},q_{0})\cdot
H_{qp}(p_{0},q_{0})-[H_{q}(p_{0},q_{0})+G^{T}(q_{0})\cdot\mathcal{D}_{0}\lambda_{1}]\cdot
H_{pp}(p_{0},q_{0})\end{split}$

By the constraint in (4.8), one has $G(q_{0})\mathcal{D}_{2}q_{1}=0$. Using
this, one can determine $\mathcal{D}_{0}\lambda_{1}$.

Similarly, we have

$\begin{split}\mathcal{D}_{2}p_{1}&=\lim
h^{-1}(H_{q}(P,q_{1})-H_{q}(P,q_{0}))-\mathcal{D}_{1}P\cdot
H_{pq}(p_{0},q_{0})-\frac{1}{2}\mathcal{D}_{1}q_{1}\cdot
H_{qq}-G^{T}(q_{0})\mathcal{D}_{1}[\lambda_{1}+\lambda_{2}]-\\\
&\mathcal{D}_{1}q_{1}\cdot\nabla_{2}g(q_{0})\cdot\lim\lambda_{2}\\\
&=H_{p}(p_{0},q_{0})\cdot
H_{qq}(p_{0},q_{0})+[H_{q}(p_{0},q_{0})+G^{T}(q_{0})\mathcal{D}_{0}\lambda_{1}]\cdot
H_{pq}(p_{0},q_{0})-G^{T}(q_{0})\mathcal{D}_{1}[\lambda_{1}+\lambda_{2}].\end{split}$

Moreover, By the constraint of $p_{1}$ in (4.9), we have

$0=H_{p}(p_{0},q_{0})\otimes\mathcal{D}_{2}q_{1}:\nabla_{qq}g(q_{0})+(\mathcal{D}_{2}p_{1}\cdot
H_{pp}(p_{0},q_{0})+\mathcal{D}_{2}q_{1}\cdot
H_{qp}(p_{0},q_{0}))\cdot\nabla_{q}g(q_{0})$

With this, one can determine $\mathcal{D}_{1}(\lambda_{1}+\lambda_{2})$.
Hence, we can find that

$\displaystyle\mathcal{D}_{2}\mu=H_{p}\otimes H_{p}:H_{qq}-H_{q}\otimes
H_{q}:H_{pp}+R,$ (4.20)

where

$R:=[G^{T}\mathcal{D}_{0}\lambda_{1}]\otimes
H_{p}:H_{pq}-H_{p}G^{T}\mathcal{D}_{1}[\lambda_{1}+\lambda_{2}]-[G^{T}\mathcal{D}_{0}\lambda_{1}]\otimes
H_{q}:H_{pp}.$

Every term in $R$ involves $G$, while the first two terms in
$\mathcal{D}_{2}\mu$ do not depend on $G$. Hence, $\mathcal{D}_{2}\mu\neq 0$
for all most all $(p,q)$. It is possible that at some special points, this can
be zero.

######  Remark 4.2.

For a more straightforward derivation, one may take derivatives on $\alpha$ in
(4.8)-(4.9) and then expand the quantities like
$\partial_{\alpha}P,\partial_{\alpha}p_{1}$ etc in terms of powers of $h$.
Then compare the coefficients, one can derive the same formulas.

To summarize, we expect that for almost all $(p_{0},q_{0})$,

$\lim_{h\to 0^{+}}\frac{1}{h^{2}}\frac{\partial\mu(0,h)}{\partial\alpha}\neq
0.$

This coincides with the fact the method for $\alpha\neq 0$ is first order.
Hence, the method should be second order by Corollary 4.1. We know the order
of convergence for $\alpha=0$ is two or $r=3$. Hence, mostly, we will have

$\alpha^{*}=O(h).$

It could be possible that at some special points, $\lim_{h\to
0^{+}}\frac{1}{h^{2}}\frac{\partial\mu(0,h)}{\partial\alpha}=0$. For these
points, there could still be solutions for $\alpha^{*}$ but one may have
larger $\alpha^{*}\gg h$, as we need to look into higher orders.

##  5 Numerical Experiments

In this section, we present numerical examples to verify our theoretical
findings. We verify that the method constructed in this paper can accurately
conserve the quadratic invariants and constraint manifolds of the system, and
at the same time, the appropriate parameters $\alpha^{*}$ can be selected to
make the system energy conservation at each step. As far as we know, this is
the first time that there is a numerical method for a Hamiltonian constrained
systems that can hold all three quantities at once. We will also verify the
order of convergence of the parameterized method and prove that the parameter
perturbation will not decrease the order of convergence of the numerical
solutions. All the experiments in this paper are implemented using Julia on
Version 1.4.2 and Packages: NLsolve v4.4.0.

###  5.1 Spherical pendulum

$\begin{array}[]{cc}\includegraphics[scale={0.35}]{WechatIMG3051.png}\includegraphics[scale={0.35}]{WechatIMG3052.png}\end{array}$

Figure 1: Numerical solutions.

$\begin{array}[]{cc}\includegraphics[scale={0.4}]{L-rattle.png}\quad\quad\includegraphics[scale={0.4}]{L-prk.png}\end{array}$

Figure 2: The errors in the third component of the angular momentum for
standard Rattle and Lobatto IIIA-IIB methods and the $\alpha$-Rattle and
$\alpha$-PRK III methods.

(a) Energy errors for standard Rattle method.

(b) Energy errors for standard Lobatto IIIA-IIIB method.

(c) $\alpha^{*}$ for $\alpha$-Rattle method.

(d) Energy errors for $\alpha$-Rattle method.

(e) $\alpha^{*}$ for $\alpha$-PRK III method.

(f) Energy errors for $\alpha$-PRK III method.

Figure 3:  $\alpha^{*}$ and energy errors with $h=0.1$.

Consider the spherical pendulum equation with
$H(p,q)=\frac{1}{2}(p_{1}^{2}+p_{2}^{2}+p_{3}^{2})+q_{3}$ and
$g(q)=(q_{1}^{2}+q_{2}^{2}+q_{3}^{2})-1$, so that

$\begin{split}H_{p}(p,q)=\begin{bmatrix}p_{1}\\\ p_{2}\\\ p_{3}\\\
\end{bmatrix},~{}H_{q}(p,q)=\begin{bmatrix}0\\\ 0\\\ 1\\\
\end{bmatrix},~{}G^{T}(q)=2\begin{bmatrix}q_{1}\\\ q_{2}\\\ q_{3}\\\
\end{bmatrix}.\end{split}$ (5.1)

For this example, the Lagrangian

$\displaystyle L(Q,\dot{Q})=\frac{1}{2}|\dot{Q}|^{2}-Q_{3},$

is invariant under rotations about $z$ axis. In particular, taking $D_{1}$ to
be a $2\times 2$ anti-symmetric matrix and setting

$\displaystyle D:=\left[\begin{array}[]{cc}D_{1}&0\\\ 0&0\end{array}\right].$
(5.4)

Then, using $D^{T}=-D$, and $(e^{sD}v)_{3}=v_{3}$, it can be verified easily
that

$\displaystyle
g(e^{sD}q)=q^{T}e^{sD^{T}}e^{sD}q-1=q^{T}q-1=0,~{}~{}~{}L(e^{sD}\dot{q},e^{sD}q)=\frac{1}{2}|e^{sD}\dot{q}|^{2}-(e^{sD}q)_{3}=L(\dot{q},q).$

Hence, the quantity $I(p,q)=q^{T}Dp$ is invariant under this Hamiltonian
dynamics. In particular, this implies that the third component of the angular
momentum $L(p,q)=q\times p$ is conservative. Physically, the forces (force due
to the constraint which points along the radius of the sphere, the gravity),
$q$ and $z$-axis are coplanar, so $z$-component of the torque is zero. Thus,
the $z$-component of the angular moment is conserved. We will verify through
numerical experiments in this example that the parameterized method we
constructed can preserve both the Hamilton function and $L_{3}$.

Table 1: The errors $e_{p}$, $e_{q}$ and the orders at $T=0.5$ for the
$\alpha$-Rattle method for Example 1

$h$ | $e_{p}$ | order | $e_{q}$ | order
---|---|---|---|---
0.25 | 3.3643e-4 | - | 3.5220e-4 | -
0.125 | 8.6813e-5 | 1.9543 | 8.9671e-5 | 1.9737
0.0625 | 2.1895e-5 | 1.9872 | 2.2535e-5 | 1.9924
0.03125 | 5.4863e-6 | 1.9967 | 5.6416e-6 | 1.9980
0.015625 | 1.7323e-6 | 1.9992 | 1.4108e-6 | 1.9995

Table 2: The errors $e_{p}$, $e_{q}$ and the orders at $T=1$ for the
$\alpha$-Rattle method for Example 1

$h$ | $e_{p}$ | order | $e_{q}$ | order
---|---|---|---|---
0.1 | 1.1216e-4 | - | 1.0774e-4 | -
0.05 | 4.5419e-5 | 1.3042 | 4.8323e-5 | 1.1569
0.025 | 1.5315e-5 | 1.5683 | 2.0176e-5 | 1.2600
0.0125 | 1.6383e-6 | 3.2247 | 1.5628e-6 | 3.6905
0.00625 | 5.7032e-7 | 1.5226 | 5.9045e-7 | 1.4043

Table 3: The errors $e_{p}$, $e_{q}$ and the orders at $T=1$ for the
$\alpha$-PRK III method for Example 1

$h$ | $e_{p}$ | order | $e_{q}$ | order
---|---|---|---|---
0.25 | 4.0025e-7 | - | 4.7611e-7 | -
0.125 | 2.5089e-8 | 3.9957 | 2.9843e-8 | 3.9958
0.0625 | 1.5692e-9 | 3.9989 | 1.8665e-9 | 3.9989
0.03125 | 9.8091e-11 | 3.9998 | 1.1666e-10 | 3.9998
0.015625 | 5.9634e-12 | 4.0399 | 6.9571e-12 | 4.0677

We take the initial value $q_{0}=(0,\sin(0.1),-\cos(0.1)),p_{0}=(0.06,0,0)$
for the simulations. In Fig. 2, we present the third component of the angular
momentum $L(p,q)$. From which we can see that both the standard Rattle and
Lobatto IIIA-IIIB methods and the parameterized methods constructed in this
paper can preserve this quadratic invariants to machine precision, that is up
to $10^{-16}$.

In Fig. 3, we show the results for the chosen parameters $\alpha$ and computed
Hamiltonians for the $\alpha$-Rattle in (4.8)-(4.9) and $\alpha$-PRK III
method in (4.11). Clearly, the standard symplectic methods (i.e. Rattle and
Lobatto IIIA-IIIB methods) have energy error range from $10^{-9}$ to
$10^{-6}$, while the parameterized symplectic methods can make the energy
conservation in the sense of that the energy error reach the machine accuracy.
At each step, we can find parameters $\alpha^{*}=\alpha(h,y_{n})$ that achieve
energy conservation. We can see from $(c)$ in Fig. 3 that the parameter is
small in most cases for $\alpha$-Rattle method, but it may take a slightly
larger value at a very few points, and such an exceptional points may lower
the order of the parameterized methods.

In the Table 3 and 3 we give the specific order of convergence for the
$\alpha$-Rattle method. We can see that the order of Table 1 is completely
consistent with the standard algorithm of second order when computed on
$[0,0.5]$, while the order of Table 2 is slightly reduced when we computed on
$[0,1]$. We found that there is some $\alpha$ that is big as shown in Figure 3
(c) between $t=0.5$ and $t=1$, and this seems to reduce the computational
accuracy as commented above. We are not clear with the deeper reason for this.
Maybe, the conditions in Theorem 4.3 might be violated at some point, or maybe
the nonlinear system is not solved with high accuracy as commented below in
Remark 5.1. For the 3 stage $\alpha$-PRK method, we can see from $(e)$ in
Figure 3 that the parameter is uniform small about $O(10^{-5})$. As shown in
Table 3, the convergence order of $\alpha$-PRK III method is preserved exactly
as the standard 3 stage PRK method. That is forth order.

######  Remark 5.1.

In the specific implementation of the parameterized methods presented in this
article, it is very important to accurately solve the corresponding nonlinear
system of equations (like (4.6)-(4.8)). If the accuracy of the numerical
solutions is not enough (for example, the constraints are not satisfied to
enough accuracy), the accuracy of the energy computation might be reduced to
like from $10^{-8}$ to $10^{-12}$, and the searched parameters $\alpha^{*}$
might not be very good. One may try to improve the accuracy of nonlinear
equations by giving a good initial value during the iteration. For example,
use the numerical results of the standard PRK methods as the initial value in
the iteration of the nonlinear equations for the corresponding parameterized
methods. For future work, one may introduce some correction field or augmented
variables to make the nonlinear system more stable to solve.

###  5.2 Satellites system

$\begin{array}[]{cc}\includegraphics[scale={0.45}]{ex2-p.png}\includegraphics[scale={0.45}]{ex2-q.png}\end{array}$

Figure 4: Numerical solutions obtained by $\alpha$-Rattle method with $h=0.1$.

We discuss a closed-loop rotating triangular tethered satellites system,
including three satellites (considered as mass-points) of masses
$m_{i},i=1,2,3,$ joined by inextensible, tight, and massless tethers, of
lengths $L_{i},i=1,2,3$, respectively. For sake of simplicity, we assume unit
masses and lengths, and normalize the gravity constant. Consequently, if
$q_{i}:=(x_{i},y_{i},z_{i})^{T}\in\mathbb{R}^{3},i=1,2,3$, are the positions
of the three satellites, the constraints $0=g(q)$ are given by

$\begin{split}\begin{cases}&g_{1}(q)=(q_{1}-q_{2})^{T}(q_{1}-q_{2})-1=0,\\\
&g_{2}(q)=(q_{2}-q_{3})^{T}(q_{2}-q_{3})-1=0,\\\
&g_{3}(q)=(q_{3}-q_{1})^{T}(q_{3}-q_{1})-1=0.\\\ \end{cases}\end{split}$ (5.5)

The Hamiltonian is specified by
$H(p,q)=\sum_{i=1}^{3}\left(\frac{1}{2}p_{i}^{T}p_{i}-\frac{1}{\sqrt{q_{i}^{T}q_{i}}}\right).$

Hence, we have that

$\begin{split}&H_{p}=p=\begin{bmatrix}p_{1}\\\ p_{2}\\\ p_{3}\\\
\end{bmatrix},~{}H_{q}=f(q):=\begin{bmatrix}(q_{1}^{T}q_{1})^{-\frac{3}{2}}q_{1}\\\
(q_{2}^{T}q_{2})^{-\frac{3}{2}}q_{2}\\\
(q_{3}^{T}q_{3})^{-\frac{3}{2}}q_{3}\\\
\end{bmatrix},~{}G^{T}(q)=2\begin{bmatrix}q_{1}-q_{2}&0&q_{1}-q_{3}\\\
q_{2}-q_{1}&q_{2}-q_{3}&0\\\ 0&q_{3}-q_{2}&q_{3}-q_{1}\\\
\end{bmatrix}.\end{split}$ (5.6)

(a) Energy errors for Rattle method with $h=0.1$.

(b) Energy errors for Lobatto IIIA-IIIB method with $h=0.2$.

(c) $\alpha^{*}$ for $\alpha$-Rattle method with $h=0.1$.

(d) Energy errors for $\alpha$-Rattle method with $h=0.1$.

(e) $\alpha^{*}$ for $\alpha$-PRK III method with $h=0.2$.

(f) Energy errors for $\alpha$-PRK III method with $h=0.2$.

Figure 5:  $\alpha^{*}$ and energy errors

For the simulation, we choose the consistent initial values given by

$\begin{split}q_{1}(0)=\begin{bmatrix}0\\\ \frac{1}{2}\\\ z_{0}\\\
\end{bmatrix},~{}q_{2}(0)=\begin{bmatrix}0\\\ -\frac{1}{2}\\\ z_{0}\\\
\end{bmatrix},~{}q_{3}(0)=\begin{bmatrix}0\\\ 0\\\ z_{0}-\frac{\sqrt{3}}{2}\\\
\end{bmatrix},~{}p_{1}(0)=p_{2}(0)=\begin{bmatrix}0\\\ 0\\\ 0\\\
\end{bmatrix},~{}p_{3}(0)=\begin{bmatrix}v_{0}\\\ 0\\\ 0\\\
\end{bmatrix}.\end{split}$

where $z_{0}=20$ and $v_{0}$ is such that the initial Hamiltonian is zero.
This provides a configuration in which the first two satellites remain
parallel to each other, moving in the planes $y=1/2$ and $y=-1/2$,
respectively, and the third one moves around the tether joining the first two,
in the plane $y=0$.

Table 4: The errors $e_{p}$, $e_{q}$and the orders at $T=1$ for the
$\alpha$-Rattle method for Example 2

$h$ | $e_{p}$ | order | $e_{q}$ | order
---|---|---|---|---
0.25 | 1.2290e-3 | - | 1.9300e-3 | -
0.125 | 3.0460e-4 | 2.0125 | 4.7835e-4 | 2.1247
0.0625 | 7.59745e-5 | 2.0033 | 1.1931e-4 | 2.0033
0.03125 | 1.8968e-5 | 2.0018 | 2.9789e-5 | 2.0018
0.015625 | 4.7270e-6 | 2.0046 | 7.4235e-6 | 2.0046

Table 5: The errors $e_{p}$, $e_{q}$and the orders at $T=1$ for the
$\alpha$-PRK III method for Example 2

$h$ | $e_{p}$ | order | $e_{q}$ | order
---|---|---|---|---
0.25 | 1.2299e-6 | - | 1.9307e-6 | -
0.125 | 7.1884e-8 | 4.0967 | 1.1284e-7 | 4.0967
0.0625 | 3.0717e-9 | 4.5485 | 4.8219e-9 | 4.5485
0.03125 | 9.7230e-11 | 4.9815 | 1.5263e-10 | 4.9815
0.015625 | 6.4283e-12 | 3.9188 | 1.0091e-11 | 3.9188

For this more complex example, our numerical results are similar to the first
example. The standard symplectic PRK methods have energy errors, but the
parameterized methods can preserves energy to machine accuracy. At each step,
we can find the parameter values that make the energy conserved. Meanwhile, as
shown in Table 5 and 5, the order of convergence of the two parameterized
methods can be completely consistent with the standard methods. By comparison,
it is found that $\alpha$-PRK III methods are often superior to the
$\alpha$-Rattle method, and the parameter values of $\alpha$-PRK III methods
are usually uniformly small, while $\alpha$-Rattle method may have a few
points with large $\alpha^{*}$ values.

## Acknowledgements

The authors would like to thank Mr. Changze Chen (Northwest University) for
his great help and support in the implementation of the numerical experiments
in this paper. The work of L. Li was partially sponsored by NSFC 11901389,
Shanghai Sailing Program 19YF1421300 and NSFC 11971314. The work of D. Wang
was partially sponsored by NSFC 11871057, 11931013 and Project for Young
Science and Technology Star of Shaanxi Province in China (2018 KJXX-070).

## References

  * [And83] H. C. Andersen.  Rattle: A ”velocity” version of the shake algorithm for molecular dynamics calculations.  J. Comput. Phys., 52(1):24–34, 1983.
  * [Arn13] Vladimir Igorevich Arnol’d.  Mathematical methods of classical mechanics, volume 60.  Springer Science & Business Media, 2013.
  * [BG] Giancarlo Benettin and Antonio Giorgilli.  On the Hamiltonian interpolation of near-to-the identity symplectic mappings with application to symplectic integration algorithms.  Journal of Statistical Physics, 74(5-6):1117–1143.
  * [BGI] Luigi Brugnano, Gianmarco Gurioli, and Felice Iavernaro.  Analysis of energy and quadratic invariant preserving (equip) methods.  Journal of Computational $\&$ Applied Mathematics, 335:51–73.
  * [BGIW] Luigi Brugnano, Gianmarco Gurioli, Felice Iavernaro, and Ewa B. Weinmüller.  Line integral solution of Hamiltonian systems with holonomic constraints.  Applied Numerical Mathematics, 127(MAY):56–77.
  * [BIT] Luigi Brugnano, Felice Iavernaro, and Donato Trigiante.  Energy- and quadratic invariants–preserving integrators based upon gauss collocation formulae.  SIAM Journal on Numerical Analysis, 50(6):2897–2916.
  * [BIT09] Luigi Brugnano, Felice Iavernaro, and Donato Trigiante.  Hamiltonian boundary value methods (energy conserving discrete line integral methods).  Mathematics, (1-2):17–37, 2009.
  * [CFM] Philippe Chartier, Erwan Faou, and Ander Murua.  An algebraic approach to invariant preserving integators: The case of quadratic and Hamiltonian invariants.  Numerische Mathematik, 103(4):575–590.
  * [CHJ20] Chuchu Chen, Jialin Hong, and Diancong Jin.  Modified averaged vector field methods preserving multiple invariants for conservative stochastic differential equations.  BIT Numer Math. https://doi.org/10.1007/s10543-020-00803-6, 2020\.
  * [Fen86] Kang Feng.  Difference schemes for Hamiltonian formalism and symplectic geometry.  Journal of Computational Mathematics, 4(3):279–289, 1986.
  * [FQ03] Kang Feng and Mengzhao Qin.  Symplectic Geometric Algorithms for Hamiltonian Systems.  Zhejiang Science and Technology Press, Hangzhou, 2003.
  * [GM] Zhong Ge and Jerrold E. Marsden.  Lie-Poisson Hamilton-Jacobi theory and Lie-Poisson integrators.  133(3):134–139.
  * [Gon96] O. Gonzalez.  Time integration and discrete hamiltonian systems.  Journal of Nonlinear Science, 6(5):449–467, 1996.
  * [HLW06] Ernst Hairer, Christian Lubich, and Gerhard Wanner.  Geometric numerical integration: structure-preserving algorithms for ordinary differential equations, volume 31.  Springer Science $\&$ Business Media, 2006.
  * [HW06] Ernst Hairer and Gerhard Wanner.  Solving Ordinary Differential Equations II stiff and differential-algebraic problems.  Sciencep, Beijing, 2006.
  * [Jay96] Laurent Jay.  Symplectic partitioned Runge–Kutta methods for constrained Hamiltonian systems.  SIAM Journal on Numerical Analysis, 33(1):368–387, 1996.
  * [Jay02] L. Jay.  Iterative solution of SPARK methods applied to DAEs.  Numerical Algorithms, 31(1-4):171–191, 2002.
  * [JWZ19] An Jing, Cao Waixiang, and Zhang Zhimin.  A class of efficient spectral methods and error analysis for nonlinear Hamiltonian systems.  Commun. Math. Sci., 18(2):395–428, 2019.
  * [MQR] R. I. Mclachlan, G. R. W. Quispel, and N. Robidoux.  Geometric integration using discrete gradients.  Philosophical Transactions of the Royal Society B Biological Sciences, 357(1754):1021–1045.
  * [MW01] Jerrold E Marsden and Matthew West.  Discrete mechanics and variational integrators.  Acta Numerica, 10:357–514, 2001.
  * [Noe71] Emmy Noether.  Invariant variation problems.  Transport Theory and Statistical Physics, 1(3):186–207, 1971.
  * [QM] G R W Quispel and D I McLaren.  A new class of energy-preserving numerical integration methods.  Journal of Physics A Mathematical $\&$ Theoretical, 41(4):045206.
  * [QW11] Mengzhao Qin and Yushun Wang.  Structure-preserving Algorithms for Partial Differential Equations.  Zhejiang Science and Technology Press, Hangzhou, 2011.
  * [Sha] Zaijiu Shang.  KAM theorem of symplectic algorithms for Hamiltonian systems.  Numerische Mathematik, 83(3):477–496.
  * [Sun93] Geng Sun.  Symplectic partitioned Runge-Kutta methods.  Journal of Computational Mathematics, 11(4):365–372, 1993.
  * [Sun95] Geng Sun.  Construction of high order symplectic PRK methods.  Journal of Computational Mathematics, 13(1):40–50, 1995.
  * [Tan] Yifa Tang.  Formal energy of a symplectic scheme for hamiltonian systems and its applications (i).  Computers $\&$ Mathematics with Applications, 27(7):31–39.
  * [WDQW] Yi Wei, Zichen Deng, L. I. Qingjun, and Bo Wang.  Projected Runge-Kutta methods for constrained Hamiltonian systems.  Applied Mathematics and Mechanics, 037(8):1077–1094.
  * [WOL17] T Wenger, S Oberblöbaum, and S Leyendecker.  Construction and analysis of higher order variational integrators for dynamical systems with holonomic constraints.  Advances in Computational Mathematics, 43(5):1–33, 2017.
  * [WXL] Dongling Wang, Aiguo Xiao, and Xueyang Li.  Parametric symplectic partitioned Runge-Kutta methods with energy-preserving properties for Hamiltonian systems.  Computer Physics Communications, 184(2):303–310.