The GGL variational principle for constrained mechanical systems

We present an extension of the Livens variational principle (sometimes also referred to as Hamilton-Pontryagin principle) to mechanical systems subject to holonomic constraints. The newly proposed principle embodies an index reduction in the spirit of the often-applied GGL stabilization and thus may be termed “GGL principle”. The Euler-Lagrange equations of the GGL principle assume the form of differential-algebraic equations (DAEs) with differentiation index two. In contrast to the original GGL-DAEs, the present formulation ﬁts into the Hamiltonian framework of mechanics. Therefore, the GGL principle facilitates the design of structure-preserving integrators. In particular, it offers the possibility to construct variational integrators. This is illustrated with the development of a new ﬁrst-order scheme which is symplectic by design. The numerical properties of the newly devised scheme are investigated for representative examples of constrained mechanical systems.


Introduction
Dynamical systems may be formulated in various ways. The well-known Lagrangian and Hamiltonian formalisms both consider descriptive energetic scalars and deploy certain operations on them to generate the system's equations of motion. Another formulation, which unifies both above-mentioned formalisms by means of independent position, velocity and momentum quantities has been proposed by Livens [1]. This Livens principle has been recently taken up by Bou-Rabee [2], Yoshimura & Marsden [3], and Holm [4] under the name of Hamilton-Pontryagin principle due to its close relation to the Pontryagin principle from the field of optimal control [5]. Livens principle allows for an advantageous universal description due to its mixed character.
A large variety of dynamical systems are subject to constraints, which reduce the degrees of freedom of the system and impose some constraint function to be satisfied. When describing the system with redundant coordinates, the equations of motion emerge as a set of differential-algebraic equations (DAEs), which combine both differential equations and algebraic constraint equations. It is to mention that the numerical treatment of DAEs requires some additional effort compared to purely differential equations (cf. Kunkel & Mehrmann [6]). As the constraints have to hold at every point in time (consistency condition), so-called secondary constraints on velocity level are induced.
In a vast majority of dynamical problems, one cannot find an analytical solution. Thus, in recent times, the focus of scientific research has become to derive numerical integration methods, which are capable of solving the equations of motion approximately. In this connection, the class of structure-preserving integrators seeks to inherit the conservation principles of dynamical systems in a discrete sense (cf. monographs such as Hairer et al. [7] or Leimkuhler & Reich [8]). The first contributions can be traced back to symplectic methods (see, e.g., de Vogelaere [9]). In the field of mechanics, structure-preserving integration schemes can be mainly divided into two different groups: variational integrators and energy-momentum integrators (see the book [10] for a summary of important developments).
Variational integrators approximate the action integral and are typically able to conserve the symplectic structure as well as the system's momentum maps in a discrete sense (cf. Lew & Mata [11]). These are consequences of the variational procedure of derivation (cf. Marsden & West [12]). The main idea to find discrete counterparts of the variational principles goes back to Maeda [13]. Based on this concept, Marsden & West [12] provided a framework of discrete Lagrangian and Hamiltonian mechanics. Until now, variational integrators have been developed for various applications, e.g., for constrained dynamical systems (see, e.g., Leyendecker et al. [14]).
Integration schemes for constrained dynamics typically only account for the constraints on configuration level (primary constraints) but not for the secondary constraints and thus have a differentiation index of 3. This may lead to numerical instabilities (cf. Yoshimura [15]). By replacing the primary constraints on position level with the secondary velocitylevel constraints, DAEs with index 2 are obtained. Thus, the numerical problems can be avoided but violations of the primary constraints induce the well-known drift phenomenon (cf. Simeon [16]). However, this issue can be alleviated by extending the system of unknowns and coupling the secondary constraints into the equations. The most famous technique, the Gear-Gupta-Leimkuhler (GGL) stabilization, traces back to Gear et al. [17] in 1985 and is widely used until today (see, e.g., Brüls et al. [18] and Arnold et al. [19]). This classical GGL formulation relies on the direct modification of the equations of motion. Yet, this procedure leads to a destruction of the Hamiltonian structure such that most GGL stabilized integration schemes are not symplectic.
To the best of the authors' knowledge, numerical integration schemes for constrained dynamics have been formulated either in a Hamiltonian or in a Lagrangian way. Moreover, variational integration schemes have not been yet constructed such that primary and secondary constraints are considered at once. Thus, this work tries to fill both gaps by introducing a novel Livens-based variational framework for the integration of dynamical systems accounting for both primary and secondary constraints. In particular, the new framework makes it possible to justify the commonly used GGL formulation in a variational sense. Contrary to the original version, the newly proposed formalism provides index 2 DAEs with a Hamiltonian structure. Moreover, the novel variational principle opens up the possibility to develop new variational integration schemes with a GGL-type stabilization. A first step in this direction is undertaken in the present work. In particular, a first-order scheme is presented which results from the direct discretization of the newly proposed GGL variational principle.

Livens principle
Consider a dynamical system with d degrees of freedom and positions q ∈ R d . From Hamilton's principle of least action, one can proceed by allowing the velocities to be independent variables v ∈ R d . Thus, the kinematic relationq = v has to be enforced by means of a Lagrange multiplier p ∈ R d . The corresponding augmented functional reads where L(q, v) is the Lagrangian. The functional (1) [1]). More recently, Marsden and coworkers [2,3] coined the name Hamilton-Pontryagin principle for this functional due to its close relation to the classical Pontryagin principle from the field of optimal control [5]. Due to its mixed character with three independent fields (q, v, p), it resembles the Hu-Washizu principle from the area of elasticity theory (cf. Washizu [21]). 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 computing the variations with respect to every independent variable, one obtains the equations of motion in the forṁ where D α (•) denotes the partial derivative with respect to the α's argument. With regard to (2c) the multiplier p can be identified as the conjugate momentum, which thus directly emanates from the principle. Within the framework of Hamiltonian dynamics, momentum variables have to be defined a priori or emerge from the Legendre transformation as a fiber derivative of L(q,q). Note that after reinserting (2c) into (2b) and making use of (2a), Livens principle traces back to the Lagrangian equations of the second kind. Consider a mechanical system with constant mass matrix M ∈ R d×d , such that the Lagrangian takes the form where V (q) is a potential function. Now (2c) yields p = Mv, so that (2a) and (2b) can be rewritten asq Note that D(•) represents the gradient operator. These equations correspond to the Hamiltonian form of the equations of motion. Making use of the phase space vector z = [q T , p T ] T and the symplectic structure matrix where I d×d ∈ R d×d denotes the d × d identity matrix, the Hamiltonian equations of motion readż where the standard Hamiltonian function corresponding to (3) is given by

Symplectic structure of dynamics
Let us firstly introduce a bilinear and skew-symmetric function : R 2d × R 2d → R acting on two elements ξ , η ∈ R 2d , which assume an ordering of components as in the phase space vector such that, exemplarily, Bilinearity refers to the fact that is linear in both arguments. Skew-symmetry implies that (ξ , η) = − (η, ξ ). The canonical structure matrix introduced in (5) gives rise to the symplectic two-form A map : R 2d → R 2d is called symplectic if it leaves the symplectic two-form invariant in the sense that where the original two-form is equal to the two-form of the transports of ξ and η under the linearization of . Figure 1 displays the mapping of elements ξ and η by and the preservation of the symplectic two-form, representing the oriented area for d = 1. Making use of the definition (9), the last equation leads to a symplecticness condition that reads The wedge product of two differential one-forms da ∈ R d and db ∈ R d acting on any two vectors ξ , η ∈ R 2d is given by (cf. Leimkuhler & Reich [8]) Thus, the symplectic two form (9) can be rewritten in terms of the wedge product as where the differential one-forms dq i , dp i extract the ith coordinate or momentum component, respectively, such that Omitting the arguments, the symplectic two-form can be rewritten more briefly in vector notation as = dq ∧ dp. (15) Note that this representation is only a briefer notation of (13) that still accounts for the summation of the wedge product of scalar one-forms. One can show that Hamiltonian flow maps are symplectic. Thus, it is equivalent to say that the symplectic two-form (15) is conserved along solutions of the Hamiltonian equations of motions (6), viz.
For differential one-forms, given in vector notation as da, db and dc ∈ R d , any scalar valued quantities α, β ∈ R, any matrix A ∈ R d×d and any symmetric matrix B = B T ∈ R d×d the wedge product defined in (12) has the following properties (see, for example, Leimkuhler & Reich [8]):

GGL method for constrained mechanical systems
Assume that the coordinates q are redundant due to the presence of m independent scleronomic, holonomic constraints g k : R d → R (k = 1, . . . , m). The constraints can be comprised in a column vector g ∈ R m such that Since all constraint functions shall be independent, the constraint Jacobian G(q) = Dg(q) is of rank m. As (18) is true for any point in time, the time derivative has to vanish accordingly (consistency condition). Thus, the constraints on velocity level or secondary constraints d dt g(q) = G(q)q = 0 (19) are induced. It is well-known that the motion of the constrained mechanical systems under consideration is governed by differential-algebraic equations (DAEs) which have differentiation index ν = 3. These equations of motion can be derived with a variational approach, which augments Livens principle (1). Accordingly, introducingŜ(q, v, p, λ) = S(q, v, p) + T 0 λ · g(q) dt , stating the stationary condition δŜ(q, v, p, λ) = 0 and eliminating the velocities as above leads to an extension of the Hamiltonian equations (4a), (4b) for constrained systems, such that the index-3 DAEs are obtained aṡ The classical GGL stabilization, which traces back to Gear et al. [17], represents an index reduction technique by minimal extension (see Kunkel & Mehrmann [6]). The main idea of the GGL stabilization is to couple the secondary constraints into the dynamics by making use of additional variables γ ∈ R m , such that the system of equations of motion is extended and the differentiation index drops to ν = 2. Correspondingly, the numerical ill-conditioning of index-3 DAEs is alleviated without having the drawback of drift phenomena. The resulting index-2 DAEs can be written in the forṁ Ever since, the GGL stabilization has been widely used and is thus of great importance (see, e.g., Brüls et al. [18]). Numerical methods can be constructed directly by discretizing the DAEs (21a)-(21d). Note, however, that due to the GGL modification of the kinematic equation (21a), the system (21a)-(21d) loses its Hamiltonian structure. For the time-continuous case, some algebra leads to γ = 0. Consequently, the GGL-DAEs boil down to the standard formulation (20a)-(20c).

Governing equations
The newly proposed GGL principle relies on a generalization of Livens principle (1) by considering Lagrange multipliers λ, γ ∈ R m to enforce the primary constraints (18) and secondary constraints (19), respectively. Imposing stationary on a corresponding augmented action integral yields the stationary conditions As one can see, relation (24b) requires some more effort in order to achieve the final Euler-Langrange equation. Therefore, δq can be replaced through integration by parts such that The endpoint conditions on admissible variations δq(0) = δq(T ) = 0 make the latter two terms vanish. Furthermore, the variation of the gradient of the constraint functions can be executed as In order to avoid the third-order expression D 2 g(q), the stationary condition (24b) can be written in terms of the individual constraint functions g k (q) for k = 1, . . . , m. Thus, the variation of the constraint gradients is given by where the constraint Hessian D 2 g k (q) ∈ R d×d . Consequently, the arbitrariness of the variations δq, δv, δp, δλ and δγ can be taken into account such that the governing DAEs are deduced asq Thus equations in the fashion of the standard GGL stabilization (21a)-(21d) are obtained with an additional term in the momentum equation. Note that (28c) represents the fiber derivative of the newly proposed variational principle, such that the Lagrange multiplier p denotes the conjugate momenta. By introducing the secondary constraints to the functional, the resulting Euler-Lagrange equations are DAEs with index ν = 2, similarly to the GGL stabilized equations of motion, whereas the standard DAEs of constrained dynamics have index ν = 3. Thus, the newly established DAEs (28a)-(28e) can be regarded as an extension to the classical GGL stabilization. Similar to the classical GGL stabilization, one obtains γ = 0 for the time-continuous case. However, the third term on the right-hand side of (28b) is of crucial importance as it maintains the Hamiltonian structure of the equations of motion, in contrast to the classical GGL method.

Hamiltonian structure
The equations of motion induced by the GGL principle have Hamiltonian structure, viz.
with a corresponding augmented Hamiltonian The above equations can be related to the Euler-Lagrange equations (28a)-(28e) of the GGL functional after elimination of the velocities by employing the Legendre transformation (28c). Next, we show that the equations of motion of the novel framework conserve the Hamiltonian exactly. For that purpose, we compute the time derivative of the augmented Hamiltonian (30) such that where the Hamiltonian equations of motion (29a)-(29d) have been considered. As the augmented Hamiltonian (30) is conserved along solutions of the equations of motion, also the Hamiltonian H itself is conserved since both constraints on configuration level and momentum level are identically zero such thatḢ = 0. In contrast to the original GGL formulation (cf. Sect. 2.3), this conservation law holds regardless of the actual value of the Lagrange multipliers γ .

Symplectic structure
We show that the GGL functional inherits the symplectic structure of Hamiltonian systems and thus the symplectic two-form is conserved along solutions of (28a)-(28e) or (29a)-(29d), respectively. We begin by deriving the total differentials based on the equations of motion as where we have introduced the distinct functions g q (q) = g(q) for the holonomic constraint on configuration level and g v (q, p) = Dg(q)M −1 p for the corresponding constraint on momentum level. It is straightforward to compute the temporal evolution of by means of the product rule such that into which the above differential equations (32a) and (32b) can be inserted. One can consider the symmetry of the Hessian of H and make use of the properties of the wedge product. Note that D 2 11 H and D 2 22 H both are symmetric matrices such that D 2 22 H (q, p) dp ∧ dp = 0, due to property (17d). Moreover, the two terms with the off-diagonal entries of the Hessian of H cancel each other out, because where in the last equation property (17c) has been used. The term stemming from the primary constraints can be written as where again (17c) and (17d) have been used along with (32c). Therefore, all terms emerging from the right-hand side of (33) cancel out except for those containing the constraint on momentum level g v (q, p). We therefore obtain where the first term vanishes in view of the total differential of the secondary constraint (32d) and the second one cancels due to the symmetry of the Hessian of g v k . Note that it has been taken into account that terms including D 2 11 g v k and D 2 22 g v k , respectively, cancel due to their symmetry. This proofs the symplecticness of the equations of motion emanating from the GGL functional. Again this property does not depend on the Lagrange multiplier γ , which is an advantage over the original GGL method by Gear et al. [17], for which γ = 0 is required in order to conserve .

GGL variational integrator
We next illustrate how the GGL principle introduced in Sect. 3 can be employed to derive a variational integrator. Subsequently, structure-preserving properties of the newly devised variational integrator will be considered.

Governing equations
Let us construct a time-stepping scheme by means of a direct discretization of the GGL functional (23). Enforcing the constraint on configuration level in the endpoint and the constraint on momentum level in an intermediate state and discretizing the velocity by means of an explicit Euler method, we obtain the discrete action integral where the configuration variableq = q n + h v n has been introduced. Stationary conditions can be applied directly to the discrete functional, yielding along with Applying an index shift in the second part of (40b) from n + 1 to n and taking into account the arbitrariness of all variations, we obtain the discrete EL equations for n = 0, . . . , N − 1. In total we have obtained a set of (3 d + 2 m) equations for the unknowns (q n+1 , p n+1 , v n , λ n , γ n+1 ) in every time step. These are discrete counterparts of the continuous EL equations given in (28a)-(28e). It is advantageous that, due to the enhancement of the discrete action integral, the secondary constraints are now taken into account as well (cf. relation (41e)). Note that relation (41c) can be interpreted as the discrete fiber derivative of the Legendre transformation, which links velocity and momentum quantities. It is worth mentioning that scheme (41a)-(41e) can be regarded as generalization to constrained mechanical systems of the symplectic Euler method (see Hairer et al. [7] and Euler-B in Leimkuhler & Reich [8]).

Conservation properties
It is clear that the primary constraints are correctly captured in every time step by design (see relation (41d)). The secondary constraints are enforced in an intermediate sense (cf. (41e)). Moreover, we can show that the integrator governed by (41a)-(41e) is symplectic. In order to demonstrate this, we calculate the differentials of (41a) to (41c). This yields where g v (q, p n+1 ) = G(q)M −1 p n+1 has been introduced in analogy to the continuous case. Moreover, D 2 12 L(q n , v n ) = D 2 21 L(q n , v n ) T = 0 as been taken into account which is valid for Lagrangians of the form (3). The differential forms of the constraint equations (41d) and (41e) read dg(q n+1 ) = G(q n+1 ) dq n+1 = 0, dg v (q, p n+1 ) = D 1 g v (q, p n+1 ) dq + D 2 g v (q, p n+1 ) dp n+1 = 0.
Now, making use of the skew-symmetry of the wedge product, property (17a), one can deduce that dp n+1 ∧ ( dq n+1 − dq n ) + ( dp n+1 − dp n ) ∧ dq n = dq n ∧ dp n − dq n+1 ∧ dp n+1 . (44) Substituting from (42a) and (42b) into the last equation, we obtain dq n ∧ dp n − dq n+1 ∧ dp n+1 = (45) We next insert dp n+1 from (42c) into the first term on the right-hand side of (45). Moreover, the third term on the right-hand side of (45) vanishes due to property (17d) of the wedge product. The fourth one vanishes in analogy to relation (37). Consequently, we obtain dq n ∧ dp n − dq n+1 ∧ dp n+1 = D 2 22 The symmetric matrix multiplication property of the wedge product (17d) can be used once more to cancel the first term on the right-hand side of (46). Moreover, it is possible to collect the second and fourth term, yielding dq n ∧ dp n − dq n+1 ∧ dp n+1 = (47) Executing the remaining differentials leads to the expression dq n ∧ dp n − dq n+1 ∧ dp n+1 = (48) In analogy to the proof given in Sect. 3.3, terms including D 2 11 g v k and D 2 22 g v k , respectively, cancel due to their symmetry. It becomes obvious that the last two terms on the right-hand side of (48) cancel each other out due to (17c) since (D 2 12 g v k ) T = D 2 12 g v k . The first two terms on the right-hand side of (48) can be collected such that (43b) can be taken into account. Eventually, the whole expression on the right-hand side of (48) vanishes and we obtain dq n ∧ dp n = dq n+1 ∧ dp n+1 , which shows that the present scheme is indeed symplectic.

Extension to non-conservative systems
Considering dynamical systems being subject to (possibly configuration-and velocitydependent) non-conservative forces F nc (q, v) ∈ R d , the motion is governed by the Lagrange-d'Alembert principle (see, for example, [22]). It augments Hamilton's principle of least action by adding the virtual work of the non-conservative forces such that δS GGL (q, v, p, λ, γ The corresponding equations of motion are obtained aṡ Obviously, the Hamiltonian of such systems is not conserved. Following the procedure from Sect. 3.2 in an analogously manner, we obtain the energy balance Thus, the change of total energy is related to the power of the non-conservative forces. In the following we want to focus on velocity-dependent viscous forces of the form which can be derived from the Rayleigh dissipation function by differentiating with respect to the velocities such that F nc (q, v) = −D 2 G(q, v). In the last equation, C(q) ∈ R d×d denotes the positive semi-definite dissipation matrix. Correspondingly, the rate of change of the Hamiltonian is strictly nonpositive and matches the dissipated power P v = −2G. We next show, how the variational integrator developed above can be extended to account for non-conservative forces of the form (53). To this end, we consider the discrete version of (50) given by Here, the discrete action sum S d is again given by (39) and F nc d represents a discrete approximation of the non-conservative forces. Evaluating the viscous forces at t n such that F nc d = C(q n )v n , one arrives at an integration scheme, which extends the previously deduced method (41a)-(41e) and obeys the discrete Euler-Lagrange equations: for n = 0, . . . , N − 1 withq = q n + h v n . For more details on variational integrators for dissipative systems, we refer to the work by Kane et al. [23].

Numerical examples
The objective of this section is to analyze the behavior of the previously derived variational integrator by means of representative numerical examples. Equations (41a)-(41e) are solved for the respective systems be means of Newton's method, where the tolerance has been set to ε tol = 10 −9 . The computation has been performed using the metis code, which is available at [24].

Steady precession of a gyroscopic top
In this first example, we investigate the motion of a rigid body which is modeled by means of a director formulation proposed by Betsch & Steinmann [25]. The director formulation can be directly linked to natural coordinates [26], as shown in [27]. This formulation describes rigidity using the orthonormality condition of three directors {d i } positioned in the center of mass ϕ. Every point of the rigid body is uniquely defined by its material coordinates X i with respect to the director frame's origin such that we are able to express the spatial placement as a function of the material coordinates and time according to where the summation convention applies. Eventually, it is possible to describe the motion of the rigid body by n = 3 + 9 = 12 redundant coordinates accumulated in the coordinate vector For simplicity, we assume that the directors {d i } coincide with the principal axes of the body. This framework consequently allows for a representation of the system's Lagrangian in the standard fashion (3) with the constant and diagonal mass matrix where principal values of the Euler tensor can be computed with the principal moments of inertia as E i = 1/2 J j + J k − J i for even permutations of the indices (i, j, k). The primary constraints enforce the directors to stay orthonormal for all times due to the rigidity of the body, viz.
Specifically in this example, a gyroscopic top, as depicted in Fig. 2, has been investigated for a total simulation time of T = 2 s with a time step size of h = 0.002 s computed with the symplectic variational integrator from Sect. 4.
The total mass of the top amounts to m = 0.7069 kg and the moments of inertia read J 1 = J 2 = J 3 = 5.3014 · 10 −4 kg m 2 . This examplarily corresponds to a symmetric cone with mass density ρ = 2700 kg/m 3 , height a = 0.1m, top radius r = a/2, and a location of the center of mass along the symmetry axis l = 3/4 a. In this case, the moments of inertia can be computed via and the total mass is given by m = 1 3 ρπr 2 a. Gravitation acts in the negative e 3 -direction with b = [0 , 0 , −9.81m/s 2 ] T such that the potential energy of the system only depends on the position of the center of mass ϕ, viz.
It is crucially important that the top is subject to an additional constraint which fixes the tip of the gyroscopic top to the origin of {e i } by enforcing that the center of mass is located on the axis of symmetry with a distance of l to the origin. The initial nutation angle is α 0 = π/3 and the gyroscopic top is subject to an initial angular velocity vector where the initial precession rate is chosen as ω p = 10 s −1 and the initial spin rate for the case of steady precession can be computed via the relation (cf. p. 221 in Goldstein [28]), which amounts to ω s = 135.6 s −1 for the present case. Note that g = 9.81 m/s 2 denotes the magnitude of gravitational acceleration here. The transformation from the global e i coordinate system to the initially inclined system can easily be done by the use of a rotation matrix R 0 ∈ SO(3) with the well-known property R T 0 = R −1 0 . For the present case, this reads which prescribes a rotation about the e 1 -axis with the angle α 0 as can be seen in Fig. 2. Thus, the initial velocities of the center of mass and the directors can be computed by taking the cross product with the initial configuratioṅ where ω 0 can be transformed to the global coordinate system first with the aid of R 0 to simplify the computation. By ensuring the condition for a steady precession (66) the center of mass is rotating on a constant height ϕ 3 around the vertical axis since gravitational forces and restabilizing effects due to the rotation are in an equilibrium. The horizontal coordinate of the center of mass ϕ 3 can thus be regarded as an analytical reference to the solutions. The results given by the symplectic integration scheme oscillates around this analytical solution, as can be seen in Fig. 10. The evolution of the energetic quantities T , V , and H is shown in Fig. 3. The total energy of the system is not conserved identically along the solutions of this integration scheme as it can be seen in Fig. 5 which displays the increments in H from one time step to another. However, as it is typical for symplectic methods (cf. Fig. 16 in Lew & Mata [11]), H (t) oscillates around its true value and the energy error remains stable. Furthermore, as  the gyroscopic top is subject to external forces acting in the e 3 -direction, the symmetry of the system reduces to a conservation of the angular momentum about the e 3 -axis such that L 3 = constant. This is correctly captured by the symplectic method and can be seen in Fig. 4. Differences in L 3 from one point in time to another are close to computer precision (cf. Fig. 6).
By design, the constraints on configuration level g q k (t) corresponding to (41d) are identically fulfilled which can be observed in Fig. 7 (each line corresponds to one specific k ∈ {1, . . . , 9}). In contrast to that, in each time step the secondary constraints g v k (t) are merely enforced in an intermediate configuration, corresponding to (41e). This leads to the results depicted in Fig. 8. We can moreover analyze the h-convergence of the symplectic variational integrator. Therefore we have investigated the relative error in the vertical coordinate of the center of mass, which is supposed to remain constantly ϕ 3,ana = 0.0375 m, after a total simulation time of t = 0.001 s for various time step sizes. In Fig. 9 we display the for different time step sizes h. It becomes visible that the present method is first-order accurate.

Double four-bar linkage
In the second example, we consider a classical benchmark system for multibody dynamics, the double four-bar linkage (cf. González et al. [29] and Bayo & Avello [30]). As depicted in Fig. 11, the closed-chain structure at hand consists of five rigid bars with uniformly distributed mass m = 1 kg and length l = 1 m and seven joints interconnecting them. The lower The configuration of the system is described by the same director-based framework as in the first example. Further details about the planar director formulation of rigid body dynamics can be found in [31]. For each body, we make use of the position vector to the center of mass ϕ (i) ∈ R 2 and two directors d (i) α ∈ R 2 (α ∈ {1, 2}) specifying the orientation of the rigid body. Thus, we obtain six redundant coordinates along with the mass matrix for each rigid body, where E (i) α = 1 2 J (i) and J (i) is the moment of inertia about the center of mass. The system's potential energy is given by To enforce rigidity each body is subject to three internal holonomic constraints Moreover, there are seven external constraints due to the revolute joints. They are exemplarily given by for the joints P 1 and P 4 , respectively. Eventually, for the double four-bar linkage, we have a set of 5 · 6 = 30 redundant coordinates and 5 · 3 + 7 · 2 = 29 scalar-valued holonomic constraints. Consequently, the system has one degree of freedom, e.g., the inclination angle θ (cf. Fig. 11). According to González et al. [29], this benchmark problem can be used to test the ability of numerical formulations to overcome singular states, which occur for horizontal configurations of the linkage (e.g., for θ = 0). With the present variational integrator there were no numerical problems throughout the simulation. Correspondingly, one obtains smooth evolutions of the position and velocity components. This is exemplarily shown for the 1-components of the position x 4 and the velocity v 4 in point P 4 in Figs. 12 and 13, respectively.
It is straightforward to compute the energetic quantities over time. The results can be found in Fig. 14, where the energy transfer between potential and kinetic energy becomes visible. Furthermore, it can be observed that the typical oscillations in the total energy of symplectic schemes are present again. Increments in the Hamiltonian are of the order of 10 −3 and remain stable (see Fig. 16). The angular momentum in two dimensions for the given director formulation can be computed according to [31] as with Due to the presence of gravity, the system's Lagrangian does not have rotational symmetry and thus the angular momentum varies over time (cf. Figs. 15 and 17). Concerning the holonomic constraints several observations can be made. The constraints on position level are identically fulfilled by the integration scheme in each time step. This is depicted both for the internal constraints (cf. Fig. 18) and the external constraints (cf.   By design of the present method, the velocity level constraints are not satisfied at the endpoints, but rather in an intermediate sense (cf. equation (41e)). Correspondingly, the in- ternal constraints on velocity level exhibit slight deviations at the endpoints (see Fig. 19), which are more pronounced in regions close to the singular configurations. Due to the linearity of the external constraints (e.g., relation (74b)), the constraint Jacobian is constant, such that the corresponding velocity constraints are exactly fulfilled at the endpoints as well (cf. Fig. 21).
Eventually, we include physical damping in the system by applying the method developed in Sect. 4.3. To this end, we add two viscous damping elements to the double four-bar linkage. Those elements with constants η 1 = η 2 = 0.5 Ns/m connect the points P 1 and P 5 as well as P 2 and P 6 . Defining the relative velocities between the respective points as v rel,1 and v rel,2 , respectively, the Rayleigh dissipation function (54) is given by Since P 1 and P 2 are fixed to the ground, the relative velocities are given by v rel,1 =φ (4) + l 2ḋ v rel,2 =φ (5) + l 2ḋ 1 .
(78b) Thus, the dissipation matrix introduced in (53) for the system at hand reads in which for i ∈ {1, 2}. As a result, the total energy H of the systems decays over time as it is dissipated by the two viscous damping elements. This can be observed in the energy plot in Fig. 22. In Fig. 23 it becomes visible that the lost amount of total energy in each time step is due to the dissipated power P n v of the damping elements.

Conclusion and outlook
In this work, a new variational principle for the analysis of constrained dynamics has been proposed. The underlying functional takes account of both primary and secondary constraints. Due to its mixed character with independent position, velocity and momentum quantities, it generalizes Livens principle [1] and thus unites Lagrangian and Hamiltonian viewpoints. By coupling constraints on position and velocity level into the equations we have obtained a set of DAEs which can be regarded as an extension of the well-known GGL stabilization [17]. Contrary to the original formulation, however, the emanating equations of motion have Hamiltonian structure. The novel GGL functional gives rise to DAEs with differentiation index 2 and thus circumvents the numerical problems of the standard index-3 DAEs pertaining to mechanical systems subject to holonomic constraints. We could show, that the formulation is symplectic and has Hamiltonian structure. The conservation principles of constrained dynamics can be carried over to the novel augmented formulation. We have demonstrated that, in analogy to the classical GGL formulation [17], in the time-continuous case the additional Lagrange multipliers need to vanish. However, in contrast to the original GGL formulation, this property is not required to retain the conservation of the Hamiltonian and the symplectic structure.
Based on the newly proposed variational principle, we have successfully derived a new first-order variational integrator. This integrator satisfies the primary constraints and is capable to conserve the angular momentum of the system. The secondary constraints have been taken into account in an intermediate sense. We could show that the method is symplectic, which is a typical property of variational integrators (cf. Marsden & West [12]). Moreover, due to an appropriate extension, a scheme which accounts for viscous damping has been derived.
The numerical properties of the present method could be demonstrated in representative examples of multibody dynamics, making use of a rotationless director-formulation.
The novel variational framework represents a promising basis for the construction of structure-preserving integration schemes. The method which has been deduced throughout this work can thus be seen as a starting point for further developments. In particular, due to the close relationship of the GGL principle to optimal control, previously developed direct methods based on the philosophy "first discretize then optimize" (see, for example, Betsch & Becker [32]) can be used to obtain higher-order variational integrators for constrained mechanical systems. These integrators are symplectic by design. Other approaches for the design of higher-order methods (see, e.g., Wenger et al. [33] or Altmann & Herzog [34]) might be of interest for the discretization of the GGL principle as well. Furthermore, slight modifications can be applied to obtain energy-momentum consistent integrators, which are second-order accurate and represent another important class of structure-preserving timestepping schemes.
We eventually note that in the present work we have put the focus on constrained mechanical systems with nonsingular mass matrices. Since singular mass matrices present additional challenges (see, for example, Udwadia & Phohomsiri [35] and García de Jalón & Gutiérrez [36]), this case would also be of interest for future investigations.
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/.