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.
2.1. 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
, where
is a symplectic manifold and
S is symplectic, such that for all
we have
. 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
that has 0 as a regular value. Then, we have
where
and we denote
. Since
S is symplectic, we can restrict the Hamiltonian
H to
S and, pulling back the symplectic form to
S, we have
for
.
Proposition 1. Each of the is a section of . Furthermore, these sections are linearly independent and, at any point, span this distribution.
Proof. We have that if and only if for all . Now, let . We then have, for each i, ; therefore, is a section of . From the independence of the , a dimension count shows that the linearly independent span , completing the proof. □
Now, define for each
a
matrix by
We then have the following.
Proposition 2. Fix . The matrix is invertible.
Proof. Fix . We know that is symplectic and, therefore, is too. By the previous proposition, we know that the ’s form a basis of the symplectic space . Therefore, the entry of the matrix is precisely the symplectic form evaluated on the vectors , which is, thus, a non-degenerate matrix since the restriction of to is non-degenerate. □
We then have the following version of the Dirac formula for the Hamiltonian vector field.
Theorem 1. With the definition of given as above, denoting its inverse by for , the following formula holds: Proof. Since
S is symplectic, we have that for all
,
. We will show that the projection
relative to this decomposition is given by the right-hand side of (
1). This is equivalent to showing that
is given by
To establish this, first observe that the right-hand side lies in
by Proposition 1. Next, if
where
, then the right-hand side of (
2) is
which shows that the operator is the identity on the subspace
. Next, suppose that
. Then, we have
for all
i and, therefore, the right-hand side of (
2) vanishes, as required. This proves that the projection
is given by the Formula (
2) and, thus,
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
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
on
S itself. Thus, the Dirac formula gives us a way to extend Hamiltonian vector fields on
S to the ambient space.
2.2. Feedback Integrators for Mechanical Systems with Holonomic Constraints
Given a Hamiltonian function
H on a
-dimensional symplectic manifold
P and its Hamiltonian vector field
on
P, consider a Hamiltonian system
on a
-dimensional symplectic submanifold
S of
P whose Hamiltonian function is the restriction
of
H to
S. So, the equations of motion of
can be written as
Assume that
S is expressed as a level of a function
. Thanks to the Dirac Formula (
1), the dynamical system (
3) extends to
P as
It is easy to verify that
S is an invariant manifold of (
4). We now wish to integrate the dynamics (
4) from a point
in
S.
Suppose that we can embed the manifold
P in some Euclidean space
with
and extend the vector field (
4) to
, not necessarily in the Dirac way. Denote the extended vector field by
X and write the corresponding dynamical system as
Suppose that both functions
f and
H also extend to the ambient Euclidean space
in such a way that
and
in
, and that the manifold
S is still a level set of the extension of
f in
. 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
of (
5) other than
H and
f, where the function
I may include part of a function on
whose level set defines
P in
. Define a function
by
, and a function
V on
by
where
is a
constant positive definite symmetric matrix. The minimum value of
V is 0 and the function
V satisfies
The set
is an invariant manifold of (
5) since
.
Modify the dynamical system (
5) as follows:
where
L is a map from
into the set of
positive definite symmetric matrices and
is computed as
. Since 0 is the minimum value of
V, the gradient
vanishes on
, which implies that the two dynamical systems (
5) and (
6) coincide on
. It follows that the set
is an invariant manifold of (
6). Along any trajectory of (
6),
is an nonincreasing function of
t since
; so, the trajectory converges to
as
or it remains close to
if it starts near
. Under some conditions on
V, the set
becomes a unique attractor of (
6) in a neighborhood of itself in
; refer to [
1] for those conditions. Due to the invariance of
and the coincidence of (
5) and (
6) on
, integrating (
5) from
and integrating (
6) from
produce the same trajectory. Numerically, however, integrating (
6) has the following advantage: if the trajectory deviates from
at some numerical integration step, then it will get pushed back toward the attractor
, 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
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.
2.3. 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
, which is a symplectic submanifold of
. This submanifold is globally defined as the
-level set of the function
with
, where
and
for
. We fix the Hamiltonian
. Restricted to
, this gives the Hamiltonian of the spherical pendulum under gravity, with its
symmetry. In order to write down the extended Hamiltonian vector field, we compute the following:
and
where
. Hence, by the Dirac formula, the spherical pendulum system extends to
as
where
for
. The extended system has four first integrals on
: the constraint functions
and
, the Hamiltonian
H, and the vertical component
of angular momentum.
Choose a point
, and define a function
by
where
,
are constants;
,
;
; and
. It is easy to verify that 0 is the minimum value of
V and
where
and
. The gradient of
V is computed as
where
,
,
, and
. Then, the feedback integrator corresponding to the function
V for the spherical pendulum is given by
or
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 second-order numerical integrators, respectively [
5]. For the splitting methods, we split the Hamiltonian function
H into the kinetic function
and the potential function
, and we note that the dynamics associated to
and
can be integrated analytically. For numerical simulation, choose the parameter values
for convenience, and the initial points
and
on
. The corresponding initial values of the first integrals
,
,
H, and
J are
,
,
, and
, respectively. We fix the time step size
and the time interval
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:
,
,
, and
.
The computational results are plotted in
Figure 1,
Figure 2,
Figure 3,
Figure 4 and
Figure 5. In
Figure 1, the trajectories
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.
Figure 2 and
Figure 3 show the plots of the deviations
and
from the constraint manifold
. The result by the RATTLE method is the best, and the trajectory produced by the feedback integrator with Euler remains close to
with the step size
taken into account. Likewise, the trajectories by the Lie–Trotter method and the Strang method stay close to
, because the flows of the split Hamiltonians
and
each preserve
, 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).
2.4. 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 by
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
and
are the coordinates and momenta of the pendulum, respectively, and
. The first integrals on
are the Hamiltonian
and the constraint functions
Similar to the spherical pendulum, we consider these three conserved quantities and for given initial values
, we define the function
given by
where
are positive constants;
;
; and
. The feedback integrator for the simple pendulum then becomes
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
. For all numerical simulations, we set
, and we use
and
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
. 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
).
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.