Line Integral Solution of Hamiltonian Systems with Holonomic Constraints

TL;DR

This paper introduces a second-order energy-conserving numerical method for Hamiltonian systems with holonomic constraints, utilizing a line integral framework to ensure accuracy and conservation properties.

Contribution

The paper presents a novel line integral-based approach for energy-conserving approximations in constrained Hamiltonian systems, advancing numerical methods in this area.

Findings

The method conserves energy accurately in numerical experiments.

Numerical results validate the theoretical properties of the proposed approach.

The approach is applicable to a range of holonomic constrained systems.

Abstract

In this paper, we propose a second-order energy-conserving approximation procedure for Hamiltonian systems with holonomic constraints. The derivation of the procedure relies on the use of the so-called line integral framework. We provide numerical experiments to illustrate theoretical findings.

Equations

$$H(q,p)=\frac{1}{2}{p}^{\top }{M}^{-1}p+U(q),q,p\in {\mathbb{R}}^m,$$

$$g(q)=0\in {\mathbb{R}}^{\nu },$$

$$\hat{H}(q,p,\lambda )=H(q,p)+{\lambda }^{\top }g(q),$$

$$\dot{q}=M^{-1}p,\dot{p}=-\nabla U(q)-\nabla g(q)\lambda ,g(q)=0,t\in [0,T],$$

$$q(0)=q_0,p(0)=p_0,$$

$$g(q_0)=0,\nabla g(q_0{)}^{\top }{M}^{-1}{p}_0=0.$$

$$\mathcal{M}=\left\{q\in {\mathbb{R}}^m:g(q)=0\right\},$$

$$\dot{g}(q)=\nabla g(q{)}^{\top }\dot{q}=\nabla g(q{)}^{\top }{M}^{-1}p=0,$$

$$q_1\approx q(h),p_1\approx p(h),$$

$$H(q_1,p_1)=H(q_0,p_0),g(q_1)=g(q_0)=0,$$

$$\nabla g(q_1{)}^{\top }{M}^{-1}{p}_1=O(h^2).$$

$$\ddot{g}(q)={\nabla }^{2}g(q)(M^{-1}p,M^{-1}p)-\nabla g(q{)}^{\top }{M}^{-1}\left[\nabla U(q)+\nabla g(q)\lambda \right].$$

$$\left[\nabla g(q{)}^{\top }{M}^{-1}\nabla g(q)\right]\lambda ={\nabla }^{2}g(q)(M^{-1}p,M^{-1}p)-\nabla g(q{)}^{\top }{M}^{-1}\nabla U(q).$$

$$\nabla g(q{)}^{\top }{M}^{-1}\nabla g(q)$$

$$\lambda$$

$$g(q(h))-\underset{=0}{\underbrace{g(q(0))}}={\int }_0^h\nabla g(q(t){)}^{\top }\dot{q}(t)\mathrm{d}t=0.$$

$$u(0)=q_0,u(h)=:q_1\approx q(h),$$

$$g(q_1)-g(q_0)\equiv g(u(h))-g(u(0))={\int }_0^h\nabla g(u(t){)}^{\top }\dot{u}(t)\mathrm{d}t=0,$$

$$\hat{H}(q_1,p_1,\lambda )=\hat{H}(q_0,p_0,\lambda ),g(q_1)=g(q_0)=0,$$

$$\dot{u}=M^{-1}v,\dot{v}=-\nabla U(u)-\nabla g(u){\lambda }_0,t\in [0,h],$$

$$u(0)=q_0,v(0)=p_0,$$

$$q_1:=u(h),p_1:=v(h),$$

$$q_n=q(nh)+O(h^2),p_n=p(nh)+O(h^2),g(q_n)=0,\nabla g(q_n{)}^{\top }{M}^{-1}{p}_n=O(h^2).$$

$${\lambda }_n=\lambda (q(nh),p(nh))+O(h),$$

$$P_j\in {\Pi }_j,{\int }_0^1{P}_i(c)P_j(c)\mathrm{d}c={\delta }_{ij},\forall i,j=0,1,\dots ,$$

$$M^{-1}v(ch)$$

$$\nabla g(u(ch))$$

$${\gamma }_j(v)$$

$${\rho }_j(u)$$

$$\dot{u}(ch)=\sum _{j\ge 0}{P}_j(c){\gamma }_j(v),\dot{v}(ch)=-\sum _{j\ge 0}{P}_j(c)[{\psi }_j(u)+{\rho }_j(u){\lambda }_0],c\in [0,1].$$

$$u(ch)=q_0+h\sum _{j\ge 0}{\int }_0^c{P}_j(x)\mathrm{d}x{\gamma }_j(v),v(ch)=p_0-h\sum _{j\ge 0}{\int }_0^c{P}_j(x)\mathrm{d}x[{\psi }_j(u)+{\rho }_j(u){\lambda }_0],c\in [0,1].$$

Full text

Keywords: constrained Hamiltonian systems; holonomic constraints; energy-conserving methods; line integral methods; Hamiltonian Boundary Value Methods; HBVMs.

MSC: 65P10, 65L80, 65L06.

1 Introduction

We consider the numerical approximation of a constrained Hamiltonian dynamics, described by the separable Hamiltonian

[TABLE]

where $M$ is a symmetric and positive-definite matrix. The problem is completed by $\nu$ holonomic constraints,

[TABLE]

where we assume that $\nu<m$ holds. Moreover, we also assume that all points are regular for the constraints, i.e., $\nabla g(q)\in\mathbb{R}^{m\times\nu}$ has full column rank or, equivalently, $\nabla g(q)^{\top}M^{-1}\nabla g(q)$ is nonsingular. For simplicity, both $U$ and $g$ are assumed to be analytic.

It is well-known that the problem defined by (1)–(2) can be cast in Hamiltonian form by defining the augmented Hamiltonian

[TABLE]

where $\lambda$ is the vector of Lagrange multipliers. The resulting constrained Hamiltonian system reads:

[TABLE]

and is subject to consistent initial conditions,

[TABLE]

such that

[TABLE]

Note that the condition $g(q_{0})=0$ ensures that $q_{0}$ belongs to the manifold

[TABLE]

as required by the constraints, whereas the condition $\nabla g(q_{0})^{\top}M^{-1}p_{0}$ means that the motion initially stays on the tangent space to ${\cal M}$ at $q_{0}$. On a continuous level, this condition is satisfied by all points on the solution trajectory, since, in order for the constraints to be conserved,

[TABLE]

holds. These latter constraints are sometimes referred to as hidden constraints.

We stress that the condition (8) can be conveniently relaxed for the numerical approximation. There, we only ask for $\nabla g(q)^{\top}M^{-1}p$ to be suitably small along the numerical solution. Consequently, when solving the problem on the interval $[0,h]$, we require that the approximations,

[TABLE]

satisfy the conservation of both, the Hamiltonian and the constraints,

[TABLE]

and that the hidden constraints are relaxed to

[TABLE]

We recall that a formal expression for the vector $\lambda$ is obtained by an additional differentiating of (8), i.e.,

[TABLE]

Imposing the vanishing of this derivative yields

[TABLE]

Consequently, the following result follows.

Theorem 1

The vector $\lambda$ exists and is uniquely determined, provided that the matrix

[TABLE]

is nonsingular. In fact, in such a case, from (13), we obtain

[TABLE]

Note that, for later use, we have introduced the notation $\lambda(q(t),p(t))$ in place of $\lambda(t)$, to explicitly underline the dependence of the Lagrange multiplier on the state variables $q$ and $p$ at time $t$.

Remark 1

We observe that an additional differentiation of (13) provides a differential equation for the Lagrange multipliers, which can be solved together with the original problem. However, this procedure is cumbersome in general, since it requires the evaluation of higher order tensors.

Numerical solution of Hamiltonian problems with holonomic constraints has been for a long time in the focus of interest. Many different approaches have been proposed such as the basic Shake-Rattle method [38, 3], which has been shown to be symplectic [31], higher order methods obtained via symplectic PRK methods [24], composition methods [34, 35], symmetric LMFs [23]. Further methods are based on discrete derivatives [33], local parametrizations of the manifold containing the solution [4], or on projection techniques [36, 39]. See also [29, 37, 25, 32] and the monographs [5, 30, 27, 28].

In this paper we pursue a different approach, utilizing the so-called line integral, which has already been used when deriving the energy-conserving Runge-Kutta methods, for unconstrained Hamiltonian systems, cf. Hamiltonian Boundary Value Methods (HBVMs) [11, 12, 13, 16, 17] and the recent monograph [10]. Such methods have also been applied in a number of applications [6, 9, 14, 15, 2, 8, 1], and here are used to cope with the constrained problem (1)–(2). Roughly speaking, the conservation of the invariant will be guaranteed by requiring that a suitable line integral vanishes. This line integral represents a discrete-time version of (8). In fact, if we fix a stepsize $h>0$, then the conservation of the constraints (2) at $h$, starting from the point $q_{0}$ defined in (5), can be recast into

[TABLE]

For the continuous solution, this integral vanishes since the integrand is identically zero due to (4) and (8). However, we can relax this requirement in the context of a numerical method describing a discrete-time dynamics. In such a case, the conservation properties have to be satisfied only on a set of discrete times which are multiples of the stepsize $h$. Consequently, we consider a local approximation to $q(t)$, say $u(t)$, such that

[TABLE]

and

[TABLE]

without requiring the integrand to be identically zero. This, in turn, enables a proper choice of the vector of the multipliers $\lambda$. As a result, we eventually obtain suitably modified HBVMs which enable to conserve both, the Hamiltonian and the constraints. We stress that the available efficient implementation of the original methods (see, e.g., [13, 7, 10]), which proved to be reliable and robust in the numerical solution of the unconstrained Hamiltonian problems, can now be adapted for dealing with the holonomic constraints.

The paper is organized as follows. In Section 2, we provide the framework for devising the method via a suitable choice of the vector $\lambda$ of the Lagrange multipliers, which we approximate by a piecewise-constant function. In Section 3 further simplification towards numerical procedure is discussed. Then, in Section 4, we present a fully discrete method, resulting in a suitable modification of the original HBVMs. In Section 5, numerical experiments are shown to illustrate how the method works for a number of constrained Hamiltonian problems. Section 6 contains the conclusions and possible future investigations.

2 Piecewise-constant approximation of $\lambda$

In this section, we show that we can approximate the solution of problem (4) on the interval $[0,h]$, $h=T/N$, by looking for a constant vector $\lambda\in\mathbb{R}^{\nu}$ such that (10) is satisfied. This is equivalent to require

[TABLE]

where the constant parameter $\lambda$ is chosen in such a way that the constraints $g(q_{1})=0$ hold. We will show that this procedure provides us with a second order approximation of the original problem, which becomes exact when the true multiplier is constant. Consequently, we approximate the problem (4)–(5), by the local problem

[TABLE]

subject to the initial conditions, cf. (5),

[TABLE]

satisfying (6). By setting

[TABLE]

the constant parameter ${\lambda_{0}}$ is chosen to guarantee the conservation of the Hamiltonian and the constraints, i.e., (10). Starting from (18), the procedure is then repeated on $[h,2h]$ and the following intervals. The convergence result is now formulated in the following theorem.

Theorem 2

For all sufficiently small stepsizes $h>0$, the above procedure defines a sequence of approximations $(q_{n},p_{n})$ such that, for all $n=1,2,\ldots:$

[TABLE]

Moreover, $(q_{n+1},p_{n+1})$ is obtained from $(q_{n},p_{n})$ using a constant vector $\lambda_{n}$ such that

[TABLE]

where $\lambda(q,p)$ is defined in (14) and, consequently,  $H(q_{n+1},p_{n+1})=H(q_{n},p_{n}).$

The aim of this section is to show (19)–(20). Let us first consider the orthonormal basis on $[0,1]$ given by the shifted and scaled Legendre polynomials $\{P_{j}\}$,

[TABLE]

along with the expansions,

[TABLE]

with

[TABLE]

Consequently, following the approach defined in [16], the differential equations in (16) can be rewritten as

[TABLE]

Moreover, using the initial conditions (5), we formally obtain

[TABLE]

The following result is now cited from [16, Lemma 1].

Lemma 1

Let $G:[0,h]\rightarrow V$, with $V$ a vector space, admit a Taylor expansion at [math]. Then

[TABLE]

As a straightforward consequence, one has the following result.

Corollary 1

All coefficients specified in (23) are $O(h^{j})$.

Concerning the conservation properties of the approximations, the following result holds.

Theorem 3

For all ${\lambda_{0}}\in\mathbb{R}^{\nu}$, the solution of (16)–(18) satisfies

[TABLE]

Proof For any given ${\lambda_{0}}\in\mathbb{R}^{\nu}$, it follows from (3), (16), and (18),

[TABLE]

[TABLE]

As observed above, the conservation of the Hamiltonian (1) is guaranteed, once the constraints are satisfied, i.e., $g(q_{1})=0$. We now apply a line integral technique to determine the vector ${\lambda_{0}}$ and formulate the following result describing the very first step of the approximation procedure.

Theorem 4

Let us consider the problem (16)–(17) and assume that $(q_{0},p_{0})$ is given such that,

• •

$\nabla g(q_{0})^{\top}M^{-1}\nabla g(q_{0})\in\mathbb{R}^{\nu\times\nu}$* is nonsingular;*

• •

$g(q_{0})=0$;

• •

$\nabla g(q_{0})^{\top}M^{-1}p_{0}=0$.

Then, for all sufficiently small $h>0$, $\exists!{\lambda_{0}}\in\mathbb{R}^{\nu}$ such that the approximations in (18) satisfy

• •

$g(q_{1})=0$* and, therefore, $H(q_{1},p_{1})=H(q_{0},p_{0})$;*

• •

$\lambda_{0}=\lambda(q_{0},p_{0})+O(h)$;

• •

$q_{1}-q(h)=O(h^{2})$, $p_{1}-p(h)=O(h^{2})$;

• •

$\nabla g(q_{1})^{\top}M^{-1}p_{1}=O(h^{2})$.

Remark 2

Clearly, Theorem 4 is the discrete counterpart of Theorem 1.

Before showing Theorem 4, we have to state the following preliminary results.

Lemma 2

Let us consider the polynomial basis (21). Then, we have

[TABLE]

where $\left(X_{s}\right)_{j+1,i+1}$ is the $(j+1,i+1)$ entry of the matrix

[TABLE]

Proof Since the integrand on the left-hand side in (26) is a polynomial of degree at most $2s-1$, the integral can be computed exactly via the Gauss-Legendre formula of order $2s$. Let $c_{1},\dots,c_{s}$ be the zeros of $P_{s}$ and $b_{1},\dots,b_{s}$ be the corresponding weights. Then, introducing the matrices

[TABLE]

and setting $e_{i}\in\mathbb{R}^{s}$, the $i$-th unit vector, we have

[TABLE]

The result follows by observing that, due to the properties of Legendre polynomials [10, Section 1.4.3],

[TABLE]

follows, where $X_{s}$ is the matrix defined in (27), and $I_{s}\in\mathbb{R}^{s\times s}$ is the identity matrix. This yields

[TABLE]

We also need the following expansions.

Lemma 3

From (16) and (23), we conclude

[TABLE]

Proof The first equality follows from Lemma 1 and Corollary 1,

[TABLE]

The second statement follows by transposition and multiplication from the right by $M^{-1}p_{0}$.  $\Box{~{}}$

We now show the results formulated in Theorem 4.

Proof (of Theorem 4).   Let us assume that $g(q_{0})=0$ holds. Then, it follows from (16)–(18),

[TABLE]

Due to (21),

[TABLE]

and according to (26)–(27), we conclude

[TABLE]

By virtue of (23) and Corollary 1,

[TABLE]

follows. Now, from (23) and Lemma 3, we have

[TABLE]

and this means that (30) tends to (12), for $h\rightarrow 0$. Consequently, ${\lambda_{0}}$ exists and is unique for all sufficiently small stepsizes $h>0$.

On the other hand, $g(q_{1})-g(q_{0})=g(q_{1})=0$, provided that (see (29))

[TABLE]

This means,

[TABLE]

and can be formally recast into the following linear system:

[TABLE]

Due to (23) and Corollary 1, the coefficient matrix reads:

[TABLE]

and the right-hand side is

[TABLE]

Consequently, (32) is consistent with (13), since  $\nabla g(q_{0})^{\top}M^{-1}p_{0}=0$, thus giving

[TABLE]

From Theorem 3, the conservation of energy follows.

The third statement of the theorem can be shown using the nonlinear variation of constants formula, by noting that for $t\in[0,h]$,

[TABLE]

The last result follows from

[TABLE]

Next, let us consider the mesh

[TABLE]

and the sequence of problems

[TABLE]

subject to initial conditions

[TABLE]

where $\lambda_{n}$ is a suitable constant vector. Then, the following result follows.

Theorem 5

Consider the IVPs (36)–(37) and let us denote by $(q(t),p(t))$ the solution of the problem (4)–(5). Moreover, let us assume that $(q_{n},p_{n})$ satisfies the following conditions:

• •

$q_{n}-q(t_{n})=O(h)$,   $p_{n}-p(t_{n})=O(h)$;

• •

$\nabla g(q_{n})^{\top}M^{-1}\nabla g(q_{n})\in\mathbb{R}^{\nu\times\nu}$* is nonsingular;*

• •

$g(q_{n})=0$;

• •

$\nabla g(q_{n})^{\top}M^{-1}p_{n}=O(h^{2})$.

Then, for all sufficiently small $h>0$, $\exists!\lambda_{n}\in\mathbb{R}^{\nu}$ such that the approximations

[TABLE]

satisfy:

• •

$g(q_{n+1})=0$* and, therefore, $H(q_{n+1},p_{n+1})=H(q_{n},p_{n})$;*

• •

$\lambda_{n}=\lambda(q(t_{n}),p(t_{n}))+O(h)$;

• •

$q_{n+1}-q(t_{n+1})=O(h)$,  $p_{n+1}-p(t_{n+1})=O(h)$;

• •

$\nabla g(q_{n+1})^{\top}M^{-1}p_{n+1}=O(h^{2})$.

Proof To show the first statement, we argue as we did to prove the first results in Theorem 4. This yields $g(q_{n+1})=0$ provided that, cf. (32)–(34),

[TABLE]

where $A(h)$ and $b(h)$ are defined as in (33)–(34) but with $q_{0}$ and $p_{0}$ replaced by $q_{n}$ and $p_{n}$, respectively. Consequently, from (14), we obtain

[TABLE]

Energy conservation follows, as before, from Theorem 3. Moreover, the nonlinear variation of constants formula, yields

[TABLE]

due to $\lambda_{n}-\lambda(q(t),p(t))=O(h)$, for $t\in[t_{n},t_{n+1}]$, and the hypothesis $q_{n}=q(t_{n})+O(h)$, $p_{n}=p(t_{n})+O(h)$.

In order to prove $\nabla g(q_{n+1})^{\top}M^{-1}p_{n+1}=O(h^{2})$, we note that from the hypothesis

[TABLE]

the existence of ${\tilde{p}}_{n}\in\mathbb{R}^{m}$ such that

[TABLE]

follows. Using $(q_{n},{\tilde{p}}_{n})$ as local initial conditions for (36), and repeating above steps to satisfy the constraints at $t_{n+1}$, we obtain

[TABLE]

and corresponding approximations ${\tilde{q}}_{n+1}$, ${\tilde{p}}_{n+1}$ such that

[TABLE]

From $p_{n}-{\tilde{p}}_{n}=O(h^{2})$ and $\lambda_{n}-{\tilde{\lambda}}_{n}=O(h)$, the nonlinear variation of constants formula yields,

[TABLE]

Consequently,

[TABLE]

By means of Theorem 5, a straightforward induction argument enables to show the following relaxed version of Theorem 2.

Corollary 2

For all sufficiently small stepsizes $h>0$, the above procedure defines a sequence of approximations $(q_{n},p_{n})$ such that, for all $n=1,2,\ldots,$

[TABLE]

Moreover, $(q_{n+1},p_{n+1})$ is obtained from $(q_{n},p_{n})$ utilizing a constant vector $\lambda_{n}$ such that

[TABLE]

where $\lambda(q,p)$ is the function defined in (14) and consequently,  $H(q_{n+1},p_{n+1})=H(q_{n},p_{n})$ holds.

The fact that (19) holds, in place of the weaker result (39), is due to the symmetry of the proposed procedure. Note that a symmetric method is necessarily of even order [27, Theorem 3.2]. The method is symmetric since we have shown that, for all sufficiently small $h>0$, there exists a unique $\lambda_{n}$ such that from the solution of (36)–(37), with

[TABLE]

we arrive at a new point, where

[TABLE]

Since $\lambda_{n}$ is unique, when we start from (41) and solve backward in time (36), we arrive at (40), i.e. the procedure is symmetric. As a consequence of the symmetry of the method, the approximation order of $(q_{n},p_{n})$ is even and therefore, (19) holds in place of (39). This completes the proof of Theorem 2.

In the next theorem, we summarize in a more comprehensive form the statements derived previously in this section.

Theorem 6

Let us consider the problem (4)–(6), the mesh (35), and the sequence of problems (36)–(37). The constant vector $\lambda_{n}$ is chosen in such a way that for the new approximations defined by (38), $g(q_{n+1})=0$ follows. Then, for all sufficiently small stepsizes $h>0$, the above procedure defines a sequence of approximations $(q_{n},p_{n},\lambda_{n})$ satisfying 111For the definition of $\lambda(q,p)$ see (14).

[TABLE]

Moreover, in case that $\lambda$ is constant, $\lambda(q(t),p(t))\equiv\bar{\lambda}$, $\forall t\in[0,T]$, the following statements hold:

[TABLE]

With other words, the discrete solution is exact, when the vector of the Lagrange multipliers is constant. Otherwise, it is second-order accurate for $(q_{n},p_{n})$ and first order accurate for $\lambda_{n}$. In the latter case, the constraints and the Hamiltonian are conserved, while the hidden constraints remain $O(h^{2})$ close to zero.

3 Polynomial approximation

The first step towards the numerical solution of (4)–(6) is to truncate the series in the right-hand side of (24),

[TABLE]

Here, the coefficients $\gamma_{j},\psi_{j}$, and $\rho_{j}$ are as those defined in (23). By imposing the initial conditions, the local approximation over the first step becomes

[TABLE]

with the new approximations which are formally still given by (18). Again, the constant vector of the multipliers is uniquely determined by requiring that the constraints are satisfied at $t_{1}=h$. Consequently, we have the following expression, in place of (2):

[TABLE]

which, in turn, yields equations which are formally similar to (32)–(34)222As before, this basic step defines a symmetric procedure.. The process is then repeated by defining the mesh (35) and considering the local problems

[TABLE]

where the coefficients $\gamma_{j}(v_{n}),\psi_{j}(u_{n}),\rho_{j}(u_{n})$ are defined in (23), with $u$ and $v$ replaced by $u_{n}$ and $v_{n}$, respectively. Consequently, we formally obtain the piecewise polynomial approximation,

[TABLE]

with the new approximations given by (see (21))

[TABLE]

As before, the constant vector $\lambda_{n}\in\mathbb{R}^{\nu}$ is chosen to satisfy the constraints $g(q_{n+1})=0$ and it is implicitly defined by the equation,

[TABLE]

This equation reduces to (3) for $n=0$. Using arguments similar to those from the previous section (see also [16]), it is possible to show the following result. This result is a counterpart to Theorem 6 for the piecewise polynomial approximation (46) to the solution $(q(t),p(t))$ of problem (4)–(6).

Theorem 7

For all sufficiently small stepsizes $h>0$, the approximation procedure (45)–(3) is well defined and provides a sequence of approximations $(q_{n},p_{n},\lambda_{n})$ such that

[TABLE]

Moreover, in case that $\lambda$ is constant, $\lambda(q(t),p(t))\equiv\bar{\lambda}$, $\forall t\in[0,T]$, the following statements hold:

[TABLE]

4 Full discretization

In order to cast the above algorithm into a computational method, the integrals defining the coefficients $\gamma_{j}(v),\,\psi_{j}(u),\,\rho_{j}(u),\,j=0,\dots,s-1,$ in (42), need to be approximated.333Since the method is a one-step method, we shall, as usual, only consider the first step. To this aim, following the discussion in [12, 16, 10], we use the Gauss-Legendre quadrature of order $2k$ (the interpolatory quadrature formula based at the zeros of $P_{k}(c)$), with nodes and weights $(\hat{c}_{i},\hat{b}_{i})$, where $k\geq s$. Consequently,

[TABLE]

Formally, this is a $k$-stage Runge-Kutta method, whose computational cost depends on $s$ rather than on $k$, since the actual unknowns are the $3s$ coefficients (49) and the vector ${\lambda_{0}}$ (see (42)). We refer, to [13, 10] for details.

Let us now formulate the discrete problem to be solved in each integration step. We first define the matrices, cf. (28),

[TABLE]

and the vectors and matrices

[TABLE]

Recall that $m$ is the dimension of the continuous problem and $\nu$ is the number of constraints.

Then, the $3s$ equations from (49), defining the discrete problem to be solved, amount to the system of equations, of (block) dimension $s$,

[TABLE]

We augment (50) by the equation (3) for ${\lambda_{0}}$ which, by taking (49) into account, can be rewritten as

[TABLE]

In (50), $\nabla U$, when evaluated in a block vector of (block) dimension $k$, stands for the block vector made up of the $k$ vectors resulting from the corresponding application of the function. The same straightforward notation is used for $\nabla g$. The new approximation is then given by, see (18) and (49),

[TABLE]

Note that the equations in (50), together with (52), formally coincide with those provided by a HBVM$(k,s)$ method444Here, $s$ is the degree of the polynomial approximation and $k$ defines the order (actually equal to $2k$) of the quadrature in the approximations (49). applied to solve the problem defined by the Hamiltonian (3), where the vector of the multiplier ${\lambda_{0}}$ is considered as a parameter,

[TABLE]

cf. [13, 10] for details. Consequently, equation (51) defines the proper extension for handling the constrained Hamiltonian problem (1)–(2). For this reason, we continue to refer to the numerical method specified in (50)–(52) as to HBVM$(k,s)$. Now, it is a ready to use numerical procedure.

The discrete problem (50)–(51) can be solved via a straightforward fixed-point iteration, which converges under regularity assumptions, for all sufficiently small stepsizes $h>0$.555We refer to [13, 7, 10] for further details on procedures for solving the involved discrete problems. They are based on suitable Newton-splitting procedures, already implemented in computational codes [18, 19, 20]. Moreover, for separable Hamiltonians, as it is the case in (1), the last two equations in (50) can be substituted into the first one, resulting in a single vector equation for $\hat{\bm{\gamma}}$. By setting

[TABLE]

we obtain

[TABLE]

plus (51) for ${\lambda_{0}}$.666Clearly, $\hat{\bm{\psi}}$ and $\hat{\bm{\rho}}$ can be computed via the last two equations in (50), once $\hat{\bm{\gamma}}$ and ${\lambda_{0}}$ are known. We skip further details, since they are exactly the same as for the original HBVMs, when applied to solve separable (unconstrained) Hamiltonian problems [13, 10].

The following result follows from Theorem 7 along with standard arguments from the analysis of HBVMs [16, 10].777For brevity, we omit the proof here.

Theorem 8

For all sufficiently small stepsizes $h>0$, the HBVM$(k,s)$ method (51)–(54) is well defined and symmetric. It provides a sequence of approximations $(q_{n},p_{n},\lambda_{n})$, $n=0,1,\dots,N$, such that

[TABLE]

and

[TABLE]

Moreover, in case $\lambda$ is constant, $\lambda(q(t),p(t))\equiv\bar{\lambda}$, $\forall t\in[0,T]$, the following statements hold:

[TABLE]

Remark 3

We stress that by (LABEL:Hg), an exact or a (at least) practical conservation of both the constraints and the Hamiltonian can always be guaranteed. In fact, by choosing sufficiently large $k$, either the quadrature becomes exact, in the polynomial case, or the quadrature error is within the round-off error level, in the non polynomial case. This feature of the method will be always exploited in the numerical tests discussed in Section 5.

Finally, we shall mention that for $k=s$, the HBVM$(s,s)$ method reduces to the $s$-stage Gauss collocation method, [12, 16, 10], which is symplectic. Moreover, in the limit $k\rightarrow\infty$, we retrieve the formulae studied in Section 3. This means that our approach can be also considered in the framework of Runge-Kutta methods with continuous stages [12, 26].

5 Numerical tests

In this section, to illustrate the theoretical properties of HBVM$(k,s)$, we apply them to numerically simulate some Hamiltonian problems of the form (1)–(2) with holonomic constraints. In the focus of our attention are properties described in Theorem 8.

5.1 Pendulum

We begin with the planar pendulum in Cartesian coordinates, where a massless rod of length $L$ connects a point of mass $m$ to a fixed point (the origin). We assume a unit mass and length, $m=1$ and $L=1$, and normalize the gravity acceleration. Then, the Hamiltonian is given by

[TABLE]

and is subject to the constraint

[TABLE]

We also prescribe the initial conditions of the form

[TABLE]

Consequently, the constrained Hamiltonian problem reads:

[TABLE]

In order to obtain a reference solution, we rewrite the problem in polar coordinates in such a way that $\theta=0$ locates the pendulum at its stable rest position, so that

[TABLE]

Thus, we arrive at the unconstrained Hamiltonian problem

[TABLE]

Once this problem is solved, the solution of (65) is recovered via the transformations (66). Moreover, the Lagrange multiplier in (63) turns out to be given by

[TABLE]

To compute the reference solution for (65), we solve (67) by means of a HBVM$(12,6)$ method 888For unconstrained Hamiltonian problems. of order 12 which is practically energy-conserving.

According to (62) and (63), the Hamiltonian and the constraint are quadratic, so we expect HBVM$(s,s)$ to conserve the energy and the constraint. In Table 1, we list the following quantities, obtained from the HBVM$(s,s)$ methods for $s=1,2,3$, the stepsizes $h=10^{-1}2^{-n}$ and the interval of integration $[0,10]$:

• •

the solution error ($e_{s}$),

• •

the multiplier error ($e_{\lambda}$),

• •

the Hamiltonian error $(e_{H})$,

• •

the constraint error $(e_{g})$;

• •

the hidden constraint error, defined by

[TABLE]

As predicted in Theorem 8, we can see that

• •

all methods are second-order accurate in the space variables, with HBVM(1,1) less accurate than the others;

• •

all methods are first-order accurate in the Lagrange multiplier;

• •

all methods exactly conserve the Hamiltonian and the constraint;

• •

all methods are second-order accurate in the hidden constraint.

5.2 Conical pendulum

Next, we consider the so-called conical pendulum, a particular case of the spherical pendulum, namely a pendulum of mass $m$, which is connected to a fixed point (i.e., the origin) by a massless rod of length $L$. For the conical pendulum, the initial condition is chosen in such a way that the motion is periodic with period $T$ and occurs in the horizontal plane $q_{3}=z_{0}$999Clearly, $0>z_{0}>-L$.. Again, assuming $m=1$ and $L=1$, and normalizing the acceleration of gravity, the Hamiltonian is

[TABLE]

with the constraint

[TABLE]

Here, we prescribe the consistent initial conditions

[TABLE]

generating a periodic motion with

[TABLE]

Moreover, in such a case, the multiplier ${\lambda_{0}}$, which has the physical meaning of the tension on the rod, has to be constant and is given by

[TABLE]

Note that the Hamiltonian and the constraint are quadratic and according to Theorem 8, any HBVM$(s,s)$ method conserves both of them and has order $2s$. In Table 2, we list the errors in

• •

the solution ($e_{s}$),

• •

the multiplier ($e_{\lambda}$),

• •

the Hamiltonian $(e_{H})$,

• •

the constraints $(e_{g})$;

• •

the hidden constraints, defined by

[TABLE]

The problem is solved over $10$ periods, with stepsizes $h=T/n$. As expected, the estimated rate of convergence for HBVM$(s,s)$, $s=1,2,3,4$, is $2s$. Also, the Hamiltonian and the constraint are conserved up to round-off errors. Remarkably, also the error in the multiplier $(e_{\lambda})$ and in the hidden constraints $(e_{hc})$ appear to be within the round-off error level, whatever stepsize is used.

In Figure 2, we plot the solution error from the computation over 100 periods using HBVM(2,2) with the stepsize $h=T/100\approx 0.053$, in Figure 2, the errors of the multiplier, Hamiltonian, constraint, and hidden constraints. One can see the linear growth of the solution error. The errors of the multiplier, constraint, Hamiltonian, and hidden constraints are negligible.

5.3 Modified pendulum

We now consider a modified version of the previous problem. With this simulation, we aim at exploiting the conservation of energy and constraints using a suitable high-order quadrature rule (49). More precisely, we consider the following non-quadratic polynomial Hamiltonian:

[TABLE]

with the non-quadratic polynomial constraint

[TABLE]

Here, we use the same initial condition as in (71) but the vector of the multipliers is no more constant, so that the order of the method reduces to $2$ (and $1$ for the vector of the multipliers). Moreover, since the constraints and the Hamiltonian are polynomials of degree not larger than $6$, any HBVM$(3s,s)$ method will conserve both quantities. This is confirmed by Table 3, where the results for HBVM(3,1), HBVM(6,2), and HBVM(9,3) are listed. The interval of integration was again $[0,10]$. Table 3 contains the errors in

• •

the solution ($e_{s}$),

• •

the multiplier ($e_{\lambda}$),

• •

the Hamiltonian $(e_{H})$,

• •

the constraints $(e_{g})$;

• •

the hidden constraints ($e_{hc}$), formally defined via (74).

Again, as predicted in Theorem 8, we see that

• •

all methods are second-order accurate in the state variables, with HBVM(3,1) less accurate than the others;

• •

all methods are first-order accurate in the Lagrange multiplier;

• •

all methods exactly conserve the Hamiltonian and the constraints;

• •

all methods are second-order accurate in the hidden constraints.

5.4 Tethered satellites system

Finally, we discuss a closed-loop rotating triangular tethered satellites system,101010This example can be found in [22, 39]. 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, which orbit around a massive body.111111The Earth, in the original example. As before, 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})^{\top}\in\mathbb{R}^{3}$, $i=1,2,3$, are the positions of the three satellites, the constraints are given by

[TABLE]

and the Hamiltonian is specified by

[TABLE]

The consistent initial conditions have the form

[TABLE]

and

[TABLE]

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=\frac{1}{2}$ and $y=-\frac{1}{2}$, respectively, and the third one moves around the tether joining the first two, in the plane $y=0$. In such a case, the Hamiltonian is non-polynomial. Nevertheless, using the HBVM(6,2) method with the stepsize $h=0.1$ over $10^{4}$ steps, we obtain a qualitatively correct solution which conserves the Hamiltonian and the constraints within the round-off error level, see Figure 3. Here, we also plot the hidden constraints errors $\|\nabla g(q_{n})^{\top}M^{-1}p_{n}\|$. At last, in Table 4, we list the following errors arising when solving the problem with HBVM$(6,s)$ methods for $s=1,2,3$ and the stepsizes $h=10^{-1}2^{-n}$, over the interval $[0,10]$:

• •

the solution error ($e_{s}$),

• •

the multipliers error ($e_{\lambda}$),

• •

the Hamiltonian error $(e_{H})$,

• •

the constraints error $(e_{g})$;

• •

the hidden constraints errors ($e_{hc}$), formally defined by (74).

Again, as shown in Theorem 8,

• •

all methods are second-order accurate in the state variables, with HBVM(6,1) much less accurate than the other two (which are almost equivalent);

• •

all methods are first-order accurate in the Lagrange multipliers;

• •

all methods exactly conserve the Hamiltonian and the constraints;

• •

all methods are second-order accurate in the hidden constraints.

To draw a general conclusion: it seems that using of HBVM$(k,s)$, with $s>1$, in context of the numerical solution of the Hamiltonian problems with holonomic constraints can be recommended, although the method is only second-order accurate.

6 Conclusions

In this paper, we have considered the numerical solution of Hamiltonian problems with holonomic constraints, by resorting to a line-integral formulation of the conservation of the constraints. This approach enables to derive an expression for the Lagrange multipliers in which second derivatives are not used. From the discretization of the resulting formulae, we have obtained a suitable variant of the Hamiltonian Boundary Value Methods (HBVMs), formerly designed as an energy-conserving Runge-Kutta methods for unconstrained Hamiltonian problems. Numerical experiments supporting the theoretical findings are enclosed.

This paper has been initiated in spring 2017, during the visit of the first author at the Institute for Analysis and Scientific Computing, Vienna University of Technology, Vienna, Austria.