Line Integral Solution of Hamiltonian Systems with Holonomic Constraints
Luigi Brugnano, Gianmarco Gurioli, Felice Iavernaro, Ewa B. Weinmueller
aa r X i v : . [ m a t h . NA ] D ec Line Integral Solution of Hamiltonian Systems withHolonomic Constraints
Luigi Brugnano a , Gianmarco Gurioli a , Felice Iavernaro b , Ewa B. Weinm¨uller ca Dipartimento di Matematica e Informatica “U. Dini”, Universit`a di FirenzeViale Morgagni 67/a, I-50134 Firenze, Italy. b Dipartimento di Matematica, Universit`a di BariVia Orabona 4, I-70125 Bari, Italy. c Institute for Analysis and Scientific Computing, Vienna University of TechnologyA-1040 Wien, Austria. – Dedicated to John Butcher, on the occasion of his 84-th birthday –
Abstract
In this paper, we propose a second-order energy-conserving approximation procedure forHamiltonian systems with holonomic constraints. The derivation of the procedure relies on theuse of the so-called line integral framework . We provide numerical experiments to illustratetheoretical findings.
Keywords: constrained Hamiltonian systems; holonomic constraints; energy-conserving meth-ods; line integral methods; Hamiltonian Boundary Value Methods; HBVMs.
MSC:
We consider the numerical approximation of a constrained Hamiltonian dynamics, described by theseparable Hamiltonian H ( q, p ) = 12 p ⊤ M − p + U ( q ) , q, p ∈ R m , (1)where M is a symmetric and positive-definite matrix. The problem is completed by ν holonomic constraints, g ( q ) = 0 ∈ R ν , (2)where we assume that ν < m holds. Moreover, we also assume that all points are regular forthe constraints, i.e., ∇ g ( q ) ∈ R m × ν has full column rank or, equivalently, ∇ g ( q ) ⊤ M − ∇ g ( q ) isnonsingular. 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 definingthe augmented Hamiltonian ˆ H ( q, p, λ ) = H ( q, p ) + λ ⊤ g ( q ) , (3)1here λ is the vector of Lagrange multipliers. The resulting constrained Hamiltonian system reads:˙ q = M − p, ˙ p = −∇ U ( q ) − ∇ g ( q ) λ, g ( q ) = 0 , t ∈ [0 , T ] , (4)and is subject to consistent initial conditions, q (0) = q , p (0) = p , (5)such that g ( q ) = 0 , ∇ g ( q ) ⊤ M − p = 0 . (6)Note that the condition g ( q ) = 0 ensures that q belongs to the manifold M = { q ∈ R m : g ( q ) = 0 } , (7)as required by the constraints, whereas the condition ∇ g ( q ) ⊤ M − p means that the motion initiallystays on the tangent space to M at q . On a continuous level, this condition is satisfied by all pointson the solution trajectory, since, in order for the constraints to be conserved,˙ g ( q ) = ∇ g ( q ) ⊤ ˙ q = ∇ g ( q ) ⊤ M − p = 0 , (8)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 ∇ g ( q ) ⊤ M − p to be suitably small along the numerical solution. Conse-quently, when solving the problem on the interval [0 , h ], we require that the approximations, q ≈ q ( h ) , p ≈ p ( h ) , (9)satisfy the conservation of both, the Hamiltonian and the constraints, H ( q , p ) = H ( q , p ) , g ( q ) = g ( q ) = 0 , (10)and that the hidden constraints are relaxed to ∇ g ( q ) ⊤ M − p = O ( h ) . (11)We recall that a formal expression for the vector λ is obtained by an additional differentiating of(8), i.e., ¨ g ( q ) = ∇ g ( q )( M − p, M − p ) − ∇ g ( q ) ⊤ M − [ ∇ U ( q ) + ∇ g ( q ) λ ] . (12)Imposing the vanishing of this derivative yields (cid:2) ∇ g ( q ) ⊤ M − ∇ g ( q ) (cid:3) λ = ∇ g ( q )( M − p, M − p ) − ∇ g ( q ) ⊤ M − ∇ U ( q ) . (13)Consequently, the following result follows. Theorem 1
The vector λ exists and is uniquely determined, provided that the matrix ∇ g ( q ) ⊤ M − ∇ g ( q ) is nonsingular. In fact, in such a case, from (13), we obtain λ = (cid:2) ∇ g ( q ) ⊤ M − ∇ g ( q ) (cid:3) − (cid:2) ∇ g ( q )( M − p, M − p ) − ∇ g ( q ) ⊤ M − ∇ U ( q ) (cid:3) = : λ ( q, p ) . (14) Note that, for later use, we have introduced the notation λ ( q ( t ) , p ( t )) in place of λ ( t ) , to explicitlyunderline the dependence of the Lagrange multiplier on the state variables q and p at time t . emark 1 We observe that an additional differentiation of (13) provides a differential equationfor the Lagrange multipliers, which can be solved together with the original problem. However, thisprocedure 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 timein 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 obtainedvia symplectic PRK methods [24], composition methods [34, 35], symmetric LMFs [23]. Furthermethods are based on discrete derivatives [33], local parametrizations of the manifold containingthe 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 hasalready been used when deriving the energy-conserving Runge-Kutta methods, for unconstrainedHamiltonian systems, cf. Hamiltonian Boundary Value Methods (HBVMs) [11, 12, 13, 16, 17] andthe 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). Roughlyspeaking, the conservation of the invariant will be guaranteed by requiring that a suitable lineintegral vanishes. This line integral represents a discrete-time version of (8). In fact, if we fix astepsize h >
0, then the conservation of the constraints (2) at h , starting from the point q definedin (5), can be recast into g ( q ( h )) − g ( q (0)) | {z } =0 = Z h ∇ g ( q ( t )) ⊤ ˙ q ( t )d t = 0 . 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 describinga discrete-time dynamics. In such a case, the conservation properties have to be satisfied only ona set of discrete times which are multiples of the stepsize h . Consequently, we consider a localapproximation to q ( t ), say u ( t ), such that u (0) = q , u ( h ) =: q ≈ q ( h ) , and g ( q ) − g ( q ) ≡ g ( u ( h )) − g ( u (0)) = Z h ∇ g ( u ( t )) ⊤ ˙ u ( t )d t = 0 , without requiring the integrand to be identically zero. This, in turn, enables a proper choice ofthe vector of the multipliers λ . As a result, we eventually obtain suitably modified HBVMs whichenable to conserve both, the Hamiltonian and the constraints. We stress that the available efficientimplementation of the original methods (see, e.g., [13, 7, 10]), which proved to be reliable androbust in the numerical solution of the unconstrained Hamiltonian problems, can now be adaptedfor dealing with the holonomic constraints.The paper is organized as follows. In Section 2, we provide the framework for devising the methodvia a suitable choice of the vector λ of the Lagrange multipliers, which we approximate by apiecewise-constant function. In Section 3 further simplification towards numerical procedure isdiscussed. Then, in Section 4, we present a fully discrete method, resulting in a suitable modificationof the original HBVMs. In Section 5, numerical experiments are shown to illustrate how the methodworks for a number of constrained Hamiltonian problems. Section 6 contains the conclusions andpossible future investigations. 3 Piecewise-constant approximation of λ 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 λ ∈ R ν such that (10) is satisfied. This is equivalent torequire ˆ H ( q , p , λ ) = ˆ H ( q , p , λ ) , g ( q ) = g ( q ) = 0 , (15)where the constant parameter λ is chosen in such a way that the constraints g ( q ) = 0 hold. We willshow 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 theproblem (4)–(5), by the local problem˙ u = M − v, ˙ v = −∇ U ( u ) − ∇ g ( u ) λ , t ∈ [0 , h ] , (16)subject to the initial conditions, cf. (5), u (0) = q , v (0) = p , (17)satisfying (6). By setting q := u ( h ) , p := v ( h ) , (18)the constant parameter λ is chosen to guarantee the conservation of the Hamiltonian and theconstraints, i.e., (10). Starting from (18), the procedure is then repeated on [ h, h ] and the followingintervals. The convergence result is now formulated in the following theorem. Theorem 2
For all sufficiently small stepsizes h > , the above procedure defines a sequence ofapproximations ( q n , p n ) such that, for all n = 1 , , . . . : q n = q ( nh ) + O ( h ) , p n = p ( nh ) + O ( h ) , g ( q n ) = 0 , ∇ g ( q n ) ⊤ M − p n = O ( h ) . (19) Moreover, ( q n +1 , p n +1 ) is obtained from ( q n , p n ) using a constant vector λ n such that λ n = λ ( q ( nh ) , p ( nh )) + O ( h ) , (20) where λ ( 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 , { P j } , P j ∈ Π j , Z P i ( c ) P j ( c )d c = δ ij , ∀ i, j = 0 , , . . . , (21)along with the expansions, M − v ( ch ) = X j ≥ P j ( c ) γ j ( v ) , ∇ U ( u ( ch )) = X j ≥ P j ( c ) ψ j ( u ) , ∇ g ( u ( ch )) = X j ≥ P j ( c ) ρ j ( u ) , c ∈ [0 , , (22)4ith γ j ( v ) = M − Z P j ( c ) v ( ch )d c, ψ j ( u ) = Z P j ( c ) ∇ U ( u ( ch ))d c,ρ j ( u ) = Z P j ( c ) ∇ g ( u ( ch ))d c, j ≥ . (23)Consequently, following the approach defined in [16], the differential equations in (16) can berewritten as˙ u ( ch ) = X j ≥ P j ( c ) γ j ( v ) , ˙ v ( ch ) = − X j ≥ P j ( c )[ ψ j ( u ) + ρ j ( u ) λ ] , c ∈ [0 , . (24)Moreover, using the initial conditions (5), we formally obtain u ( ch ) = q + h X j ≥ Z c P j ( x )d x γ j ( v ) , v ( ch ) = p − h X j ≥ Z c P j ( x )d x [ ψ j ( u ) + ρ j ( u ) λ ] , c ∈ [0 , . (25)The following result is now cited from [16, Lemma 1]. Lemma 1
Let G : [0 , h ] → V , with V a vector space, admit a Taylor expansion at . Then Z P j ( c ) G ( ch )d c = O ( h j ) , j ≥ . 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 λ ∈ R ν , the solution of (16)–(18) satisfies ˆ H ( q , p , λ ) = ˆ H ( q , p , λ ) . Proof For any given λ ∈ R ν , it follows from (3), (16), and (18),ˆ H ( q , p , λ ) − ˆ H ( q , p , λ ) = ˆ H ( u ( h ) , v ( h ) , λ ) − ˆ H ( u (0) , v (0) , λ )= Z h dd t ˆ H ( u ( t ) , v ( t ) , λ )d t = Z h n ˆ H q ( u ( t ) , v ( t ) , λ ) ⊤ ˙ u ( t ) + ˆ H p ( u ( t ) , v ( t ) , λ ) ⊤ ˙ v ( t ) o d t = h Z (cid:8) [ ∇ U ( u ( ch )) + ∇ g ( u ( ch )) λ ] ⊤ ˙ u ( ch ) + [ M − v ( ch )] ⊤ ˙ v ( ch ) (cid:9) d c = h Z [ ∇ U ( u ( ch )) + ∇ g ( u ( ch )) λ ] ⊤ X j ≥ P j ( c ) γ j ( v ) − [ M − v ( ch )] ⊤ X j ≥ P j ( c )[ ψ j ( u ) + ρ j ( u ) λ ] d c h X j ≥ ((cid:18)Z P j ( c )[ ∇ U ( u ( ch )) + ∇ g ( u ( ch )) λ ]d c (cid:19) ⊤ γ j ( v ) − (cid:18) M − Z P j ( c ) v ( ch )d c (cid:19) ⊤ [ ψ j ( u ) + ρ j ( u ) λ ] ) = h X j ≥ (cid:8) [ ψ j ( u ) + ρ j ( u ) λ ] ⊤ γ j ( v ) − γ j ( v ) ⊤ [ ψ j ( u ) + ρ j ( u ) λ ] (cid:9) = 0 . (cid:3) As observed above, the conservation of the Hamiltonian (1) is guaranteed, once the constraints aresatisfied, i.e., g ( q ) = 0. We now apply a line integral technique to determine the vector λ andformulate 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 , p ) is given such that, • ∇ g ( q ) ⊤ M − ∇ g ( q ) ∈ R ν × ν is nonsingular; • g ( q ) = 0 ; • ∇ g ( q ) ⊤ M − p = 0 .Then, for all sufficiently small h > , ∃ ! λ ∈ R ν such that the approximations in (18) satisfy • g ( q ) = 0 and, therefore, H ( q , p ) = H ( q , p ) ; • λ = λ ( q , p ) + O ( h ) ; • q − q ( h ) = O ( h ) , p − p ( h ) = O ( h ) ; • ∇ g ( q ) ⊤ M − p = O ( h ) . 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 Z P j ( c ) Z c P i ( x )d x d c = ( X s ) j +1 ,i +1 , i, j = 0 , . . . , s − , (26) where ( X s ) j +1 ,i +1 is the ( j + 1 , i + 1) entry of the matrix X s := ξ − ξ ξ . . .. . . . . . − ξ s − ξ s − , ξ j = 12 p | j − | , j = 0 , . . . , s − . (27)6roof Since the integrand on the left-hand side in (26) is a polynomial of degree at most 2 s − s . Let c , . . . , c s bethe zeros of P s and b , . . . , b s be the corresponding weights. Then, introducing the matrices P s = ( P j − ( c i )) , I s = (cid:18)Z c i P j − ( x )d x (cid:19) , Ω = diag( b , . . . , b s ) ∈ R s × s , (28)and setting e i ∈ R s , the i -th unit vector, we have Z P j ( c ) Z c P i ( x )d x d c = s X ℓ =1 b ℓ P j ( c ℓ ) Z c ℓ P i ( x )d x ≡ e ⊤ j +1 P ⊤ s Ω I s e i +1 . The result follows by observing that, due to the properties of Legendre polynomials [10, Sec-tion 1.4.3], I s = P s X s , P ⊤ s Ω P s = I s follows, where X s is the matrix defined in (27), and I s ∈ R s × s is the identity matrix. This yields e ⊤ j +1 P ⊤ s Ω I s e i +1 = e ⊤ j +1 P ⊤ s Ω P s X s e i +1 = e ⊤ j +1 X s e i +1 . (cid:3) We also need the following expansions.
Lemma 3
From (16) and (23), we conclude ρ ( u ) = ∇ g ( q ) + h ∇ g ( q ) M − p + O ( h ) ,ρ ( u ) ⊤ M − p = ∇ g ( q ) ⊤ M − p + h ∇ g ( q )( M − p , M − p ) + O ( h ) . Proof The first equality follows from Lemma 1 and Corollary 1, ρ ( u ) = Z ∇ g ( u ( ch ))d c = Z (cid:2) ∇ g ( u (0)) + ch ∇ g ( u (0)) ˙ u (0) + O (( ch ) ) (cid:3) d c = ∇ g ( u (0)) + h ∇ g ( u (0)) ˙ u (0) + O ( h )= ∇ g ( q ) + h ∇ g ( q ) X j ≥ P j (0) γ j ( v ) + O ( h )= ∇ g ( q ) + h ∇ g ( q ) γ ( v ) + O ( h )= ∇ g ( q ) + h ∇ g ( q ) M − [ p + O ( h )] + O ( h )= ∇ g ( q ) + h ∇ g ( q ) M − p + O ( h ) . The second statement follows by transposition and multiplication from the right by M − p . (cid:3) We now show the results formulated in Theorem 4.7roof (of Theorem 4). Let us assume that g ( q ) = 0 holds. Then, it follows from (16)–(18), g ( q ) = g ( q ) − g ( q ) = g ( u ( h )) − g ( u (0)) = Z h dd t g ( u ( t ))d t = Z h ∇ g ( u ( t )) ⊤ ˙ u ( t )d t = h Z ∇ g ( u ( ch )) ⊤ ˙ u ( ch )d c = h Z ∇ g ( u ( ch )) ⊤ X j ≥ P j ( c ) γ j ( v )d c = h X j ≥ ρ j ( u ) ⊤ γ j ( v ) = h X j ≥ ρ j ( u ) ⊤ M − Z P j ( c ) v ( ch )d c = h X j ≥ ρ j ( u ) ⊤ M − Z P j ( c ) p − h X i ≥ Z c P i ( x )d x [ ψ i ( u ) + ρ i ( u ) λ ] d c = h X j ≥ ρ j ( u ) ⊤ M − p Z P j ( c )d c − h X i,j ≥ ρ j ( u ) ⊤ M − [ ψ i ( u ) + ρ i ( u ) λ ] Z P j ( c ) Z c P i ( x )d x d c. Due to (21), Z P j ( c )d c = δ j , and according to (26)–(27), we conclude g ( u ( h )) − g ( u (0)) = hρ ( u ) ⊤ M − { p − h [ ξ ( ψ ( u ) + ρ ( u ) λ ) − ξ ( ψ ( u ) + ρ ( u ) λ )] }− h X j ≥ ρ j ( u ) ⊤ M − { [ ξ j ( ψ j − ( u ) + ρ j − ( u ) λ ) − ξ j +1 ( ψ j +1 ( u ) + ρ j +1 ( u ) λ )] } =: ˆΓ( u, v, λ , h ) . (29)By virtue of (23) and Corollary 1, g ( u ( h )) − g ( u (0)) − hρ ( u ) ⊤ M − p h = − (cid:8)(cid:2) ρ ( u ) ⊤ M − ρ ( u ) + O ( h ) (cid:3) λ + ρ ( u ) ⊤ M − ψ ( u ) + O ( h ) (cid:9) (30)follows. Now, from (23) and Lemma 3, we have g ( u ( h )) − g ( u (0)) − hρ ( u ) ⊤ M − p h = 12 (cid:8) ¨ g ( q ) − ∇ g ( q )( M − p , M − p ) + O ( h ) (cid:9) ,ρ ( u ) ⊤ M − ρ ( u ) = ∇ g ( q ) ⊤ M − ∇ g ( q ) + O ( h ) ,ρ ( u ) ⊤ M − ψ ( u ) = ∇ g ( q ) ⊤ M − ∇ U ( q ) + O ( h ) , and this means that (30) tends to (12), for h →
0. Consequently, λ exists and is unique for allsufficiently small stepsizes h >
0. 8n the other hand, g ( q ) − g ( q ) = g ( q ) = 0, provided that (see (29))ˆΓ( u, v, λ , h ) = 0 . This means, ρ ( u ) ⊤ M − p (31)= h X j ≥ ρ j ( u ) ⊤ M − (cid:8) ξ j [ ψ j − δ j ( u ) + ρ j − δ j ( u ) λ ] − ξ j +1 [ ψ j +1 ( u ) + ρ j +1 ( u ) λ ] (cid:9) = h ξ ρ ( u ) ⊤ M − ρ ( u ) + X j ≥ ξ j (cid:2) ρ j ( u ) ⊤ M − ρ j − ( u ) − ρ j − ( u ) ⊤ M − ρ j ( u ) (cid:3) λ + ξ ρ ( u ) ⊤ M − ψ ( u ) + X j ≥ ξ j (cid:2) ρ j ( u ) ⊤ M − ψ j − ( u ) − ρ j − ( u ) ⊤ M − ψ j ( u ) (cid:3) , and can be formally recast into the following linear system: A ( h ) λ = b ( h ) . (32)Due to (23) and Corollary 1, the coefficient matrix reads: A ( h ) = hξ ρ ( u ) ⊤ M − ρ ( u ) + O ( h ) ≡ h ∇ g ( q ) ⊤ M − ∇ g ( q ) + O ( h ) , (33)and the right-hand side is b ( h ) = ρ ( u ) ⊤ M − p − ξ hρ ( u ) ⊤ M − ψ ( u ) + O ( h ) ≡ ∇ g ( q ) ⊤ M − p + h (cid:2) ∇ g ( q )( M − p , M − p ) − ∇ g ( q ) ⊤ M − ∇ U ( q ) (cid:3) + O ( h ) . (34)Consequently, (32) is consistent with (13), since ∇ g ( q ) ⊤ M − p = 0, thus giving λ == (cid:2) ∇ g ( q ) ⊤ M − ∇ g ( q ) + O ( h ) (cid:3) − (cid:2) ∇ g ( q )( M − p , M − p ) − ∇ g ( q ) ⊤ M − ∇ U ( q ) + O ( h ) (cid:3) ≡ λ ( q , p ) + O ( h ) . From Theorem 3, the conservation of energy follows.The third statement of the theorem can be shown using the nonlinear variation of constantsformula, by noting that for t ∈ [0 , h ], λ ( t ) − λ ≡ λ ( q ( t ) , p ( t )) − λ ( q (0) , p (0)) | {z } = O ( h ) + O ( h ) = O ( h ) . The last result follows from ∇ g ( q ) ⊤ M − p = ∇ g ( q ( h ) + O ( h )) ⊤ M − ( p ( h ) + O ( h ))= ∇ g ( q ( h )) ⊤ M − p ( h ) | {z } =0 + O ( h ) = O ( h ) . (cid:3) t n = nh, n = 0 , . . . , N, (35)and the sequence of problems˙ u = M − v, ˙ v = −∇ U ( u ) − ∇ g ( u ) λ n , t ∈ [ t n , t n +1 ] , (36)subject to initial conditions u ( t n ) = q n , v ( t n ) = p n , (37)where λ 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 prob-lem (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 ) ; • ∇ g ( q n ) ⊤ M − ∇ g ( q n ) ∈ R ν × ν is nonsingular; • g ( q n ) = 0 ; • ∇ g ( q n ) ⊤ M − p n = O ( h ) .Then, for all sufficiently small h > , ∃ ! λ n ∈ R ν such that the approximations q n +1 := u ( t n +1 ) , p n +1 := v ( t n +1 ) , (38) satisfy: • g ( q n +1 ) = 0 and, therefore, H ( q n +1 , p n +1 ) = H ( q n , p n ) ; • λ n = λ ( 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 ) ; • ∇ g ( q n +1 ) ⊤ M − p n +1 = O ( h ) . 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), A ( h ) λ n = b ( h )where A ( h ) and b ( h ) are defined as in (33)–(34) but with q and p replaced by q n and p n , respec-tively. Consequently, from (14), we obtain λ n = (cid:2) ∇ g ( q n ) ⊤ M − ∇ g ( q n ) + O ( h ) (cid:3) − ∇ g ( q n )( M − p n , M − p n ) − ∇ g ( q n ) ⊤ M − ∇ U ( q n ) − h = O ( h ) z }| { ∇ g ( q n ) ⊤ M − p n + O ( h ) = (cid:2) ∇ g ( q n ) ⊤ M − ∇ g ( q n ) + O ( h ) (cid:3) − (cid:2) ∇ g ( q n )( M − p n , M − p n ) − ∇ g ( q n ) ⊤ M − ∇ U ( q n ) + O ( h ) (cid:3) ≡ λ ( q n , p n ) + O ( h ) = λ ( q ( t n ) + O ( h ) , p ( t n ) + O ( h )) + O ( h ) = λ ( q ( t n ) , p ( t n )) + O ( h ) . (cid:18) q n +1 − q ( t n +1 ) p n +1 − p ( t n +1 ) (cid:19) = O ( h ) + O ( h ) = O ( h ) , due to λ n − λ ( q ( t ) , p ( t )) = O ( h ), for t ∈ [ 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 ∇ g ( q n +1 ) ⊤ M − p n +1 = O ( h ), we note that from the hypothesis ∇ g ( q n ) ⊤ M − p n = O ( h ) , the existence of ˜ p n ∈ R m such that p n − ˜ p n = O ( h ) , ∇ g ( q n ) ⊤ M − ˜ p n = 0follows. Using ( q n , ˜ p n ) as local initial conditions for (36), and repeating above steps to satisfy theconstraints at t n +1 , we obtain˜ λ n = λ ( q n , ˜ p n ) + O ( h ) ≡ λ ( q n , p n + O ( h )) + O ( h ) = λ ( q n , p n ) + O ( h ) ≡ λ n + O ( h ) , and corresponding approximations ˜ q n +1 , ˜ p n +1 such that g (˜ q n +1 ) = 0 , ∇ g (˜ q n +1 ) ⊤ M − ˜ p n +1 = O ( h ) . From p n − ˜ p n = O ( h ) and λ n − ˜ λ n = O ( h ), the nonlinear variation of constants formula yields, q n +1 − ˜ q n +1 = O ( h ) , p n +1 − ˜ p n +1 = O ( h ) . Consequently, ∇ g ( q n +1 ) ⊤ M − p n +1 = ∇ g (˜ q n +1 + O ( h )) ⊤ M − (˜ p n +1 + O ( h ))= ∇ g (˜ q n +1 ) ⊤ M − ˜ p n +1 | {z } = O ( h ) + O ( h ) = O ( h ) . (cid:3) By means of Theorem 5, a straightforward induction argument enables to show the followingrelaxed version of Theorem 2.
Corollary 2
For all sufficiently small stepsizes h > , the above procedure defines a sequence ofapproximations ( q n , p n ) such that, for all n = 1 , , . . . ,q n = q ( nh ) + O ( h ) , p n = p ( nh ) + O ( h ) , g ( q n ) = 0 , ∇ g ( q n ) ⊤ M − p n = O ( h ) . (39) Moreover, ( q n +1 , p n +1 ) is obtained from ( q n , p n ) utilizing a constant vector λ n such that λ n = λ ( q ( nh ) , p ( nh )) + O ( h ) , where λ ( q, p ) is the function defined in (14) and consequently, H ( q n +1 , p n +1 ) = H ( q n , p n ) holds. symmetry of theproposed 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 aunique λ n such that from the solution of (36)–(37), with( q n , p n ) , g ( q n ) = 0 , (40)we arrive at a new point, where ( q n +1 , p n +1 ) , g ( q n +1 ) = 0 , (41)Since λ 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 approximationorder of ( q n , p n ) is even and therefore, (19) holds in place of (39). This completes the proof ofTheorem 2.In the next theorem, we summarize in a more comprehensive form the statements derivedpreviously 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 λ n is chosen in such a way that for the new approximations definedby (38), g ( q n +1 ) = 0 follows. Then, for all sufficiently small stepsizes h > , the above proceduredefines a sequence of approximations ( q n , p n , λ n ) satisfying q n = q ( nh ) + O ( h ) ,p n = p ( nh ) + O ( h ) ,g ( q n ) = 0 , ∇ g ( q n ) ⊤ M − p n = O ( h ) ,λ n = λ ( q ( nh ) , p ( nh )) + O ( h ) ,H ( q n , p n ) = H ( q , p ) , n = 0 , , . . . , N. Moreover, in case that λ is constant, λ ( q ( t ) , p ( t )) ≡ ¯ λ , ∀ t ∈ [0 , T ] , the following statements hold: q n = q ( nh ) ,p n = p ( nh ) ,g ( q n ) = 0 , ∇ g ( q n ) ⊤ M − p n = 0 ,λ n = ¯ λ,H ( q n , p n ) = H ( q , p ) , n = 0 , , . . . , N. With other words, the discrete solution is exact , when the vector of the Lagrange multipliers isconstant. Otherwise, it is second-order accurate for ( q n , p n ) and first order accurate for λ n . In thelatter case, the constraints and the Hamiltonian are conserved, while the hidden constraints remain O ( h ) close to zero. For the definition of λ ( q, p ) see (14). Polynomial approximation
The first step towards the numerical solution of (4)–(6) is to truncate the series in the right-handside of (24),˙ u ( ch ) = s − X j =0 P j ( c ) γ j ( v ) , ˙ v ( ch ) = − s − X j =0 P j ( c )[ ψ j ( u ) + ρ j ( u ) λ ] , c ∈ [0 , . (42)Here, the coefficients γ j , ψ j , and ρ j are as those defined in (23). By imposing the initial conditions,the local approximation over the first step becomes u ( ch ) = q + h s − X j =0 Z c P j ( x )d x γ j ( v ) , v ( ch ) = p − h s − X j =0 Z c P j ( x )d x [ ψ j ( u ) + ρ j ( u ) λ ] , c ∈ [0 , , (43)with the new approximations which are formally still given by (18). Again, the constant vectorof the multipliers is uniquely determined by requiring that the constraints are satisfied at t = h .Consequently, we have the following expression, in place of (31): ρ ( u ) ⊤ M − p (44)= h s − X j =0 ρ j ( u ) ⊤ M − (cid:8) ξ j [ ψ j − δ j ( u ) + ρ j − δ j ( u ) λ ] − ξ j +1 [ ψ j +1 ( u ) + ρ j +1 ( u ) λ ] (cid:9) + hξ s − ρ s − ( u ) ⊤ M − [ ψ s − ( u ) + ρ s − ( u ) λ ]= h ξ ρ ( u ) ⊤ M − ρ ( u ) + s − X j =1 ξ j (cid:2) ρ j ( u ) ⊤ M − ρ j − ( u ) − ρ j − ( u ) ⊤ M − ρ j ( u ) (cid:3) λ + ξ ρ ( u ) ⊤ M − ψ ( u ) + s − X j =1 ξ j (cid:2) ρ j ( u ) ⊤ M − ψ j − ( u ) − ρ j − ( u ) ⊤ M − ψ j ( u ) (cid:3) , which, in turn, yields equations which are formally similar to (32)–(34) . The process is thenrepeated by defining the mesh (35) and considering the local problems˙ u n ( ch ) = s − X j =0 P j ( c ) γ j ( v n ) , ˙ v n ( ch ) = − s − X j =0 P j ( c )[ ψ j ( u n ) + ρ j ( u n ) λ n ] , c ∈ [0 , ,u n (0) = q n , v n (0) = p n , n = 0 , . . . , N − , (45)where the coefficients γ j ( v n ) , ψ j ( u n ) , ρ 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, u n ( ch ) = q n + h s − X j =0 Z c P j ( x )d xγ j ( v n ) , (46) v n ( ch ) = p n − h s − X j =0 Z c P j ( x )d x [ ψ j ( u n ) + ρ j ( u n ) λ n ] , c ∈ [0 , , As before, this basic step defines a symmetric procedure. q n +1 := u n ( h ) ≡ q n + hγ ( v n ) , p n +1 := v n ( h ) ≡ p n − h [ ψ ( u n ) + ρ ( u n ) λ n ] . (47)As before, the constant vector λ n ∈ R ν is chosen to satisfy the constraints g ( q n +1 ) = 0 and it isimplicitly defined by the equation, ρ ( u n ) ⊤ M − p n (48)= h ξ ρ ( u n ) ⊤ M − ρ ( u n ) + s − X j =1 ξ j (cid:2) ρ j ( u n ) ⊤ M − ρ j − ( u n ) − ρ j − ( u n ) ⊤ M − ρ j ( u n ) (cid:3) λ n + ξ ρ ( u n ) ⊤ M − ψ ( u n ) + s − X j =1 ξ j (cid:2) ρ j ( u n ) ⊤ M − ψ j − ( u n ) − ρ j − ( u n ) ⊤ M − ψ j ( u n ) (cid:3) . This equation reduces to (44) 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 6for the piecewise polynomial approximation (46) to the solution ( q ( t ) , p ( t )) of problem (4)–(6). Theorem 7
For all sufficiently small stepsizes h > , the approximation procedure (45)–(48) iswell defined and provides a sequence of approximations ( q n , p n , λ n ) such that q n = q ( nh ) + O ( h ) ,p n = p ( nh ) + O ( h ) ,g ( q n ) = 0 , ∇ g ( q n ) ⊤ M − p n = O ( h ) ,λ n = λ ( q ( nh ) , p ( nh )) + O ( h ) ,H ( q n , p n ) = H ( q , p ) , n = 0 , , . . . , N. Moreover, in case that λ is constant, λ ( q ( t ) , p ( t )) ≡ ¯ λ , ∀ t ∈ [0 , T ] , the following statements hold: q n = q ( nh ) + O ( h s ) ,p n = p ( nh ) + O ( h s ) ,g ( q n ) = 0 , ∇ g ( q n ) ⊤ M − p n = O ( h s ) ,λ n = ¯ λ + O ( h s ) ,H ( q n , p n ) = H ( q , p ) , n = 0 , , . . . , N. In order to cast the above algorithm into a computational method, the integrals defining the coeffi-cients γ j ( v ) , ψ j ( u ) , ρ j ( u ) , j = 0 , . . . , s − , in (42), need to be approximated. To this aim, followingthe discussion in [12, 16, 10], we use the Gauss-Legendre quadrature of order 2 k (the interpolatory Since the method is a one-step method, we shall, as usual, only consider the first step. P k ( c )), with nodes and weights (ˆ c i , ˆ b i ), where k ≥ s .Consequently, γ j ( v ) ≈ ˆ γ j := M − k X ℓ =1 ˆ b ℓ P j (ˆ c ℓ ) v (ˆ c ℓ h ) , ψ j ( u ) ≈ ˆ ψ j := k X ℓ =1 ˆ b ℓ P j (ˆ c ℓ ) ∇ U ( u (ˆ c ℓ h )) ,ρ j ( u ) ≈ ˆ ρ j := k X ℓ =1 ˆ b ℓ P j (ˆ c ℓ ) ∇ g ( u (ˆ c ℓ h )) , j = 0 , . . . , s − . (49)Formally, this is a k -stage Runge-Kutta method, whose computational cost depends on s ratherthan on k , since the actual unknowns are the 3 s coefficients (49) and the vector λ (see (42)). Werefer, to [13, 10] for details.Let us now formulate the discrete problem to be solved in each integration step. We first definethe matrices, cf. (28),ˆ P s = ( P j − (ˆ c i )) , ˆ I s = Z ˆ c i P j − ( x )d x ! ∈ R k × s , ˆΩ = diag(ˆ b , . . . , ˆ b k ) ∈ R k × k , and the vectors and matrices e = ∈ R k , ˆ γ = ˆ γ ...ˆ γ s − , ˆ ψ = ˆ ψ ...ˆ ψ s − ∈ R sm , ˆ ρ = ˆ ρ ...ˆ ρ s − ∈ R sm × ν . Recall that m is the dimension of the continuous problem and ν is the number of constraints.Then, the 3 s equations from (49), defining the discrete problem to be solved, amount to thesystem of equations, of (block) dimension s ,ˆ γ = ˆ P ⊤ s ˆΩ ⊗ M − h e ⊗ p − h ˆ I s ⊗ I m (cid:16) ˆ ψ + ˆ ρ λ (cid:17)i , ˆ ψ = ˆ P ⊤ s ˆΩ ⊗ I m ∇ U (cid:16) e ⊗ q + h ˆ I s ⊗ I m ˆ γ (cid:17) , (50)ˆ ρ = ˆ P ⊤ s ˆΩ ⊗ I m ∇ g (cid:16) e ⊗ q + h ˆ I s ⊗ I m ˆ γ (cid:17) . We augment (50) by the equation (48) for λ which, by taking (49) into account, can be rewrittenas h ξ ˆ ρ ⊤ M − ˆ ρ + s − X j =1 ξ j (cid:0) ˆ ρ ⊤ j M − ˆ ρ j − − ˆ ρ ⊤ j − M − ˆ ρ j (cid:1) λ = ˆ ρ ⊤ M − (cid:16) p − hξ ˆ ψ (cid:17) − h s − X j =1 ξ j (cid:16) ˆ ρ ⊤ j M − ˆ ψ j − − ˆ ρ ⊤ j − M − ˆ ψ j (cid:17) . (51)In (50), ∇ U , when evaluated in a block vector of (block) dimension k , stands for the block vectormade up of the k vectors resulting from the corresponding application of the function. The same15traightforward notation is used for ∇ g . The new approximation is then given by, see (18) and(49), q = q + h ˆ γ , p = p − h [ ˆ ψ + ˆ ρ λ ] . (52)Note that the equations in (50), together with (52), formally coincide with those provided by aHBVM( k, s ) method applied to solve the problem defined by the Hamiltonian (3), where thevector of the multiplier λ is considered as a parameter,˙ q = M − p, ˙ p = −∇ U ( q ) − ∇ g ( q ) λ , t ≥ , q (0) = q , p (0) = p , cf. [13, 10] for details. Consequently, equation (51) defines the proper extension for handling theconstrained Hamiltonian problem (1)–(2). For this reason, we continue to refer to the numericalmethod 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, whichconverges under regularity assumptions, for all sufficiently small stepsizes h > Moreover, forseparable Hamiltonians, as it is the case in (1), the last two equations in (50) can be substitutedinto the first one, resulting in a single vector equation for ˆ γ . By settingΘ λ ( q ) := ∇ U ( q ) + ∇ g ( q ) λ , (53)we obtain ˆ γ = ˆ P ⊤ s ˆΩ ⊗ M − h e ⊗ p − h ˆ I s ˆ P ⊤ s ˆΩ ⊗ I m Θ λ (cid:16) e ⊗ q + h ˆ I s ⊗ I m ˆ γ (cid:17)i , (54)plus (51) for λ . We skip further details, since they are exactly the same as for the original HB-VMs, when applied to solve separable (unconstrained) Hamiltonian problems [13, 10].The following result follows from Theorem 7 along with standard arguments from the analysisof HBVMs [16, 10]. Theorem 8
For all sufficiently small stepsizes h > , the HBVM ( k, s ) method (51)–(54) is welldefined and symmetric. It provides a sequence of approximations ( q n , p n , λ n ) , n = 0 , , . . . , N , suchthat q n = q ( nh ) + O ( h ) ,p n = p ( nh ) + O ( h ) , ∇ g ( q n ) ⊤ M − p n = O ( h ) ,λ n = λ ( q ( nh ) , p ( nh )) + O ( h ) , Here, s is the degree of the polynomial approximation and k defines the order (actually equal to 2 k ) of thequadrature in the approximations (49). We refer to [13, 7, 10] for further details on procedures for solving the involved discrete problems. They arebased on suitable Newton-splitting procedures, already implemented in computational codes [18, 19, 20]. Clearly, ˆ ψ and ˆ ρ can be computed via the last two equations in (50), once ˆ γ and λ are known. For brevity, we omit the proof here. nd g ( q n ) = ( , if g is a polynomial of degree not larger than k/s , O ( h k ) , otherwise , (55) H ( q n , p n ) − H ( q , p ) = ( , if H is a polynomial of degree not larger than k/s , O ( h k ) , otherwise . Moreover, in case λ is constant, λ ( q ( t ) , p ( t )) ≡ ¯ λ , ∀ t ∈ [0 , T ] , the following statements hold: q n = q ( nh ) + O ( h s ) ,p n = p ( nh ) + O ( h s ) , ∇ g ( q n ) ⊤ M − p n = O ( h s ) ,λ n = ¯ λ + O ( h s ) . Remark 3
We stress that by (55), an exact or a (at least) practical conservation of both theconstraints 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 theround-off error level, in the non polynomial case. This feature of the method will be always exploitedin 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 Gausscollocation method, [12, 16, 10], which is symplectic. Moreover, in the limit k → ∞ , we retrievethe formulae studied in Section 3. This means that our approach can be also considered in theframework of Runge-Kutta methods with continuous stages [12, 26]. In this section, to illustrate the theoretical properties of HBVM( k, s ), we apply them to numericallysimulate some Hamiltonian problems of the form (1)–(2) with holonomic constraints. In the focusof our attention are properties described in Theorem 8.
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 = 1and L = 1, and normalize the gravity acceleration. Then, the Hamiltonian is given by H ( q, p ) = 12 p ⊤ p + e ⊤ q, e := (cid:18) (cid:19) , q := (cid:18) xy (cid:19) , p := ˙ q ∈ R , (56)and is subject to the constraint g ( q ) ≡ q ⊤ q − . (57)We also prescribe the initial conditions of the form q (0) = (0 , − ⊤ , p (0) = (1 , ⊤ . (58)17onsequently, the constrained Hamiltonian problem reads:¨ x = − xλ, ¨ y = − − yλ, x + y = 1 , (59) x (0) = 0 , y (0) = − , ˙ x (0) = 1 , ˙ y (0) = 0 . In order to obtain a reference solution, we rewrite the problem in polar coordinates in such a waythat θ = 0 locates the pendulum at its stable rest position, so that x = sin θ, y = − cos θ. (60)Thus, we arrive at the unconstrained Hamiltonian problem¨ θ + sin θ = 0 , θ (0) = 0 , ˙ θ (0) = 1 . (61)Once this problem is solved, the solution of (59) is recovered via the transformations (60). Moreover,the Lagrange multiplier in (57) turns out to be given by λ = 12 (cid:16) ˙ θ + cos θ (cid:17) . (62)To compute the reference solution for (59), we solve (61) by means of a HBVM(12 ,
6) method oforder 12 which is practically energy-conserving.According to (56) and (57), the Hamiltonian and the constraint are quadratic, so we expectHBVM( 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 , ,
3, the stepsizes h = 10 − − n and the intervalof integration [0 , • the solution error ( e s ), • the multiplier error ( e λ ), • the Hamiltonian error ( e H ), • the constraint error ( e g ); • the hidden constraint error, defined by e hc := max n | x n ˙ x n + y n ˙ y n | . As predicted in Theorem 8, we can see that • all methods are second-order accurate in the space variables, with HBVM(1,1) less accuratethan 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. For unconstrained Hamiltonian problems. .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 masslessrod of length L . For the conical pendulum, the initial condition is chosen in such a way that themotion is periodic with period T and occurs in the horizontal plane q = z . Again, assuming m = 1 and L = 1, and normalizing the acceleration of gravity, the Hamiltonian is H ( q, p ) := 12 p ⊤ p + e ⊤ q, e := , q := xyz , p := ˙ q ∈ R , (63)with the constraint g ( q ) := q ⊤ q − . (64)Here, we prescribe the consistent initial conditions q (0) = − − − , p (0) = − , (65)generating a periodic motion with T = 2 π, z = − − . (66)Moreover, in such a case, the multiplier λ , which has the physical meaning of the tension on therod, has to be constant and is given by λ = 2 − . (67)Note that the Hamiltonian and the constraint are quadratic and according to Theorem 8, anyHBVM( s, s ) method conserves both of them and has order 2 s . In Table 2, we list the errors in • the solution ( e s ), • the multiplier ( e λ ), • the Hamiltonian ( e H ), • the constraints ( e g ); • the hidden constraints , defined by e hc := max n k∇ g ( q n ) ⊤ M − p n k . (68)The problem is solved over 10 periods, with stepsizes h = T /n . As expected, the estimated rateof convergence for HBVM( s, s ), s = 1 , , ,
4, is 2 s . Also, the Hamiltonian and the constraint areconserved up to round-off errors. Remarkably, also the error in the multiplier ( e λ ) and in the hiddenconstraints ( e hc ) appear to be within the round-off error level, whatever stepsize is used.In Figure 1, we plot the solution error from the computation over 100 periods using HBVM(2,2)with the stepsize h = T / ≈ . Clearly, 0 > z > − L .
10 20 30 40 50 60 70 80 90 100
PERIODS E RR O R -5 Figure 1: Conical pendulum (63)–(67). Linear growth of the solution error. The solution wascomputed over 100 periods using HBVM(2,2) and the stepsize h = T / ≈ . TIME -17 -16 -15 -14 -13 -12 -11 E RR O R MULTIPLIERHIDDEN CONSTRAINT HAMILTONIANCONSTRAINT
Figure 2: Conical pendulum (63)–(67). Multiplier error (dashed line), hidden constraint error(dash-dotted line), Hamiltonian error (solid line), and constraint error (dotted line) versus time.The simulation was was carried our using HBVM(2,2) with the stepsize h = T / ≈ . .3 Modified pendulum We now consider a modified version of the previous problem. With this simulation, we aim atexploiting the conservation of energy and constraints using a suitable high-order quadrature rule(49). More precisely, we consider the following non-quadratic polynomial Hamiltonian: H ( q, p ) := 12 p ⊤ p + ( e ⊤ q ) , q, p ∈ R , (69)with the non-quadratic polynomial constraint g ( q ) := X i =1 ( e ⊤ i q ) − i ) − .
625 = 0 . (70)Here, we use the same initial condition as in (65) but the vector of the multipliers is no moreconstant, 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(3 s, s ) method will conserve both quantities. This is confirmed by Table 3, where theresults for HBVM(3,1), HBVM(6,2), and HBVM(9,3) are listed. The interval of integration wasagain [0 , • the solution ( e s ), • the multiplier ( e λ ), • the Hamiltonian ( e H ), • the constraints ( e g ); • the hidden constraints ( e hc ), formally defined via (68).Again, as predicted in Theorem 8, we see that • all methods are second-order accurate in the state variables, with HBVM(3,1) less accuratethan 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. Finally, we discuss a closed-loop rotating triangular tethered satellites system, including threesatellites (considered as mass-points) of masses m i , i = 1 , ,
3, joined by inextensible, tight, andmassless tethers, of lengths L i , i = 1 , ,
3, respectively, which orbit around a massive body. As before, for sake of simplicity, we assume unit masses and lengths, and normalize the gravity This example can be found in [22, 39]. The Earth, in the original example. q i := ( x i , y i , z i ) ⊤ ∈ R , i = 1 , ,
3, are the positions of the threesatellites, the constraints are given by g ( q ) := ( q − q ) ⊤ ( q − q ) − q − q ) ⊤ ( q − q ) − q − q ) ⊤ ( q − q ) − = 0 ∈ R , (71)and the Hamiltonian is specified by H ( q, p ) = X i =1 p ⊤ i p i − p q ⊤ i q i ! . (72)The consistent initial conditions have the form q (0) = z , q (0) = − z , q (0) = z − √ , (73)and p (0) = p (0) = , p (0) = v , (74)where z = 20 and v is such that the initial Hamiltonian is zero. This provides a configuration inwhich the first two satellites remain parallel to each other, moving in the planes y = and y = − ,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) methodwith the stepsize h = 0 . steps, we obtain a qualitatively correct solution which conservesthe Hamiltonian and the constraints within the round-off error level, see Figure 3. Here, we alsoplot the hidden constraints errors k∇ g ( q n ) ⊤ M − p n k . At last, in Table 4, we list the followingerrors arising when solving the problem with HBVM(6 , s ) methods for s = 1 , , h = 10 − − n , over the interval [0 , • the solution error ( e s ), • the multipliers error ( e λ ), • the Hamiltonian error ( e H ), • the constraints error ( e g ); • the hidden constraints errors ( e hc ), formally defined by (68).Again, as shown in Theorem 8, • all methods are second-order accurate in the state variables, with HBVM(6,1) much lessaccurate 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;22 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 thenumerical solution of the Hamiltonian problems with holonomic constraints can be recommended,although the method is only second-order accurate.
In this paper, we have considered the numerical solution of Hamiltonian problems with holonomicconstraints, by resorting to a line-integral formulation of the conservation of the constraints. Thisapproach enables to derive an expression for the Lagrange multipliers in which second derivativesare not used . From the discretization of the resulting formulae, we have obtained a suitable variantof the Hamiltonian Boundary Value Methods (HBVMs), formerly designed as an energy-conservingRunge-Kutta methods for unconstrained Hamiltonian problems. Numerical experiments support-ing the theoretical findings are enclosed.This paper has been initiated in spring 2017, during the visit of the first author at the Institutefor Analysis and Scientific Computing, Vienna University of Technology, Vienna, Austria.
References [1] L. Barletti, L. Brugnano, G. Frasca Caccia, F. Iavernaro. Energy-conserving methods for thenonlinear Schr¨odinger equation.
Appl. Math.Comput. (2018) 3–18.[2] P. Amodio, L. Brugnano, F. Iavernaro. Energy-conserving methods for Hamiltonian BoundaryValue Problems and applications in astrodynamics.
Adv. Comput. Math. (2015) 881–905.[3] H.C. Andersen. Rattle: a “velocity” version of the Shake algorithm for molecular dynamicscalculations. J. Comp. Phys. (1983) 24–34.[4] G. Benettin, A.M. Cherubini, F. Fass`o. A changing chart symplectic algorithm for rigid bodiesand other Hamiltonian systems on manifolds. SIAM J. Sci. Comput. , No. 4 (2001) 1189–1203.[5] K.E. Brenan, S.L. Campbell, L.R. Petzold. Numerical solution of initial-value problems indifferential-algebraic equations . SIAM, Philadelphia, PA, 1996.[6] L. Brugnano, M.Calvo, J.I.Montijano, L.R`andez. Energy preserving methods for Poisson sys-tems.
J. Comput. Appl. Math. (2012) 3890–3904.[7] L. Brugnano, G. Frasca Caccia, F. Iavernaro. Efficient implementation of Gauss collocation andHamiltonian Boundary Value Methods.
Numer. Algorithms (2014) 633–650.[8] L. Brugnano, G. Frasca Caccia, F. Iavernaro. Energy conservation issues in the numerical solu-tion of the semilinear wave equation. Appl. Math.Comput. (2015) 842–870.[9] L. Brugnano, F. Iavernaro. Line Integral Methods which preserve all invariants of conservativeproblems.
J. Comput. Appl. Math. (2012) 3905–3919.23
100 200 300 400 500 600 700 800 900 100010 -10 -5 -16 -16 -14 -14 Figure 3: Tethered satellites system (71)–(74). The problem was solved using HBVM(6,2) withthe stepsize h = 0 . steps. From top to bottom: hidden constraints error (first plot),Hamiltonian error (second plot), constraints errors (third to fifth plots).2410] L. Brugnano, F. Iavernaro. Line Integral Methods for Conservative Problems . Chapman etHall/CRC, Boca Raton, FL, 2016.[11] L. Brugnano, F. Iavernaro, D. Trigiante. Hamiltonian BVMs (HBVMs): a family of “drift-free”methods for integrating polynomial Hamiltonian systems.
AIP Conf. Proc. (2009) 715–718.[12] L. Brugnano, F. Iavernaro, D. Trigiante. Hamiltonian Boundary Value Methods (Energy Pre-serving Discrete Line Integral Methods).
JNAIAM. J. Numer. Anal. Ind. Appl. Math. , No. 1-2(2010), 17–37.[13] L. Brugnano, F. Iavernaro, D. Trigiante. A note on the efficient implementation of HamiltonianBVMs. J. Comput. Appl. Math. (2011) 375–383.[14] L. Brugnano, F. Iavernaro, D. Trigiante. A two-step, fourth-order method with energy preserv-ing properties.
Comput. Phys. Commun. (2012) 1860–1868.[15] L. Brugnano, F. Iavernaro, D. Trigiante. Energy and QUadratic Invariants Preserving integra-tors based upon Gauss collocation formulae.
SIAM J. Numer. Anal. , No. 6 (2012) 2897–2916.[16] L. Brugnano, F. Iavernaro, D. Trigiante. A simple framework for the derivation and analysis ofeffective one-step methods for ODEs. Appl. Math.Comput. (2012) 8475–8485.[17] L. Brugnano, F. Iavernaro, D. Trigiante. Analisys of Hamiltonian Boundary Value Methods(HBVMs): a class of energy-preserving Runge-Kutta methods for the numerical solution ofpolynomial Hamiltonian systems.
Commun. Nonlinear Sci. Numer. Simul. (2015) 650–667.[18] L. Brugnano, C. Magherini. Recent advances in linear analysis of convergence for splittings forsolving ODE problems. Appl. Numer. Math. (2009) 542–557.[19] L. Brugnano, C. Magherini. The BiM code for the numerical solution of ODEs. J. Comput.Appl. Math. (2004) 145–158.[20] L. Brugnano, C. Magherini. Blended Implicit Methods for solving ODE and DAE problems,and their extension for second order problems.
J. Comput. Appl. Math. (2007) 777–790.[21] L. Brugnano, Y. Sun. Multiple invariants conserving Runge-Kutta type methods for Hamilto-nian problems.
Numer. Algorithms (2014) 611–632.[22] Z.Q. Cai, X.F. Li, X. F., H. Zhou. Nonlinear dynamics of a rotating triangular tethered satelliteformation near libration points. Aerospace Science and Technology (2015) 384–391.[23] P. Console, E. Hairer, C. Lubich. Symmetric multistep methods for constrained Hamiltoniansystems. Numer. Math. (2013) 517–539.[24] L. Jay. Symplectic partitioned Runge-Kutta methods for constrained Hamiltonian systems.
SIAM J. Numer Anal. , No. 1 (1996) 368–387.[25] E. Hairer. Global modified Hamiltonian for constrained symplectic integrators. Numer. Math. (2003) 325–336. 2526] E. Hairer. Energy preserving variant of collocation methods. JNAIAM. J. Numer. Anal. Ind.Appl. Math. , No. 1-2 (2010), 73–84.[27] E. Hairer, C. Lubich, G. Wanner. Geometric Numerical Integration, Structure-Preserving Al-gorithms for Ordinary Differential Equations, Second Edition. Springer-Verlag, Berlin, 2006.[28] E. Hairer, G. Wanner. Solving Ordinary Differential Equations II, Stiff and Differential-A lge-braic Problems, Second Revised Edition . Springer-Verlag, Berlin, 2010.[29] B. Leimkhuler, S. Reich. Symplectic integration of constrained Hamiltonian systems.
Math.Comp. , No. 208 (1994) 589–605.[30] B. Leimkhuler, S. Reich. Simulating Hamiltonian Dynamics . Cambridge University Press, 2004.[31] B. Leimkhuler, R.D. Skeel. Symplectic numerical integrators in constrained Hamiltonian sys-tems.
J. Comput. Phys. (1994) 117–125.[32] S. Leyendecker, P. Betsch, P. Steinmann. Energy-conserving integration of constrained Hamil-tonian systems – a comparison of approaches.
Computational Mechanics (2004) 174–185.[33] O. Gonzalez. Mechanical systems subject to holonomic constraints: Differential–algebraic for-mulations and conservative integration. Physica D (1999) 165–174.[34] S. Reich. Symplectic integration of constrained Hamiltonian systems by composition methods.
SIAM J. Numer. Anal. No. 2 (1996) 475–491.[35] S. Reich. On higher-order semi-explicit symplectic partitioned Runge-Kutta methods for con-strained Hamiltonian systems.
Numer. Math. (1997) 231–247.[36] W.M. Seiler. Position versus momentum projections for constrained Hamiltonian systems. Nu-mer. Algorithms (1998) 223–234.[37] W.M. Seiler. Numerical integration of constrained Hamiltonian systems using Dirac brackets. Math. Comp. , No. ˙226 (1999) 661–681.[38] J.P. Ryckaert, G. Ciccotti, H.J.C. Berendsen. Numerical integration of the cartesian equationsof motion of a system with constraints: molecular dynamics of n -alkanes. J. Comput. Phys. (1977) 327–341.[39] Y. Wei, Z. Deng, Q. Li, B. Wang. Projected Runge-Kutta methods for constrained Hamiltoniansystems. Appl. Math. Mech. – Engl. Ed. , No. 8 (2016) 1077–1094.26able 1: Planar pendulum (56)–(57). Errors from HBVM( s, s ) method for s = 1 , ,
3, when solvingthe problem over the interval [0 ,
10] with stepsizes h = 10 − − n . s = 1 n e s rate e λ rate e H e g e hc rate0 2.5700e-02 – 3.4253e-02 – 5.5511e-17 9.9920e-16 2.3487e-03 –1 6.4260e-03 2.00 1.7386e-02 0.98 1.1102e-16 6.6613e-16 5.8639e-04 2.002 1.6070e-03 2.00 8.7406e-03 0.99 1.1102e-16 5.5511e-16 1.4654e-04 2.003 4.0181e-04 2.00 4.3835e-03 1.00 1.1102e-16 1.9984e-15 3.6633e-05 2.004 1.0045e-04 2.00 2.1948e-03 1.00 1.1102e-16 1.8874e-15 9.1580e-06 2.005 2.5114e-05 2.00 1.0982e-03 1.00 1.1102e-16 2.4425e-15 2.2895e-06 2.006 6.2785e-06 2.00 5.4929e-04 1.00 1.1102e-16 3.4417e-15 5.7238e-07 2.007 1.5696e-06 2.00 2.7470e-04 1.00 1.1102e-16 6.1062e-15 1.4311e-07 2.008 3.9248e-07 2.00 1.3743e-04 1.00 1.1102e-16 8.5487e-15 3.5902e-08 2.00 s = 2 n e s rate e λ rate e H e g e hc rate0 1.6695e-03 – 3.5176e-02 – 1.1102e-16 8.8818e-16 2.3539e-03 –1 4.1412e-04 2.01 1.7585e-02 1.00 1.1102e-16 8.8818e-16 5.8670e-04 2.002 1.0332e-04 2.00 8.7919e-03 1.00 1.1102e-16 1.1102e-15 1.4656e-04 2.003 2.5816e-05 2.00 4.3958e-03 1.00 1.1102e-16 9.9920e-16 3.6634e-05 2.004 6.4533e-06 2.00 2.1979e-03 1.00 1.1102e-16 1.8874e-15 9.1581e-06 2.005 1.6133e-06 2.00 1.0990e-03 1.00 1.1102e-16 2.7756e-15 2.2895e-06 2.006 4.0332e-07 2.00 5.4948e-04 1.00 1.1102e-16 3.6637e-15 5.7238e-07 2.007 1.0086e-07 2.00 2.7477e-04 1.00 1.1102e-16 5.5511e-15 1.4314e-07 2.008 2.5286e-08 2.00 1.3751e-04 1.00 1.1102e-16 1.0547e-14 3.5884e-08 2.00 s = 3 n e s rate e λ rate e H e g e hc rate0 1.6658e-03 – 3.5178e-02 – 1.1102e-16 1.1102e-15 2.3539e-03 –1 4.1386e-04 2.01 1.7585e-02 1.00 1.1102e-16 8.8818e-16 5.8670e-04 2.002 1.0331e-04 2.00 8.7919e-03 1.00 1.1102e-16 7.7716e-16 1.4656e-04 2.003 2.5815e-05 2.00 4.3958e-03 1.00 1.1102e-16 8.8818e-16 3.6634e-05 2.004 6.4532e-06 2.00 2.1979e-03 1.00 1.1102e-16 1.6653e-15 9.1581e-06 2.005 1.6133e-06 2.00 1.0990e-03 1.00 1.1102e-16 2.9976e-15 2.2895e-06 2.006 4.0331e-07 2.00 5.4948e-04 1.00 1.1102e-16 3.4417e-15 5.7238e-07 2.007 1.0086e-07 2.00 2.7477e-04 1.00 1.1102e-16 5.2180e-15 1.4314e-07 2.008 2.5220e-08 2.00 1.3739e-04 1.00 1.1102e-16 8.9928e-15 3.5791e-08 2.0027able 2: Conical pendulum (63)–(67). Errors from HBVM( s, s ) method for s = 1 , , ,
4, whensolving the problem over 10 periods with stepsizes h = T /n . s = 1 n e s rate e λ e H e g e hc
10 1.1543e 00 – 6.5503e-15 1.1102e-16 1.5543e-15 6.9435e-1520 4.0996e-01 1.49 1.2212e-14 1.1102e-16 6.6613e-16 7.3344e-1530 1.9021e-01 1.89 7.2831e-14 1.1102e-16 8.8818e-16 2.6870e-1440 1.0794e-01 1.97 2.3959e-13 1.1102e-16 6.6613e-16 6.8291e-1450 6.9285e-02 1.99 2.6112e-13 1.1102e-16 4.4409e-16 5.5241e-1460 4.8178e-02 1.99 1.8818e-13 1.1102e-16 1.2212e-15 4.0510e-1470 3.5420e-02 2.00 6.1351e-13 1.1102e-16 6.6613e-16 9.4355e-1480 2.7130e-02 2.00 8.1246e-13 1.1102e-16 5.5511e-16 1.0761e-1390 2.1441e-02 2.00 7.5329e-13 1.1102e-16 5.5511e-16 8.8662e-14100 1.7371e-02 2.00 1.4311e-12 1.1102e-16 5.5511e-16 1.5112e-13 s = 2 n e s rate e λ e H e g e hc
10 1.1168e-02 – 6.5503e-15 1.1102e-16 6.6613e-16 7.7539e-1520 7.1061e-04 3.97 1.8208e-14 5.5511e-17 3.3307e-16 1.1374e-1430 1.4083e-04 3.99 6.4060e-14 1.1102e-16 3.3307e-16 2.2125e-1440 4.4610e-05 4.00 2.2171e-13 1.1102e-16 4.4409e-16 6.2087e-1450 1.8282e-05 4.00 9.2704e-14 1.1102e-16 4.4409e-16 2.1613e-1460 8.8190e-06 4.00 3.6859e-13 1.1102e-16 4.4409e-16 6.4763e-1470 4.7611e-06 4.00 8.5076e-13 1.1102e-16 4.4409e-16 1.3128e-1380 2.7912e-06 4.00 1.2552e-12 1.1102e-16 2.2204e-16 1.6921e-13 s = 3 n e s rate e λ e H e g e hc
10 3.1758e-05 – 6.1062e-15 1.1102e-16 1.5543e-15 9.8364e-1520 5.0199e-07 5.98 1.7319e-14 1.1102e-16 6.6613e-16 1.0819e-1430 4.4164e-08 5.99 9.8921e-14 1.1102e-16 2.2204e-16 3.4922e-1440 7.8663e-09 6.00 2.0650e-13 1.1102e-16 4.4409e-16 5.4473e-1450 2.0628e-09 6.00 2.3292e-13 1.1102e-16 4.4409e-16 4.9311e-1460 6.9103e-10 6.00 4.7151e-13 1.1102e-16 4.4409e-16 8.3119e-14 s = 4 n e s rate e λ e H e g e hc
10 4.9944e-08 – 1.1768e-14 1.1102e-16 1.5543e-15 1.5603e-1420 1.9676e-10 7.99 1.7431e-14 1.1102e-16 6.6613e-16 1.1910e-1430 7.6630e-12 8.00 5.7399e-14 1.1102e-16 3.3307e-16 2.1564e-1440 7.3944e-13 8.13 4.1411e-14 1.1102e-16 4.4409e-16 1.0999e-1428able 3: Modified pendulum (69)–(70). Errors from HBVM(3 s, s ) method for s = 1 , ,
3, whensolving the problem over the interval [0 ,
10] with stepsizes h = 10 − − n . s = 1 n e s rate e λ rate e H e g e hc rate0 2.0539e-02 – 1.0864e-01 – 1.1102e-16 1.6431e-14 1.5279e-02 –1 4.9675e-03 2.05 6.4279e-02 0.76 2.2204e-16 7.5495e-15 3.9290e-03 1.962 1.2365e-03 2.01 3.5621e-02 0.85 2.2204e-16 3.6637e-15 9.7072e-04 2.023 3.0878e-04 2.00 1.8727e-02 0.93 2.2204e-16 2.1094e-15 2.4193e-04 2.004 7.7173e-05 2.00 9.5965e-03 0.96 2.2204e-16 8.8818e-16 6.0436e-05 2.005 1.9292e-05 2.00 4.8569e-03 0.98 2.2204e-16 5.5511e-16 1.5106e-05 2.006 4.8229e-06 2.00 2.4432e-03 0.99 2.2204e-16 6.6613e-16 3.7764e-06 2.007 1.2057e-06 2.00 1.2251e-03 1.00 2.2204e-16 6.6613e-16 9.4417e-07 2.008 3.0139e-07 2.00 6.1472e-04 1.00 2.2204e-16 6.6613e-16 2.3608e-07 2.00 s = 2 n e s rate e λ rate e H e g e hc rate0 6.0495e-03 – 1.5224e-01 – 1.1102e-16 5.7732e-15 1.7516e-02 –1 1.4027e-03 2.11 7.6627e-02 0.99 2.2204e-16 3.8858e-15 4.6710e-03 1.912 3.4600e-04 2.02 3.8806e-02 0.98 1.1102e-16 2.3315e-15 1.1666e-03 2.003 8.6197e-05 2.01 1.9532e-02 0.99 2.2204e-16 1.6653e-15 2.9091e-04 2.004 2.1530e-05 2.00 9.7988e-03 1.00 2.2204e-16 5.5511e-16 7.2716e-05 2.005 5.3813e-06 2.00 4.9076e-03 1.00 2.2204e-16 5.5511e-16 1.8175e-05 2.006 1.3453e-06 2.00 2.4559e-03 1.00 2.2204e-16 6.6613e-16 4.5440e-06 2.007 3.3632e-07 2.00 1.2285e-03 1.00 2.2204e-16 6.6613e-16 1.1360e-06 2.008 8.4002e-08 2.00 6.1386e-04 1.00 2.2204e-16 6.6613e-16 2.8414e-07 2.00 s = 3 n e s rate e λ rate e H e g e hc rate0 6.0698e-03 – 1.5231e-01 – 2.2204e-16 4.6629e-15 1.7532e-02 –1 1.4040e-03 2.11 7.6632e-02 0.99 1.1102e-16 3.9968e-15 4.6715e-03 1.912 3.4608e-04 2.02 3.8806e-02 0.98 1.1102e-16 1.7764e-15 1.1666e-03 2.003 8.6202e-05 2.01 1.9532e-02 0.99 2.2204e-16 1.3323e-15 2.9091e-04 2.004 2.1530e-05 2.00 9.7988e-03 1.00 2.2204e-16 5.5511e-16 7.2716e-05 2.005 5.3814e-06 2.00 4.9076e-03 1.00 2.2204e-16 6.6613e-16 1.8175e-05 2.006 1.3453e-06 2.00 2.4560e-03 1.00 2.2204e-16 6.6613e-16 4.5439e-06 2.007 3.3623e-07 2.00 1.2283e-03 1.00 2.2204e-16 6.6613e-16 1.1361e-06 2.008 8.4116e-08 2.00 6.1452e-04 1.00 2.2204e-16 6.6613e-16 2.8410e-07 2.0029able 4: Tethered satellite system (71)–(74). Errors from the HBVM(6 , s ) method for s = 1 , , h = 10 − − n . s = 1 n e s rate e λ rate e H e g e hc rate0 9.2893e-04 – 2.1218e-06 – 5.5511e-17 1.2212e-14 9.6503e-07 –1 2.3234e-04 2.00 1.0689e-06 0.99 5.5511e-17 1.2212e-14 2.4108e-07 2.002 5.8093e-05 2.00 5.3623e-07 1.00 6.9389e-17 1.5765e-14 6.0272e-08 2.003 1.4524e-05 2.00 2.6834e-07 1.00 6.9389e-17 1.5099e-14 1.5090e-08 2.00 s = 2 n e s rate e λ rate e H e g e hc rate0 1.8586e-07 – 2.1635e-06 – 6.2450e-17 1.1990e-14 1.3053e-06 –1 3.9859e-08 2.22 1.0812e-06 1.00 4.8572e-17 1.4877e-14 3.2584e-07 2.002 9.5501e-09 2.06 5.4039e-07 1.00 6.9389e-17 1.4433e-14 8.1466e-08 2.003 2.3636e-09 2.01 2.7052e-07 1.00 6.9389e-17 1.4655e-14 2.0374e-08 2.00 s = 3 n e s rate e λ rate e H e g e hchc