Feedback Integrators for Mechanical Systems with Holonomic Constraints

The feedback integrators method is improved, via the celebrated Dirac formula, to integrate the equations of motion for mechanical systems with holonomic constraints so as to produce numerical trajectories that remain in the constraint set and preserve the values of quantities, such as energy, that are theoretically known to be conserved. A feedback integrator is concretely implemented in conjunction with the first-order Euler scheme on the spherical pendulum system and its excellent performance is demonstrated in comparison with the RATTLE method, the Lie–Trotter splitting method, and the Strang splitting method.


Introduction
The method of feedback integrators was proposed in [1,2] to numerically integrate the equations of motion for dynamical systems in order to preserve their domain manifolds and first integrals. The method is summarized as follows. Suppose that there is a given invariant set Λ of a continuous-time dynamical system on a manifold P, where Λ can be the intersection of level sets of the first integrals of the system. Embed P into some Euclidean space, extend the system to the ambient Euclidean space, and modify it outside Λ to turn Λ into a local attractor of the resulting dynamical system in the ambient space. Then, trajectories originating from points in Λ generated by integration of the modified dynamics remain in Λ theoretically and near Λ numerically, irrespective of the choice of numerical integration schemes. It is rigorously shown [1,3,4] that the discrete-time dynamical system derived from any one-step numerical integrator with uniform step size h for the modified continuous-time system has an attractor Λ h that contains Λ in its interior and converges to Λ as h → 0+. In this procedure, the set of equations of motion of the modified dynamical system is called a feedback integrator. Feedback integrators can be implemented by any usual integration scheme such as Euler, Runge-Kutta, or Matlab ode45 in one single global Cartesian coordinate system for the ambient Euclidean space and do not require projecting numerical trajectories to a certain set or solving algebraic equations during integration.
Here, we propose a way to apply feedback integrators to mechanical systems with holonomic constraints, which is not addressed in [1,2]. Consider a symplectic manifold P, a symplectic submanifold S of P, and a Hamiltonian function H on P, where it is often the case that P = T * R n for some n and S is the cotangent bundle of a set of holonomic constraints in R n . The Hamiltonian function H defines a Hamiltonian vector field X H on P and the restriction H|S of H to S defines a Hamiltonian vector field, denoted X H|S , on S. In general, X H does not coincide with X H|S on S; so, we employ the celebrated Dirac formula to extend the vector field X H|S from S to P, such that the dynamical system extends from S to P. The manifold S is an invariant set of the extended dynamical system on P. Thus, we can apply the usual method of feedback integrators to the extended system on P such that S is numerically well-preserved. This paper is organized as follows. The Dirac formula is first reviewed for the sake of completeness; then, the procedure for implementing feedback integrators for mechanical systems with holonomic constraints is presented. A simulation study is carried out on the spherical pendulum system, which is a mechanical system with a holonomic constraint, to demonstrate the excellent performance of feedback integrators in preserving the constraint set, the total energy, and the vertical component of the angular momentum vector, in comparison with the RATTLE method, the Lie-Trotter method, and the Strang method. Refer to [5] for more information on the three methods. We conclude with a small-scale simulation of the planar pendulum to show that the feedback method generally gives rise to integrators that are computationally more efficient as well.

Related Work
The numerical integration of mechanical systems with holonomic constraints has been an area of active research interest over the past decades. In this section, we briefly compare our approach with existing methods from the literature. Holonomic constraints can be viewed as index-1 differential-algebraic equations (DAEs) [6,7]. These systems can be integrated by direct discretization or by reformulating the DAE as an equivalent set of ODEs to which a standard numerical integration scheme may be applied. The resulting integration scheme typically involves the solution of a nonlinear equation representing the constraint. Another approach starts from the Hamiltonian or variational nature of mechanical systems to come up with discretizations that preserve the symplectic structure and the constraint manifold [5,[8][9][10]. Such discretizations, the RATTLE algorithm [11] and its generalizations [12][13][14] chiefly among them, typically exhibit superior long-term integration properties compared with standard, nonsymplectic integration algorithms, and have found wide application in the numerical integration of mechanical and control systems [15][16][17]. This approach has also been extended to the case of classical field theories with constraints [18,19], or to systems with dissipation [20]. What all these methods have in common is that they seek to enforce the constraint equations directly, which requires specialpurpose integration algorithms. By contrast, the feedback integrator method described in this paper modifies the equations of motion directly, to approximately conserve the constraint equations and other integrals of motion. As we pointed out before, the advantage of this approach is that it allows for standard, off-the-shelf numerical integrators to be used, such as the Euler or Runge-Kutta method.

Main Results
We first review the Dirac formula [21,22], explain how to construct feedback integrators for mechanical systems with holonomic constraints, and design a feedback integrator for the planar and spherical pendulum systems to demonstrate its excellent integration performance and versatility in comparison with the RATTLE method, the Lie-Trotter method, and the Strang method.

Review of the Dirac Formula
In this section, we recall the construction of Dirac for the decomposition of the Hamiltonian vector field along a symplectic submanifold S ⊂ P, where (P, ω) is a symplectic manifold and S is symplectic, such that for all z ∈ S we have T z S ⊕ T z S ω = T z P. Suppose that we are given a Hamiltonian function H on P. We will work semiglobally as follows. Suppose that S locally is expressed as the zero level set of a function f : P → R 2k that has 0 as a regular value. Then, we have dim S = 2n − 2k where dim P = 2n and we denote f = ( f 1 , . . . , f 2k ). Since S is symplectic, we can restrict the Hamiltonian H to S and, pulling back the symplectic form to S, we have for v z ∈ T z S.

Proposition 1.
Each of the X f i | S is a section of TS ω . Furthermore, these sections are linearly independent and, at any point, span this distribution.
From the independence of the f i , a dimension count shows that the linearly independent X f i (z) span T z S ω , completing the proof. Now, define for each z ∈ P a (2k × 2k) matrix by We then have the following.
Proof. Fix z ∈ S. We know that T z S is symplectic and, therefore, T z S ω is too. By the previous proposition, we know that the X f i (z)'s form a basis of the symplectic space T z S ω . Therefore, the entry of the matrix C ij (z) is precisely the symplectic form evaluated on the vectors X f i (z), X f j (z), which is, thus, a non-degenerate matrix since the restriction of ω to T z S ω is non-degenerate.
We then have the following version of the Dirac formula for the Hamiltonian vector field.

Theorem 1.
With the definition of C ij (z) given as above, denoting its inverse by C ij (z) for z ∈ S, the following formula holds: Proof. Since S is symplectic, we have that for all z ∈ S, T z P = T z S ⊕ T z S ω . We will show that the projection π z : T z P → T z S relative to this decomposition is given by the right-hand side of (1). This is equivalent to showing that I − π z : T z P → T z S ω is given by To establish this, first observe that the right-hand side lies in T z S ω by Proposition 1.
which shows that the operator is the identity on the subspace T z S ω . Next, suppose that X H (z) ∈ T z S. Then, we have d f i (z) · X H (z) = 0 for all i and, therefore, the right-hand side of (2) vanishes, as required. This proves that the projection I − π z is given by the Formula (2) and, thus, π z is given by the right-hand side of (1).
Our main use of this theorem is to extend a Hamiltonian system on S to all of P. That is, given a Hamiltonian function H defined on P, we compute the right-hand side of Equation (1) for arbitrary z ∈ P to obtain a vector field defined on a neighborhood of the submanifold S (or level surface in this case) that will automatically be equal to X H|S on S itself. Thus, the Dirac formula gives us a way to extend Hamiltonian vector fields on S to the ambient space.

Feedback Integrators for Mechanical Systems with Holonomic Constraints
Given a Hamiltonian function H on a 2n-dimensional symplectic manifold P and its Hamiltonian vector field X H on P, consider a Hamiltonian system Σ S on a (2n − 2k)dimensional symplectic submanifold S of P whose Hamiltonian function is the restriction H|S of H to S. So, the equations of motion of Σ S can be written aṡ Assume that S is expressed as a level of a function f = ( f 1 , . . . , f 2k ) : P → R 2k . Thanks to the Dirac Formula (1), the dynamical system (3) extends to P aṡ It is easy to verify that S is an invariant manifold of (4). We now wish to integrate the dynamics (4) from a point x 0 in S.
Suppose that we can embed the manifold P in some Euclidean space R m with m ≥ 2n and extend the vector field (4) to R m , not necessarily in the Dirac way. Denote the extended vector field by X and write the corresponding dynamical system aṡ Suppose that both functions f and H also extend to the ambient Euclidean space R m in such a way that D f · X = 0 and ∇H · X = 0 in R m , and that the manifold S is still a level set of the extension of f in R m . As a result, the functions f and H are first integrals of (5) and the manifold S is an invariant manifold of (5).
Suppose that there are first integrals I = (I 1 , . . . , I ) of (5) other than H and f , where the function I may include part of a function on R m whose level set defines P in R m . Define a function F : R m → R 2k+1+ by F = ( f , H, I), and a function V on R m by The minimum value of V is 0 and the function V satisfies Modify the dynamical system (5) as follows: where L is a map from R m into the set of m × m positive definite symmetric matrices and ∇V is computed as ). Since 0 is the minimum value of V, the gradient ∇V vanishes on V −1 (0), which implies that the two dynamical systems (5) and (6) coincide on V −1 (0). It follows that the set V −1 (0) is an invariant manifold of (6). Along any trajectory of (6), V(t) is an nonincreasing function of t since . Under some conditions on V, the set V −1 (0) becomes a unique attractor of (6) in a neighborhood of itself in R m ; refer to [1] for those conditions. Due to the invariance of V −1 (0) and the coincidence of (5) and (6) on V −1 (0), integrating (5) from x 0 and integrating (6) from x 0 produce the same trajectory. Numerically, however, integrating (6) has the following advantage: if the trajectory deviates from V −1 (0) at some numerical integration step, then it will get pushed back toward the attractor V −1 (0), thus remaining on the manifold S and preserving all the first integrals; refer to [1,23] for a rigorous explanation. Although [1,23] provides some sufficient conditions for V −1 (0) to be a local attractor of (6), in practice, it is not necessary to verify them. The procedure outlined in this section is good enough for applications, which will be illustrated with the planar and spherical pendulum in the next sections.

The Spherical Pendulum
We build a feedback integrator for the spherical pendulum [24], which is a Hamiltonian system with a holonomic constraint, and compare its performance with that of such geometric numerical integration methods as the RATTLE method, the Lie-Trotter splitting method, and the Strang splitting method. The phase space of the spherical pendulum is T * S 2 , which is a symplectic submanifold of T * R 3 = R 3 × R 3 . This submanifold is globally defined as the ( 2 , 0)-level set of the function f = ( f 1 , f 2 ) : T * R 3 → R 2 with > 0, where f 1 (q, p) = ||q|| 2 and f 2 (q, p) = q · p for (q, p) ∈ T * R 3 . We fix the Hamiltonian H(q, p) = ||p|| 2 /2m + mgq 3 . Restricted to T * S 2 , this gives the Hamiltonian of the spherical pendulum under gravity, with its S 1 symmetry. In order to write down the extended Hamiltonian vector field, we compute the following: where e 3 = (0, 0, 1). Hence, by the Dirac formula, the spherical pendulum system extends to T * R 3 asẋ = X(x) where The extended system has four first integrals on T * R 3 : the constraint functions f 1 and f 2 , the Hamiltonian H, and the vertical component J(q, p) = q 1 p 2 − q 2 p 1 of angular momentum. Choose a point (q 0 , p 0 ) ∈ T * S 2 , and define a function V : where k i > 0, i = 1, . . . , 4 are constants; . It is easy to verify that 0 is the minimum value of V and where ∇ f 1 = (2q, 0), ∇ f 2 = (p, q), ∇H = (mge 3 , p/m), and ∇J = (p × e 3 , e 3 × q). Then, the feedback integrator corresponding to the function V for the spherical pendulum is given byẋ We now compare the feedback integrator with the RATTLE method, the Lie-Trotter method, and the Strang method on the spherical pendulum system. The RATTLE algorithm is given on p. 246 in [5] and is known to be symplectic and convergent of order two [5]. The Lie-Trotter method and the Strang method are so-called splitting methods and are explained on pp. 253-254 in [5], where the two splitting methods yield first-and secondorder numerical integrators, respectively [5]. For the splitting methods, we split the Hamiltonian function H into the kinetic function H [1] (q, p) = ||p|| 2 /2m and the potential function H [2] (q, p) = mgq 3 , and we note that the dynamics associated to H [1] and H [2] can be integrated analytically. For numerical simulation, choose the parameter values m = g = = 1 for convenience, and the initial points q(0) = (0, 1, 0) and p(0) = (1, 0, −1) on T * S 2 . The corresponding initial values of the first integrals f 1 , f 2 , H, and J are f 10 = 1, f 20 = 0, H 0 = 2, and J 0 = −1, respectively. We fix the time step size h = 10 −3 and the time interval [0, 100] for integration by all four methods. For the feedback integrator, the usual Euler scheme is used to integrate (7) and (8) with the following gain values: k 1 = 50, k 2 = 50, k 3 = 50, and k 4 = 50.
The computational results are plotted in Figures 1-5. In Figure 1, the trajectories q(t) = (q 1 (t), q 2 (t), q 3 (t)) of the configuration variables generated by the four methods are plotted. The feedback integrator with the Euler, RATTLE, Lie-Trotter, and Strang splitting methods all generate similar trajectories. Figures 2 and 3 show the plots of the deviations |∆ f 1 (t)| and |∆ f 2 (t)| from the constraint manifold T * S 2 . The result by the RATTLE method is the best, and the trajectory produced by the feedback integrator with Euler remains close to T * S 2 with the step size h = 10 −3 taken into account. Likewise, the trajectories by the Lie-Trotter method and the Strang method stay close to T * S 2 , because the flows of the split Hamiltonians H [1] and H [2] each preserve T * S 2 , as does their composition.
The feedback integrator with Euler and the RATTLE method perform well in preservation of the values of the Hamiltonian H as do the other two methods, as shown in Figure 4. The vertical component J of angular momentum is well-preserved by all four methods as shown in Figure 5, where it is noticeable that the Lie-Trotter method and the Strang method perform very well in preservation of J.
The computational results imply that the feedback integrator with the first-order Euler scheme performs well on the spherical pendulum system in comparison with the well-known RATTLE method that is of second order, and to the Lie-Trotter method and the Strang method. An advantage of the feedback integrator over the other three methods is that it does not require any projection or solving of algebraic equations to stay on the holonomic constraint manifold. Further, it does not require any special integration algorithms, and simply employs well-known integration algorithms available such as Euler, Runge-Kutta, or Matlab ode45. Moreover, unlike the Lie-Trotter and Strang methods, which require a particular splitting of the dynamics into two parts that are separately integrable, the feedback integrator can be made to work for any constrained dynamical system by modifying the vector field as in (6).

The Simple Pendulum
One potential drawback of the feedback integration method is that the feedback vector field is more complex due to the presence of the stabilizing forces. For example, the right-hand side of (7) and (8) is a sum of five gradients (one gradient of the Hamiltonian and four gradients of the feedback potential) and this cost compounds for higher-order numerical methods, which typically require several force evaluations per integration step. This added computational cost must be taken into account when comparing the error profile of feedback integrators with that of standard methods, such as RATTLE, which only evaluate the gradient of the Hamiltonian (but possibly multiple times per integration step).
In this section, we show that the increase in complexity is compensated by the approximate conservation properties of the integrator, and in particular, we show that feedback integrators are at least as effective as standard integrators when the computational budget is taken into account. We compare the performance of feedback integrators of different orders with that of the RATTLE method and show that the increase in computational cost is balanced by the fact that feedback integrators require fewer force evaluations overall to achieve a given accuracy.
To assess the global error of the various integrators, we take recourse to a simpler mechanical system, for which exact solutions are known and can be approximated with arbitrary precision: the simple gravitational pendulum. The simple pendulum consists of a mass m that is free to swing at the end of a rigid rod of length under the influence of gravity. The motion of the pendulum takes place entirely in a fixed plane, denoting the angle between the horizontal and the position of the pendulum by θ, is determined bÿ where g is the gravitational acceleration. The dynamics of the pendulum as a constrained system can be derived directly via a calculation as in the previous section, or by observing that the spherical pendulum naturally moves in a fixed plane if the initial position, momentum, and direction of gravity are all coplanar. Either way, the extended Hamiltonian vector field is readily seen to be where q = (x, y) and p = (p x , p y ) are the coordinates and momenta of the pendulum, respectively, and e y = (0, 1). The first integrals on T * R 2 are the Hamiltonian H(q, p) = p 2 /2m + mgy and the constraint functions Similar to the spherical pendulum, we consider these three conserved quantities and for given initial values (q 0 , p 0 ) ∈ T * S 2 , we define the function V : T * R 2 → R given by where k 1 , k 2 , k 3 are positive constants; ∆ f i = f i (q, p) − f i (q 0 , p 0 ); i = 1, 2; and ∆H = H(q, p) − H(q 0 , p 0 ). The feedback integrator for the simple pendulum then becomeṡ While the pendulum can, in principle, be integrated exactly, obtaining the solution as a function of time is not straightforward and requires inverting the elliptic integral of the first kind. To avoid this difficulty, we integrate the Equation (9) for θ using a high-order Runge-Kutta method with tolerance set to 10 −13 . For all numerical simulations, we set m = g = = 1, and we use q(0) = (1, 0) and p(0) = (0, 0) as the initial conditions. For the feedback method, we use three off-the-shelf numerical schemes to integrate (10) and (11): (a) the forward Euler method, (b) the explicit 4th-order Runge-Kutta method (RK4), and (c) the Dormand-Prince method of 8th order (DOP853). Note that the Dormand-Prince method is a variable step-size integrator while all the others use a fixed step size. Figure 6 (left) shows the trajectory error after one period of the pendulum between the standard RATTLE integrator on the one hand, and the three feedback integrators on the other hand, as a function of the step size. The global error decreases as the step size decreases, in line with the order of the method. This is to be expected, since the feedback approach merely modifies the vector field to be integrated, but does not otherwise alter the underlying numerical method.
The vector field integrated by the feedback methods is more complex; thus, one can ask what the impact is on the total execution time. To investigate this question, we give each method a fixed computational budget and modify the integrator code so that each evaluation of the vector field reduces the computational budget by a certain amount. For evaluations of the feedback vector field (10) and (11) the cost per evaluation is set to 4, since the vector field consists of 4 gradients, whereas for the evaluation of the gradient of the Hamiltonian (needed e.g., by RATTLE), the cost per evaluation is 1. We then compare the performance of the integrators over one period of the pendulum and adjust the step size (for RATTLE and the feedback Euler and RK4 methods) or the tolerance (for the feedback-DOP853 method, which uses a variable step size) so that the computational budget is exhausted over one period.
The result is shown in Figure 6 (right), which shows the global error as a function of the computational budget. Note that the error for the feedback-DOP853 method stabilizes somewhat below 10 −12 . This is roughly the point where we encounter the limits in the accuracy of the exact trajectories (which were obtained by numerical integration of (9), where the tolerance was set to 10 −13 ).
We see that feedback integrators are able to integrate the dynamics of the underlying system accurately (i.e., with low error) and efficiently (using comparable or lower numbers of force evaluations), compared with specific holonomic integrators. In terms of computational efficiency, higher-order integrators such as the feedback-DOP853 integrator clearly achieve better results than others (lower-order feedback methods and RATTLE). This again demonstrates one of the key benefits of the feedback integrator method, showing that it is possible to use any standard numerical integration scheme to obtain approximate constraint preservation, without loss of accuracy or computational efficacy. size (DOP853 is not included as it is a variable step size method). As the step size decreases, the global error decreases at a rate proportional to the error of the method. Right: The global error, but now as a function of the computational budget (number of force evaluations). Larger computational budgets correspond to smaller step sizes and, hence, lower errors, but take into account the fact that the feedback methods involve more force evaluations. Despite the overhead, feedback integrators are able to do at least as well as, or better than, RATTLE.

Conclusions
We have presented a general framework to extend the feedback integrators [1] to systems with holonomic constraints. Beginning with a symplectic submanifold S in the symplectic manifold P, where S-the holonomic constraint-is the level set of f 1 . . . f 2k , on which we have a Hamiltonian H|S. We use the Dirac formula to extend the vector field X H|S to a vector field X on P. We then apply the feedback integrator to X using an embedding of P in Euclidean space.
More specifically, in the case that the symplectic manifold S is of the form T * Q, where Q embeds in R n and T * Q embeds in T * R n as the level set of functions f 1 , . . . f 2k , on which we have an extended Hamiltonian H, we can compute the extension of the vector field X H|T * Q directly from the Dirac formula (1). With this vector field, now defined on a Euclidean space, T * R n , we construct the function V : T * R n → R ≥0 , whose 0 level set contains the dynamic invariants for a given initial point in S, where the set of dynamic invariants includes the functions f 1 , . . . f 2k . The feedback-integrator-modified vector field, constructed from the extended vector field X by adding the negative gradient of V, is fed into any integrator, for example, first-order Euler. In addition to the dynamic invariants, the integrator automatically respects, from the construction of the modified vector field, the holonomic constraints. The resulting integrator has superior performance to even symplectic integrators, which depend, typically, on implicit solvers. As future work, we plan to examine the problem of optimal control on manifolds and its relationship with Hamiltonian mechanics from the viewpoint of feedback integrators about which some preliminary works have been carried out in [25][26][27]. We also plan to examine the effect of the magnitude of feedback integrator term on the precision of numerical integration and to extend feedback integrators to field theory.