Structure-preserving integrators based on a new variational principle for constrained mechanical systems

A new variational principle for mechanical systems subject to holonomic constraints is presented. The newly proposed GGL principle is closely related to the often used Gear-Gupta-Leimkuhler (GGL) stabilization of the differential–algebraic equations governing the motion of constrained mechanical systems. The GGL variational principle relies on an extension of the Livens principle (sometimes also referred to as Hamilton–Pontryagin principle) to mechanical systems subject to holonomic constraints. In contrast to the original GGL stabilization, the new approach facilitates the design of structure-preserving integrators. In particular, new variational integrators are presented, which result from the direct discretization of the GGL variational principle. These variational integrators are symplectic and conserve momentum maps in the case of systems with symmetry. In addition to that, a new energy–momentum scheme is developed, which results from the discretization of the Euler–Lagrange equations pertaining to the GGL variational principle. The numerical properties of the newly devised schemes are investigated in representative examples of constrained mechanical systems.


Introduction
The motion of mechanical systems subject to holonomic constraints is governed by differential-algebraic equations (DAEs) with differentiation index three (cf. [1,2]). These DAEs are known to be prone to numerical instabilities. The GGL stabilization due to Gear, Leimkuhler and Gupta [3] yields an index reduction by taking into account the hidden constraints on velocity level. In particular, the GGL stabilization achieves an index reduction by minimal extension [1] and facilitates a stable numerical integration of the resulting index 2 DAEs.
Since, apart from the index reduction property, the GGL stabilization does not suffer from the drift off phenomenon [2], it gained high popularity over the last decades. In fact, the GGL stabilization has been widely used in the design of numerical methods for constrained mechanical systems. For example, the GGL stabilization is frequently used in the computer simulation of flexible multibody dynamics [2,4,5]. It has been shown to be advantageous for the numerical solution of optimal control problems of multibody systems [6,7]. Moreover, it has been employed in the design of DAE Lie group methods [8]. Related works on constrained mechanical systems benefiting from the usage of velocity constraints have been developed in [9,10].
In the past, a number of structure-preserving schemes has been developed for the numerical integration of constrained mechanical systems. Representative examples are symplectic integrators [11][12][13] and energymomentum schemes [14][15][16][17]. These methods typically rely on a direct discretization of the constrained equations of motion in the Hamiltonian framework. In the present work, we follow another path towards the design of structure-preserving integrators, which is based on Livens principle. This variational principle can be traced back to Livens [18] (see also Pars [19]), where it has been formulated in the context of mechanical systems without constraints. In recent times, Livens principle has also been termed Hamilton-Pontryagin principle (see [20][21][22]) due to its close relation with optimal control problems.
In the present work, we propose an extension of Livens principle to mechanical systems subject to holonomic constraints. The Euler-Lagrange equations of the newly devised variational principle yield an index reduction of the underlying DAEs to index 2 in the spirit of the GGL stabilization. Therefore, we choose the name GGL principle for the novel variational principle, which accounts for both position constraints and hidden velocity constraints. In contrast to the original GGL stabilization, the availability of the GGL variational principle makes possible the design of variational integrators (VIs). VIs rely on the direct discretization of the underlying variational functional [23][24][25]. We stress that previously developed VIs for constrained mechanical systems such as [23] do not rely on GGL-type stabilization, which hinges on a suitable modification of the underlying continuous formulation. VIs are typically symplectic and capable of conserving momentum maps in the case of systems with symmetry. In addition to that, we show that the Euler-Lagrange equations of the GGL principle facilitate the design of energy-momentum (EM) schemes. EM schemes typically rely on the notion of a discrete derivative [26][27][28] and are capable of conserving the system's total energy and momentum maps.
The remainder of this work is structured as follows. Section 2 contains an outline of Livens principle and a summary of the original GGL stabilization approach to constrained mechanical systems. In addition to that, structural properties of interest are summarized. In Sect. 3, the newly proposed GGL variational principle for constrained mechanical systems is introduced. Novel VIs based on the GGL principle are devised in Sect. 4. In addition to that, a new EM scheme resulting from the discretization of the Euler-Lagrange equations of the GGL principle is presented in Sect. 5. In Sect. 6, the properties of the time-stepping schemes at hand are investigated in the context of representative numerical examples. Eventually, conclusions are drawn in Sect. 7.

Fundamentals
This section contains an outline of important relationships required for the subsequent developments.

Livens principle
Consider a dynamical system with d degrees of freedom and generalized coordinates q ∈Q. Within the time interval I = [0, T ], the trajectory on the configuration spaceQ is characterized by q : I →Q.
The corresponding velocities v ∈ T qQ can be used as independent variables in Hamilton's principle of least action. To this end, the kinematic relationq = v is appended to the Lagrangian by using Lagrange multipliers p ∈ T * qQ . The corresponding augmented action integral reads where L(q, v) : TQ → R is the Lagrangian. The functional (1) can be linked to Livens theorem (cf. Section 26.2 in Pars [19]) which goes back to Livens [18]. More recently, the term Hamilton-Pontryagin principle has been used [20][21][22] in allusion to the close relation with optimization problems for dynamical systems (cf. Bryson and Ho [29]). 1 Livens principle unifies both Lagrangian and Hamiltonian viewpoints on mechanics and automatically accounts for the Legendre transformation. By stating the stationary condition δS(q, v, p) = 0 and performing the variations with respect to every independent variable, one obtains the Euler-Lagrange equations in the forṁ Note that integration by parts and the usual endpoint conditions δq(0) = δq(T ) = 0 have been employed. Moreover, D j L(q, v) stands for the partial derivative of function L with respect to the j's argument.
With regard to (2c), the multiplier p can be identified as the conjugate momentum. Accordingly, the Legendre transformation is built into the variational principle via its fiber derivative. Note that after substituting (2c) for p in (2b) and making use of (2a), Livens principle traces back to the Lagrange equations.
In the present work, we focus on mechanical systems whose Lagrangian takes the form Here, M ∈ R d×d is a symmetric, positive-definite mass matrix and V :Q → R is a potential function. Introducing the total energy function variational functional (1) can be recast in the form Remark 1 Making use of (2c), the velocities can be expressed in terms of the coordinates q and conjugate momenta p, i.e., v =ṽ(q, p). Thus, it is possible to eliminate the velocities from the formulation such that the generalized energy (4) can be identified as the Hamiltonian In the case of mechanical systems with Lagrangian of the form (3), one obtains Substituting (6) into (5), the two-field variational principlẽ remains, which is sometimes referred to as the modified Hamilton's principle (e.g., Sect. 8-5 in Goldstein [31]). The Euler-Lagrange equations corresponding to the modified Hamilton's principle are given bẏ and thus take the form of Hamilton's canonical equations of motion.

Structural properties
Next, we outline main structural properties of mechanical systems in the context of the equations of motion (2) emanating from Livens principle. These properties will also play a crucial role in the subsequent treatment of constrained mechanical systems.

Symmetries and momentum maps
Consider a one-parameter family of curves q α (t) in configuration spaceQ with q 0 (t) = q(t). Conforming with (2a), we introduce v α (t) =q α (t). The corresponding infinitesimal generators are defined by The mechanical system has symmetry if at all times t the condition holds. Correspondingly, the infinitesimal symmetry condition is given by According to Noether's theorem, each symmetry of the system leads to the conservation of an associated momentum map. This can be verified by resorting to the equations of motion (2). Scalar multiplying (2b) by ξ Q and (2c) by ξ V and subsequently adding both resulting equations yielḋ where symmetry condition (12) has been used. Since (2a) together with (10) implies ξ V =ξ Q , (13) gives rise tȯ so that the momentum map is a conserved quantity. We are particularly interested on curves resulting from the action of a matrix group G on the configuration spaceQ, where A α ∈ G ⊂ GL(d, R), the general linear group of R d . As before, q 0 (t) = q(t), which implies A 0 = I, where I is the d × d identity matrix. Now, the infinitesimal generator (10a) assumes the specific form where ξ ∈ g, the Lie algebra of G. With regard to (2a), we further obtain together with Moreover, (2c) implies p α = D 2 L(q α , v α ). Thus, making use of symmetry property (11), one obtains from which follows the corresponding infinitesimal generator We eventually remark that the symmetry properties considered above imply the invariance of the energy function (4). That is, E(q α , v α , p α ) = E(q, v, p).

Conservation of energy
In the case of an autonomous Lagrangian, the total energy of the mechanical system is conserved. Differentiating energy function (4) with respect to time yields Taking into account the Euler-Lagrange equations (2), one obtainsĖ = 0, so that the total energy is a conserved quantity.

Symplecticness
The flow generated by the equations of motion (2) is symplectic. It is convenient to define symplecticness by expressing the symplectic two-form Ω in terms of the exterior or wedge product (Sanz-Serna and Calvo [32], Leimkuhler and Reich [13]). Accordingly, where dq and dp represent differentials of the coordinates q and momenta p, respectively. Important properties of the wedge product are summarized in Appendix A.1. Correspondingly, the property of symplecticness is equivalent to the fact that the symplectic two-form Ω is conserved along solution trajectories of the equations of motion (2), i.e.,Ω = 0. For completeness, this is verified in Appendix A.2.

Constrained dynamics and classical GGL stabilization
Assume that the coordinates q are redundant due to the presence of m independent scleronomic, holonomic constraints. Arranging the constraint functions g k : R d → R (k = 1, . . . , m) in a column vector g ∈ R m , the constraints are given by The constraint functions are assumed to be independent, so that the constraint Jacobian G(q) = Dg(q) has rank m. Note that the constraints give rise to the configuration manifold To get the equations of motion, Livens principle can be augmented to account for the constraints (24). Accordingly, the corresponding functional assumes the form where λ ∈ R m represents Lagrange multipliers for the enforcement of the constraints. Imposing the stationary conditions on functional (26), the resulting Euler-Lagrange equations are given bẏ If we take into account the specific form of the Lagrangian (3), the above equations can be recast in the forṁ Note that (28c) can be used to eliminate either the velocity v or the momentum p from the formulation. It is well known that the above differential algebraic equations (DAEs) have differentiation index 3 (see, for example, [1,2]). One established way to stabilize numerical integration schemes for constrained mechanical systems is to apply the GGL method [3]. The GGL method also takes into account the secondary (or hidden) constraints on velocity level arising from differentiating the constraints (24) with respect to time and making use of (28a) and (28c). Accordingly, the secondary constraints can be written as g v (q, p) = 0, where The GGL method relies on the inclusion of the secondary constraints as additional algebraic constraints. This procedure leads to the DAEṡ The new variables γ ∈ R m can be viewed as additional multipliers for the enforcement of the secondary constraints (30d). It is well known [1][2][3] that the DAEs (30) have differentiation index 2. It is important to realize that the inclusion of both the constraints (30c) and (30d) confines the dynamics of the constrained mechanical system to the phase space manifold This is a prerequisite for the design of symplectic integrators which preserve the restriction to M of the symplectic two-form (23) in R d × R d , see Leimkuhler and Skeel [11]. However, the presence of the extra variable γ in (30a) essentially hinders the design of structurepreserving integrators. For example, concerning conservation of energy, differentiating the Hamiltonian (7) with respect to time yields Accordingly, energy conservation holds under the provision that γ = 0, which is true for the time-continuous case [1][2][3]. However, in the discrete setting the discrete counterpart of γ is in general non-zero, since it is required to impose the velocity constraints. Accordingly, the presence of the extra γ term in (32) indicates that the design of energy conserving time-stepping schemes is not feasible. Similar conclusions can be made regarding the conservation of the symplectic twoform (see Appendix A.3 for details) as well as the conservation of momentum maps.

GGL principle
In this section, we present a new variational principle which achieves an index reduction of the index 3 DAEs (27) in the spirit of the GGL method. We show that the newly introduced GGL variational principle is particularly well suited for the design of structure-preserving numerical methods.

Governing equations
The GGL variational principle can be viewed as an extension of Livens principle (1). Similar to the GGL stabilization, the GGL principle relies on the introduction of extra Lagrange multipliers γ ∈ R m which serve the purpose to enforce the secondary constraints (29) in addition to the primary constraints (24). In particular, the GGL principle relies on the augmented action integral where the velocity constraint function g v has previously been introduced in (29). It can be easily observed that functional (33) boils down to Livens augmented action integral (1) if no constraints are present. Imposing the stationary condition With regard to (35b), we apply integration by parts to obtain Consequently, taking into account the arbitrariness of the variations δq, δv, δp, δλ and δγ , as well as considering the usual endpoint conditions δq(0) = δq(T ) = 0, we obtain the Euler-Lagrange equationṡ It can be concluded from the fiber derivative (37c) that p represents the conjugate momenta in analogy to the Livens principle. Moreover, it can be observed from (37d) and (37e) that both the primary and the secondary constraints are enforced. The presence of the secondary constraints entails an index reduction in the spirit of the GGL stabilization. In particular, we show in Appendix A.4 that the DAEs (37) (i) have differentiation index 2, and (ii) yield γ = 0 in analogy to the classical GGL stabilization. Thus, in the timecontinuous setting, the DAEs (37) boil down to the underlying DAEs (27).
Remark 2 Similar to the original Livens principle, p in the GGL principle (33) plays the role of a Lagrange multiplier for the imposition of the side conditioṅ where cf. (37a).

Remark 4
There seems to be some ambiguity in the representation of the secondary constraints in functional (33). An alternative form relies oñ where the secondary constraints are given by g v (q, v) = G(q)v. This version, however, does not yield the correct equations of motion for constrained mechanical systems. For further details, we refer to Appendix A.5.

Structural properties
This section contains an outline of the main structural properties of the Euler-Lagrange equations (37) emanating from the GGL principle. These properties are not affected by the additional multiplier γ as will become obvious in the sequel.

Symmetries and momentum maps
Assume that the Lagrangian has invariance property (11), i.e., L(q α , v α ) = L(q, v), where the oneparameter curves (q α , v α ) have to comply with the constraints. That is, q α ∈ Q and v α ∈ T q Q. Correspondingly, (q α , p α ) ∈ M. Now the symmetry of the constrained mechanical system at hand is assumed to be characterized by the infinitesimal symmetry conditions Note that (41a) follows from invariance property (11) as has been shown in (12). Similarly, (41b) and (41c) typically follow from invariance properties g v (q α , p α ) = g v (q, p) and g(q α ) = g(q), respectively. The symmetry properties (41) give rise to the conservation of an associated momentum map as shown next.
Scalar multiplying (37a) by (−ξ P ), (37b) by ξ Q , (37c) by ξ V , and subsequently adding the resulting equations yields where symmetry conditions (41) have been taken into account. Inserting the infinitesimal generators (17), (19) and (21) into the last equation yields and thus Accordingly, the momentum map is a conserved quantity.

Conservation of energy
Similar to Sect. 2.2.2, we consider the time derivative of the total energy (4) to obtain The last term on the right-hand side vanishes due to (37c). Inserting from (37a) for v and from (37b) for D 1 L into the last equation yields The last equation holds independently from the Lagrange multipliers λ and γ due to the consistency conditions dg/ dt = 0 and dg v / dt = 0 which result from (37d) and (37e), respectively.
To show this, consider the differential one-forms resulting from (37): These relationships can be inserted into the expression for the time derivative of the symplectic two-form (23) leading to Missing terms in the last equation vanish (cf. Appendix A.2) and thus have been omitted. The second term on the right-hand side of the last equation vanishes as well as shown in Appendix A.3, (135). Accordingly, we arrive at d dt Note that the first term on the right-hand side of the above equation vanishes since (37e) implies dg v (q, p) = 0. The second term also vanishes due to the equality of mixed derivatives. Accordingly, the symplectic two-form is preserved, independently from multiplier γ .

GGL variational integrators
Variational integrators (VIs) can be designed by directly discretizing the underlying functional. Subsequently, the discretized functional can be extremized to obtain the discrete Euler-Lagrange (DEL) equations. In this section, this procedure shall be applied in the framework of the newly devised GGL functional. We refer to Lew and Mata [24] and Hairer et al. [33] for more background on variational integrators.

First-order scheme
To illustrate the design of VIs emanating from the GGL functional, we start with a simple scheme, which has already been presented in the preliminary work [34].
with t 0 = 0 and t N = T . For simplicity, we consider time steps of constant size h = t n+1 − t n . The approximation of a function (•)(t) at time t n is denoted by (•) n , so that, for example, q n ≈ q(t n ). The GGL functional (33) may now be approximated by the discrete action sum has been introduced. Note that in each time interval [t n t n+1 ], side condition (38) is discretized similarly to the explicit Euler method with associated Lagrange multiplier p n+1 . Moreover, t n+1 t n L(q, v)dt is approximated by hL(q n , v n ), while multiplier λ n+1 enforces the holonomic contraints at the end of each time interval.
Stationary conditions can be applied directly to the discrete functional. Taking into account the independence of the respective variables, this yields Next, we apply an index shift in the second part of (54b) from n + 1 to n. This procedure is also referred to as discrete integration by parts (see Marsden and West [23]). In addition to that, we apply the usual endpoint conditions δq 0 = δq N = 0. Taking into account the arbitrariness of the variations, we obtain the discrete Euler-Lagrange (DEL) equations for n = 0, . . . , N − 1. In total, (55) provide (3 d + 2 m) equations for the determination of the unknowns (q n+1 , p n+1 , v n , λ n , γ n+1 ) in every time step. Scheme (55) can be viewed as discrete counterpart of the continuous Euler-Lagrange equations (37). Note that relation (55c) can be interpreted as the discrete fiber derivative of the Legendre transformation, which links velocity and momentum variables. Moreover, the primary constraints (55d) are enforced at the end of each time step, while the secondary constraints are enforced in an intermediate sense (cf. (55e)).
It is worth mentioning that scheme (55) can be regarded as generalization to constrained mechanical systems of the symplectic Euler method (see Hairer et al. [33]). For the sake of completeness, we verify in Appendix A.6 that scheme (55) is symplectic and capable of conserving momentum maps of the form (45).

One-step GGL variational integrators with improved accuracy
The specific GGL variational integrator devised in the previous section is first-order accurate. We next aim at one-step GGL variational integrators with enhanced numerical accuracy. To this end, we borrow ideas from the discretization approach developed in Betsch and Becker [35] within the context of optimal control problems. The developments therein follow the philosophy first discretize then optimize, which is closely related to the present task of discretizing the GGL functional. Accordingly, we introduce the discrete GGL action sum Here, auxiliary coordinates Q n and momenta P n have been introduced which will be determined below. Furthermore, the functions and have been introduced, where L d , g d and G d will be specified in the sequel. In this connection, represents the discrete Lagrangian. Similarly, g d and G d stand for discrete versions of the constraint func-tion and the constraint Jacobian, respectively. The specific discretization (56) of the GGL functional can be motivated in the following way. In each time interval [t n t n+1 ], side condition (38) is discretized in the spirit of a one-stage Runge-Kutta method (see Hager [36]). In this regard, Q n is closely related to the stage value in a Runge-Kutta method. Moreover, p n+1 and P n play the role of Lagrange multipliers. Evaluating the discrete stationary conditions emanating from (56), we obtain for n = 0, ..., N − 1. Due to the arbitrariness of δp n+1 and δP n , (60a) and (60b) imply that Performing again an index shift in (60c) and taking into account the arbitrariness of the variations, the discrete equations of motion can be obtained as together with the discrete constraints for n = 0, ..., N − 1. In (62f), the discrete constraint function has been introduced. Note that (62d) serves as a discrete version of the fiber derivative (37c) and links velocity to momentum variables. In total, the above DEL equations constitute a set of (4 d + 2 m) equations for the determination of the unknowns (q n+1 , p n+1 , P n , v n+1 , λ n , γ n ).
Depending on the specific choice for L d , G d and g d , the above DEL equations generate a family of symplectic integrators (cf. Appendix A.7.1). In addition to that the above family of integrators is capable of conserving momentum maps of the form (45), provided that the following discrete symmetry conditions are satisfied by the specific choice for L d , G d and g d (cf. Appendix A.7.2): Of course, the discrete symmetry conditions are closely connected to the continuous symmetry conditions (41).

One-stage theta methods
To particularize specific integrators emanating from (62), we introduce the convex combination for θ ∈ [0, 1]. We now are in a position to choose It can be verified that this choice satisfies the discrete symmetry conditions (64a) and (64b), provided that the continuous conditions (41a) and (41b) hold. Based on the choice (66), functions (57) and (58) can be recast in the form and Now the partial derivatives appearing in (62) can be written as Moreover, (62d) yields Mv n+1 = p n+1 + P n , ( 7 0 ) while (62c) leads to Adding up (62b) and (62c) yields Inserting D 1 L(q n+θ , v n+1 ) from the last equation into (71), one obtains Substituting P n from the last equation into (70) eventually yields where has been introduced.
To summarize, the θ -version of the DEL (62) can be written in the form Note that we still have to specify the discrete constraint function g d (q n , q n+1 ) to obtain a specific integrator.
For that purpose, we shall consider the ensuing two options.

Option A: Intermediate constraint evaluation
The first option is based on the choice It can be shown that this choice satisfies the discrete symmetry condition (64c) under the assumption that the continuous condition (41c) holds. Choice (77) gives rise to the partial derivatives Inserting these relations into (76), the resulting DEL are given by The above DEL can be used for θ ∈ (0, 1).

Option B: Trapezoidal constraint evaluation
The second option relies on the introduction of a second parameter ϑ ∈ (0, 1] to control the evaluation of the constraint on configuration level. Specifically, we choose Since consistent initial conditions for q 0 = q(t = 0) ensure that g(q 0 ) = 0, the discrete constraint equation (80) enforces the configuration constraint at t n+1 for arbitrary ϑ ∈ (0, 1]. It can be shown that (80) satisfies the discrete symmetry condition (64c) under the assumption that the continuous condition (41c) holds. The partial derivatives of (80) yield Hence, (76) gives rise to the DEL for θ ∈ [0, 1] and ϑ ∈ (0, 1]. For example, the choice θ = 1 and ϑ = 0.5 yields the specific integrator

Energy-momentum scheme
In this section, we present a new energy-momentum (EM) conserving scheme which conserves both primary and secondary constraints. This scheme results from a direct discretization of the Euler-Lagrange equations (37) emanating from the GGL principle. Specifically, we consider where the derivatives in (37) have been replaced by corresponding discrete derivatives in the sense of Gonzalez [26]. In (84) Note that in the underlying symmetry conditions (41), the infinitesimal generators (17), (19) and (21) have been evaluated in the midpoint to arrive at (85). 2. Directionality conditions (cf. Gonzalez [26]) For systems with a constant mass matrix, the EM scheme (84) can be related to the method presented by Gonzalez [14] in a Hamiltonian framework. The present scheme, however, additionally preserves the secondary constraints and does not require to set up a Hamiltonian as it solely relies on the Lagrangian L.

Conservation of energy
To show that scheme (84) conserves the energy (4), take the scalar product of (84a) with (p n+1 − p n ), the scalar product of (84b) with −(q n+1 − q n ), the scalar product of (84c) with −h(v n+1 − v n ), and add the resulting equations to obtain Taking into account the directionality properties (86) of the (partitioned) discrete derivatives, we get The last two lines vanish due to (84d) and (84e). Accordingly, algorithmic conservation of the total energy holds in the sense that

Conservation of momentum maps
To show that scheme (84) is capable of conserving momentum maps of the form (45), take the scalar product of (84a) with ξ T p n+1/2 , the scalar product of (84b) with ξ q n+1/2 , the scalar product of (84c) with hξ v n+1/2 , and add the resulting equations to get provided that the orthogonality conditions (85) of the discrete derivatives are satisfied.

Numerical investigations
Throughout this section, the previously developed integration schemes are analyzed by means of numerical examples. In particular, we focus on the schemes summarized in Table 1. Newton's method has been applied to solve the algebraic system of equations corresponding to the respective scheme. If not otherwise specified, the known quantities at t n have been taken as initial guess for the computation of the unknowns at t n+1 . The tolerance in Newton's method has been set to tol = 10 −9 for all of the subsequent examples. The computations have been performed using the open source package metis, which is available at [37] together with detailed simulation data.

Mathematical pendulum in 3D
As a first example, we investigate the mathematical pendulum in three dimensions (Fig. 1), consisting of a massless rod of length l = 1 which forces the point mass with mass m = 1 onto its configuration manifold S 2 , the unit sphere. Gravitation acts in the negative e 3 -direction leading to the potential function V (q) = −mb · q, where b = −9.81e 3 is the gravitational acceleration. Moreover, the kinetic energy is given by T (v) = 1 2 mv · v leading to the mass matrix M = diag(m, m, m). Correspondingly, Lagrangian (3) is given by The configuration constraint assumes the form giving rise to the secondary constraint The present problem has rotational symmetry about the e 3 -axis. Accordingly, matrix A α in (16) is orthogonal (i.e., G = SO(3), the special orthogonal group in R 3 ) and can be written in the form whereê 3 is a skew-symmetric matrix with associated axial vector e 3 satisfyingê 3 a = e 3 × a for any a ∈ R 3 . To show that invariance property (11) holds, consider where orthogonality property (A α ) T = (A α ) −1 along with (A α ) T b = −9.81 exp(−αê 3 )e 3 = −9.81e 3 = b have been used. Similarly, it can be easily shown that invariance properties g(q α ) = g(q) and g v (q α , p α ) = g v (q, p) are satisfied. Moreover, equation (95) leads to Now, it can be easily verified that the symmetry conditions (41) are satisfied. The corresponding momentum map (45) takes the form Accordingly, the angular momentum J 3 of the pendulum is a first integral of the motion. The discrete derivatives featuring in the EM scheme (84) need to be applied to the functions L(q, v), g(q) and g v (q, p). Since these functions are at most quadratic, the discrete derivatives coincide with the respective midpoint derivatives. That is, It can be easily checked that orthogonality conditions (85) are satisfied.
Being initially released from position q 0 = [1, 0, 0] T with a velocity of v 0 = [0, 1, 0] T , the pendulum undergoes a motion which is exemplarily shown in Fig. 6, where the initial configuration and the configuration at t = T is depicted. The momenta have been initialized with p 0 = Mv 0 . The total simulation period is T = 10, and the time step size has been chosen as h = 0.05.
The 1-component of the positions and momenta are depicted in Figs. 2 and 3, respectively, to show that slight deviations between the different schemes occur. Those differences can also be viewed with respect to    The results for kinetic energy T , potential energy V and Hamiltonian H are exemplarily displayed for the VI-S scheme in Fig. 7. All the VIs (VI-S, VI-A and VI-B) exhibit fluctuations of the energy function E that are typical for symplectic methods (cf. Fig. 8). In contrast to that, as expected, the EM scheme conserves the energy of the pendulum up to numerical round-off. This can be seen from Figs. 8 and 9.
Furthermore, the evolution of the angular momentum components is exemplarily displayed for the VI-S method in Fig. 10. The other schemes yield only The primary constraint on position level is correctly captured by VI-S, VI-B and EM, as can be seen from Fig. 12. In contrast to that, VI-A violates the primary constraint by design (Fig. 13).  Similarly, the secondary constraint on velocity level is correctly captured by VI-S, VI-B and EM (Fig. 14), while VI-A violates the secondary constraint by design (Fig. 15). Although VI-S shows the correct behavior in this simple example, this is not the case anymore for more involved examples (see below). This example has been taken from Betsch and Leyendecker [38] (see also [16,17]). In particular, the rigid body is described as constrained mechanical system, see [39,40] for further details.
Accordingly, the configuration of the rigid body is characterized by the coordinate vector where ϕ ∈ R 3 is the position vector of the center of mass and d i ∈ R 3 (i = 1, 2, 3) constitute a body-fixed director frame (Fig. 16). For simplicity, the director axis are assumed to coincide with the principal axis of the rigid body. In the Lagrangian (3), the potential energy is given by V (q) = −mb · ϕ, where b = −9.81e 3 is the gravitational acceleration. The mass matrix takes the form The inertia quantities are given by the total mass m of the rigid body along with the principal moments of inertia I i (i = 1, 2, 3) about its center of mass, from which E i in (101) follow as E i = 1/2 I j + I k − I i for even permutations of the indices (i, j, k).
There are two types of holonomic constraints leading to the constraint function To enforce the orthonormality of the director frame, internal constraints of the form are imposed. To fix the tip of the gyroscopic top to the origin of the inertial frame (Fig. 16), the external constraint is used. Accordingly, to describe the mechanical system at hand, the n = 3 + 9 = 12 redundant coordinates are subjected to m = 6 + 3 = 9 constraints. Similar to the last example, the present problem has rotational symmetry about the e 3 -axis. Accordingly, matrix A α in (16) is given by leading to It is straightforward to verify that symmetry conditions (41) are satisfied. 2 The corresponding momentum map (45) takes the form so that the 3-component of the total angular momentum of the rigid top is a first integral of the motion.
Since the functions L(q, v), g(q), and g v (q, p) are at most quadratic, the discrete derivatives required in the EM scheme (84) coincide with the respective midpoint derivatives.
The data used in the numerical example are as follows. The total mass of the top amounts to m = 0.7069, the moments of inertia are I 1 = I 2 = I 3 = 5.3014 · 10 −4 and the center of mass is located along the symmetry axis at l = 0.075.
The initial nutation angle is α 0 = π/3. Correspondingly, the initial coordinates in (100) are specified by using the rotation matrix such that and The gyroscopic top is subject to an initial angular velocity vector of where the initial precession rate is chosen as ω p = 10.
The initial velocities corresponding to the coordinates The conjugate momenta have been set to p 0 = Mv 0 . In order to achieve the case of steady precession, the initial spin rate has to fulfill (cf. p. 221 in Goldstein [31]). In the present example, it amounts to ω s = 135.6. By ensuring the condition for steady precession (112), the center of mass has to rotate on a constant height of ϕ 3,ana = 0.0375 about the vertical axis. In addition to the present integration schemes (cf. Table 1), we apply the scheme proposed in Bobenko and Suris [41] (see also [42]), which has been developed specifically for the Lagrange top. This scheme is summarized in Appendix A.8 and labeled "BS99".
None of the schemes under investigation satisfies the condition of steady precession (i.e., ϕ 3 = ϕ 3,ana = 0.0375) as shown for h = 0.002 in Fig. 17 due to reasons of numerical accuracy.
To check the order of accuracy, the relative error in the vertical coordinate of the center of mass is calculated according to As shown in Fig. 18, VI-S, BS99 and VI-B are firstorder accurate, whereas VI-A and EM are second-order accurate. Figure 19 depicts the error for the momentum Again, as expected, only the EM scheme is capable of conserving the total energy of the system. However, all integrators at hand conserve the total angular momentum of the top about the vertical axis. In contrast to the pendulum example from Sect. 6.1, VI-S is not able to fulfill the velocity constraints exactly. Representative constraints are depicted in Fig. 20. Note, however, that the linear external constraints (i.e., k ∈ {7, 8, 9}) are still exactly fulfilled down to computer precision. As already pointed out in the previous example, VI-A does not fulfill both types of constraints in the endpoints, while VI-B and EM, by design, exactly capture the constraints both on configuration and on velocity level.

Four-particle system
The goal of this next example is to analyze the newly devised schemes with respect to robustness and the required computational effort. Moreover, we compare the EM method with the original scheme by Gonzalez [14]. The example of the four-particle system depicted in Fig. 21 has been taken from Gonzalez [14]. The configuration of the system is characterized by the configuration vector The two nonlinear springs give rise to the potential function where the spring stiffnesses are set to k 13 = 50 and k 24 = 500. The natural lengths of the two springs are l 13 = l 24 = 1, respectively. Obviously, the mass matrix is block diagonal and contains the masses m i (i = 1, . . . , 4). Specifically, m 1 = 1, m 2 = 3, m 3 = 2.3, m 4 = 1.7. There are two holonomic constraints with associated constraint functions where l 12 = l 34 = 1, respectively.
To summarize, the system has n = 12 redundant coordinates subject to m = 2 constraints. Moreover, the system has two symmetries. The first symmetry is related to the translational invariance of the system characterized by q α i = q i + αe for arbitrary e ∈ R 3 . Correspondingly, (10a) yields ξ T Q = [e T , e T , e T , e T ], so that momentum map (15) results in Since e is arbitrary, the total linear momentum L:= 4 i=1 p i is a first integral of the motion. The second symmetry is related to the rotational invariance of the system. With regard to (16), rotational invariance is characterized by for arbitrary e ∈ R 3 . The last equation yields It is straightforward to verify that symmetry conditions (41) are satisfied. The corresponding momentum map (45) takes the form Since e is arbitrary, the total angular momentum J:= 4 i=1 q i × p i is a first integral of the motion. Since the kinetic energy, g(q) and g v (q, p) are at most quadratic, the corresponding discrete derivatives required in the EM scheme (84) coincide with the midpoint derivatives. In contrast to that, the higher nonlinearity of potential function (116) demands a proper design of the corresponding discrete derivative (cf. Gonzalez [26]). To this end, we introduce the quadratic invariants Note that these invariants account for both symmetry properties. Now, potential function (116) can be recast in the form wherẽ The discrete derivative of V (q) can be viewed as discrete counterpart of the chain rule result in which the classical Greenspan's formula [43] is applied to obtain It is straightforward to verify that both orthogonality condition (85) and directionality condition (86) are satisfied.
In the numerical simulations, the following initial conditions have been applied: and v 0 The momenta have been initialized with p 0 i = m i v 0 i . The simulations have been conducted with h = 0.01 for a total time of T = 10.
The motion of the four-particle system is illustrated with some snapshots in Fig. 22.
For the sake of comparability, we have also included results of the original EM scheme developed by Gonzalez [14]. This scheme only takes into account the primary constraints and is labeled as "EM-G". Although EM-G does not exactly satisfy the secondary constraints (Fig. 23), which is in contrast to the present EM scheme (Fig. 24), both schemes yield comparable results concerning conservation of total energy (Fig. 25) conservation of angular momentum (Figs. 26 and 27) and accuracy.
Next we investigate the numerical effort caused by the schemes at hand. To this end, we compare the computation time along with the mean number of Newton iterationsn iter required per time step. In particular, the computation time t comp is normalized with respect to the time required by VI-S. The results in Table 2 refer to the use of two alternative initial guesses for q n+1 in Newton's method. In the first two columns of Table 2, q n has been used as initial guess (" simple initial guess"), while the results in the second column rely   on the initial guess q n +h v n (" advanced initial guess"). It can be observed that EM causes the largest numerical effort due to the computation of discrete gradients. Moreover, an advanced initial guess for the Newton method serves the purpose to decrease the necessary amount of iterations. Figure 28 shows that both EM and EM-G are secondorder accurate in the coordinates. Specifically, the relative error for the position of the fourth particle at t = 0.1 has been calculated via where q 4,ref represents the reference solution calculated with EM for h = 10 −5 . Similarly, both schemes yield second-order accuracy in the momenta (cf. Figure 29) and first-order accuracy in the Lagrange multipliers λ as can be observed from Fig. 30. Eventually, we focus on the numerical stability of the schemes under investigation. Since only EM is capable of conserving the total energy, it exhibits superior numerical stability properties (see Gonzalez and Simo [44] for a related stability analysis). This is confirmed by the results depicted in Figs. 31, 32 and 33. For h = 0.01, all of the schemes under investigation remain stable (cf. Fig. 31). In fact, for h ≤ 0.04 all schemes under investigation keep stable at least until T = 1000. The lower-order schemes VI-S and VI-B experience numerical energy blow-up already for h = 0.05 as shown in Fig. 32. For h = 0.25 (Fig. 33), VI-A exhibits an energy blow-up at about t = 55. Interestingly, EM keeps stable even for quite large time steps of h = 0.675.
It can be concluded that the superior numerical stability of EM comes at the expense of a larger numerical effort per time step.

Conclusion
The present approach to the numerical integration of the DAEs governing the motion of constrained mechanical systems relies on an index reduction in the spirit of the well-known GGL stabilization. We have presented the GGL variational principle which can be viewed as an extension of the Livens principle to constrained mechanical systems. The judicious inclusion of the secondary constraints into the GGL principle yields the desired index reduction in the spirit of the GGL stabilization. Accordingly, the Euler-Lagrange equations of the GGL principle assume the form of index 2 DAEs. These DAEs have been shown to preserve main structural properties of the underlying mechanical system independently of the multiplier associated with the secondary constraints. This feature is in sharp contrast to the original GGL stabilization and makes possible the design of structure-preserving discretization methods. In particular, the availability of the GGL functional has facilitated the design of novel variational integrators for the constrained mechanical systems at hand. These integrators have been shown to be symplectic and capable of conserving angular momentum maps. In addition to that, a new energy-momentum scheme has been developed by discretizing the Euler-Lagrange index 2 DAEs pertaining to the GGL principle.
While the GGL functional in the form (33) is valid for configuration-dependent (or non-constant) mass matrices, the development of variational integrators has been restricted to mechanical systems with constant mass matrix. The inclusion of non-constant mass matrices would be of interest for future work. Similar observations hold for the EM scheme (84) presented in Sect. 5. Apart from the inclusion of the velocity constraints, scheme (84) in general differs from the EM scheme due to Gonzalez [14] in the case of non-constant mass matrices. This holds true even for the unconstrained case dealt with in [26].
As the numerical investigations have shown, the present integrators are at most second-order accurate in the coordinates and the momenta. The development of higher-order variational integrators should be possible by following the ideas presented in Wenger et al. [45]. In this connection, the approach recently developed in Altmann and Herzog [46] should also be applicable. Furthermore, it would be of interest to investigate constrained symplectic partitioned Runge-Kutta methods (see Jay [12] and Marsden and West [23]) in the framework of the GGL principle.
Funding Open Access funding enabled and organized by Projekt DEAL. This work has been supported financially by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -Projektnummern 227928419 and 442997215.

Data availability
The generated and analyzed data sets as well as the corresponding source code are available at [37] Declarations Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

A.1 Properties of the wedge product
The following properties of the wedge product are of importance (see, for example, Leimkuhler and Reich [13]): Here, da, db and dc ∈ R d are arbitrary vector-valued differential one-forms, β 1 , β 2 ∈ R are scalar-valued quantities, A ∈ R d×d is an arbitrary matrix and B = B T ∈ R d×d is symmetric.
A.2 Conservation of the symplectic two-form along solutions of (2) The conservation of the symplectic two-form (23) can be verified by calculating The equations of motion (2) give rise to the differential one-forms where D 2 i j (•) denote the second-order partial derivatives. Accordingly, the temporal evolution of Ω becomes Since both D 2 11 L and D 2 22 L are symmetric matrices, the second and third term on the right-hand side of the above equation vanish due to property (129d) of the wedge product. Taking into account properties (129c) and (129a) of the wedge product, one arrives at the result due to the equality of mixed derivatives.

A.3 Original GGL stabilization and the symplectic twoform
To investigate how the GGL stabilization affects conservation of the symplectic two-form, we consider where (30a) and (30b) have been used. Due to property (129d), the first and third term vanish. The fourth term becomes zero as well since where (129c) and (129d) have been used along with (30c), which implies that Dg(q) dq = 0. Eventually, we arrive at Therefore, the conservation of the symplectic two-form Ω requires γ = 0, which holds for the continuous case.

A.4 Euler-Lagrange equations of the GGL princple
For the time-continuous case, the DAEs (37) yield γ = 0. This can be shown along the lines of the original GGL stabilization [3]. Differentiating (37d) with respect to time yields Substituting (37a) into the above equation leads to Assuming a system with Lagrangian (3), the first term on the left-hand side of the last equation vanishes due to (37e) and (37c). Thus, it remains which leads to γ = 0 since M is positive definite and G(q) has full rank. Moreover, the DAEs (37) have differentiation index 2. To show this we differentiate (37e) with respect to time to obtain Substituting from (37b) forṗ into the last equation and taking into account the above result γ = 0 yields where remaining terms have been collected in function g. Since GM −1 G is non-singular, the above equation can be solved for λ. One more time differentiation yields an explicit expression forλ and the index is thus equal to two.

A.5 Invalid form of the variational principle
The Euler-Lagrange equations associated with functional (40) are given bẏ Assuming a Lagrangian (3), relation (142c) can be used and inserted into (142e) such that Noting that D 2 g v (q, v) = G(q), we can find an explicit expression for the Lagrange multiplier as with M(q) = G(q)M −1 G(q) T . Correspondingly, the Euler-Lagrange equations (142) cannot be traced back to the standard formulation (27), since γ does in general not vanish.
A.6 Structure-preservation of scheme (55) For the purposes of this section, it is convenient to recast scheme (55) in the form p n+1 − p n = hD 1 L(q n , v n ) − hG(q n ) T λ n (145b) g(q n+1 ) = 0 , (145d) where g v (q, p n+1 ) = G(q)M −1 p n+1 has been used, which is in line with (29). Next, we demonstrate the structure-preserving properties of this scheme.
Substituting from (146a) and (146b) into the last equation, we obtain dq n ∧ dp n − dq n+1 ∧ dp n+1 = dp n+1 ∧ h dv n + h dp n+1 ∧ d D 2 g v (q, p n+1 ) T γ n+1 We next insert dp n+1 from (146c) into the first term on the right-hand side of (149). Moreover, the third term on the right-hand side of (149) vanishes due to property (129d) of the wedge product. The fourth one vanishes in analogy to relation (135). Consequently, we obtain dq n ∧ dp n − dq n+1 ∧ dp n+1 Next, it is possible to cancel the first and fifth term on the right-hand side of (150) due to properties (129a) and (129c) of the wedge product since (D 2 12 L) T = D 2 21 L. The symmetric matrix multiplication property of the wedge product (129d) can be used once more to cancel the second term. Moreover, it is possible to collect the third and last term by using (53), yielding dq n ∧ dp n − dq n+1 ∧ dp n+1 Executing the remaining differentials leads to the expression dq n ∧ dp n − dq n+1 ∧ dp n+1 = The last two terms on the right-hand side of (152) cancel each other due to (129c) since (D 2 12 g v k ) T = D 2 21 g v k . The first two terms on the right-hand side of (152) can be collected such that (147b) can be taken into account. Accordingly, (152) gives rise to dq n ∧ dp n = dq n+1 ∧ dp n+1 , which shows that the present scheme is indeed symplectic.
The differentials of (62a) through (62c) take the form dq n+1 − dq n = h dv n+1 + d (D 3 g v d ) T γ n , (158a) dp n+1 − dp n + dP n (158b) The last three equations can be inserted consecutively into the right-hand side of (157). Now, a tedious but straightforward calculation along the lines of Appendix A.6, taking into account the differential of (62d) given by 0 = D 2 31 L d dq n + D 2 32 L d dq n+1 + D 2 33 L d dv n+1 − h dp n+1 + dP n , leads to dq n ∧ dp n − dq n+1 ∧ dp n+1 = − dλ n ∧ D 1 g d dq n + D 2 g d dq n+1 − dγ n ∧ D 1 g v d dq n + D 2 g v d dq n+1 + D 3 g v d ( dp n+1 + dP n ) .
The right-hand side of the last equation is zero, which can be seen by taking the differentials of the constraints (62e) and (62f). Accordingly, dq n+1 ∧ dp n+1 = dq n ∧ dp n .