Lagrangian Mechanics and Variational Integrators on Two-Spheres
aa r X i v : . [ m a t h . NA ] J un LAGRANGIAN MECHANICS AND VARIATIONAL INTEGRATORS ONTWO-SPHERES ∗ TAEYOUNG LEE † , MELVIN LEOK ‡ , AND
N. HARRIS MCCLAMROCH § Abstract.
Euler-Lagrange equations and variational integrators are developed for Lagrangian mechanical systems evolvingon a product of two-spheres. The geometric structure of a product of two-spheres is carefully considered in order to obtain globalequations of motion. Both continuous equations of motion and variational integrators completely avoid the singularities andcomplexities introduced by local parameterizations or explicit constraints. We derive global expressions for the Euler-Lagrangeequations on two-spheres which are more compact than existing equations written in terms of angles. Since the variationalintegrators are derived from Hamilton’s principle, they preserve the geometric features of the dynamics such as symplecticity,momentum maps, or total energy, as well as the structure of the configuration manifold. Computational properties of thevariational integrators are illustrated for several mechanical systems.
Key words.
Lagrangian mechanics, geometric integrator, variational integrator, two sphere, homogeneous manifold
AMS subject classifications.
1. Introduction.
The two-sphere S is the set of all points in the Euclidean space R which are aunit distance from the origin. It is a two dimensional manifold that is locally diffeomorphic to R . Manyclassical and interesting mechanical systems, such as a spherical pendulum, a double spherical pendulum,and magnetic models, evolve on the two-sphere or on a product of two-spheres. In this paper, we deriveEuler-Lagrange equations on configuration spaces of the form ( S ) n , for a positive integer n . We also developgeometric numerical integrators referred to as discrete Euler-Lagrange equations or variational integratorson ( S ) n .In most of the literature that treats dynamic systems on ( S ) n , either 2 n angles or n explicit equalityconstraints enforcing unit length are used to describe the configuration of the system [2, 15]. These descrip-tions involve complicated trigonometric expressions and introduce additional complexity in analysis andcomputations. In this paper, we focus on developing continuous equations of motion and discrete equationsof motion directly on ( S ) n , without need of local parameterizations, constraints, or reprojections. Thisprovides a remarkably compact form of the equations of motion.Geometric numerical integrators are numerical integration algorithms that preserve the geometric struc-ture of the continuous dynamics, such as invariants, symplecticity, and the configuration manifold [5]. Con-ventional numerical integrators construct a discrete approximation to the flow using only information aboutthe vector field, and ignore the physical laws and the geometric properties inherent in the differential equa-tions [17]. Consequently, they do not preserve important characteristics of the dynamics of the continuousequations of motion. In contrast, variational integrators are constructed by discretizing Hamilton’s principle,rather than discretizing the continuous Euler-Lagrange equation [18, 16]. Since they are developed by usinga discrete version of a physical principle, the resulting integrators have the desirable property that they aresymplectic and momentum preserving, and they exhibit good energy behavior for exponentially long times.Geometric numerical integration on S has been studied in [19, 14, 13]. The two-sphere is a homogeneous ∗ TL and ML have been supported in part by NSF Grant DMS-0504747 and DMS-0726263. TL and NHM have beensupported in part by NSF Grant ECS-0244977 and CMS-0555797. † Department of Aerospace Engineering, The University of Michigan, Ann Arbor, MI 48109 ( [email protected] ). ‡ Department of Mathematics, Purdue University, West Lafayette, IN 47907 ( [email protected] ). § Department of Aerospace Engineering, The University of Michigan, Ann Arbor, MI 48109 ( [email protected] ).1
T. Lee, M. Leok, and N. H. McClamroch manifold; the special orthogonal group SO(3) acts transitively on S , and Lie group methods [7] can beadapted to generate numerical flows on S .In this paper, we study Lagrangian mechanical systems on ( S ) n . Thus, it is desirable to preserve thegeometric properties of the dynamics, such as momentum map, symplecticity, and total energy, in additionto the structure of the configuration manifold [11]. We combine the approaches of geometric integrators onhomogeneous manifolds and variational integrators to obtain variational integrators on ( S ) n that preservethe geometric properties of the dynamics as well as the homogeneous structure of the configuration manifold( S ) n concurrently.The contributions of this paper can be summarized in two aspects. In the continuous time setting, theglobal Euler-Lagrange equations on ( S ) n are developed in a compact form without local parameterizationor constraints. This provides insight into the global dynamics on ( S ) n , which is desirable for theoreticalstudies. As a geometric numerical integrator, the discrete Euler-Lagrange equations on ( S ) n are unique in thesense that they conserve both the geometric properties of the dynamics and the manifold structure of ( S ) n simultaneously. The exact geometric properties of the discrete flow not only generate improved qualitativebehavior, but they also provide accurate and reliable computational results in long-time simulation.This paper is organized as follows. Lagrangian mechanics on ( S ) n is described in Section 2. Variationalintegrators on ( S ) n are developed in Section 3. Computational properties are illustrated for several me-chanical systems, namely a double spherical pendulum, an n -body problem on a sphere, an interconnectedsystem of spherical pendula, pure bending of an elastic rod, a spatial array of magnetic dipoles, and moleculardynamics that evolves on a sphere.
2. Lagrangian mechanics on ( S ) n . In this section, continuous equations of motion for a mechanicalsystem defined on ( S ) n are developed in the context of Lagrangian mechanics. It is common in the publishedliterature that the equations of motion are developed by using either two angles or a unit length constraintto characterize S . Any description with two angles has singularities, and any trajectory near a singularityexperiences numerical ill-conditioning. The unit length constraint leads to additional complexity in numericalcomputations. We develop global continuous equations of motion without resorting to local parameterizationsor constraints. To achieve this, it is critical to understand the global characteristics of a mechanical systemon ( S ) n . This section provides a good background for understanding the theory of discrete Lagrangianmechanics on ( S ) n to be introduced in the next section.The two-sphere is the set of points that have the unit length from the origin of R , i.e. S = { q ∈ R | q · q = 1 } . The tangent space T q S for q ∈ S is a plane tangent to the two-sphere at the point q . Thus,a curve q : R → S and its time derivative satisfy q · ˙ q = 0. The time-derivative of a curve can be written as˙ q = ω × q, (2.1)where the angular velocity ω ∈ R is constrained to be orthogonal to q , i.e. q · ω = 0. The time derivativeof the angular velocity is also orthogonal to q , i.e. q · ˙ ω = 0. ( S ) n . We consider a mechanical system on the configurationmanifold S × · · · × S = ( S ) n . We assume that the Lagrangian L : T ( S ) n → R is given by the differencebetween a quadratic kinetic energy and a configuration-dependent potential energy as follows. L ( q , . . . , q n , ˙ q , . . . , ˙ q n ) = 12 n X i,j =1 M ij ˙ q i · ˙ q j − V ( q , . . . , q n ) , (2.2) agrangian Mechanics and Variational Integrators on Two-Spheres q i , ˙ q i ) ∈ T S for i ∈ { , . . . , n } , and M ij ∈ R is the i, j -th element of a symmetric positive definiteinertia matrix M ∈ R n × n for i, j ∈ { , . . . , n } . The configuration dependent potential is denoted by V :( S ) n → R .The action integral is defined as the time integral of the Lagrangian, and the variation of the actionintegral leads to continuous equations of motion by applying Hamilton’s principle. These are standardprocedures to derive the Euler-Lagrange equations. The expression for the infinitesimal variation of q i ∈ S should be carefully developed, since the configuration manifold is not a linear vector space. As in (2.1), theinfinitesimal variation of q i can be written as a vector cross product, δq i = ξ i × q i , (2.3)where ξ i ∈ R is constrained to be orthogonal to q i , i.e. ξ i · q i = 0. From this, the expression for theinfinitesimal variation of ˙ q i is given by δ ˙ q i = ˙ ξ i × q i + ξ i × ˙ q i . (2.4)These expressions are the key elements to obtaining global equations of motion on ( S ) n .The variation of the Lagrangian can be written as δL = n X i,j =1 δ ˙ q i · M ij ˙ q j − n X i =1 δq i · ∂V∂q i , where the symmetric property M ij = M ji is used. Substituting (2.3) and (2.4) into this, and using the vectoridentity ( a × b ) · c = a · ( b × c ) for any a, b, c ∈ R , we obtain δL = n X i,j =1 ˙ ξ i · ( q i × M ij ˙ q j ) + ξ i · ( ˙ q i × M ij ˙ q j ) − n X i =1 ξ i · (cid:18) q i × ∂V∂q i (cid:19) . Let G be the action integral defined as G = R T L ( q , . . . , q n , ˙ q , . . . , ˙ q n ) dt . Using the above equation andintegrating by parts, the variation of the action integral is given by δ G = n X i,j =1 ξ i · ( ˙ q i × M ij ˙ q j + q i × M ij ¨ q j ) (cid:12)(cid:12)(cid:12)(cid:12) T − n X i =1 Z T ξ i · ( q i × n X j =1 M ij ¨ q j ) + q i × ∂V∂q i . From Hamilton’s principle, δ G = 0 for any ξ i vanishing at t = 0 , T . Since ξ i is orthogonal to q i , thecontinuous equations of motion satisfy( q i × n X j =1 M ij ¨ q j ) + q i × ∂V∂q i = c i ( t ) q i (2.5)for some scalar valued functions c i ( t ) for i ∈ { , . . . , n } . Taking the dot product of (2.5) and q i implies that0 = c i ( t ) k q i k = c i ( t ), which is to say that the scalar valued functions are uniformly zero. Now we find anexpression for ¨ q i . Since the left hand side expression is perpendicular to q i , it is zero if and only if its crossproduct with q i is zero. Thus, we obtain q i × ( q i × n X j =1 M ij ¨ q j ) + q i × (cid:18) q i × ∂V∂q i (cid:19) = 0 . (2.6) T. Lee, M. Leok, and N. H. McClamroch
From the vector identity a × ( b × c ) = ( a · c ) b − ( a · b ) c for any a, b, c ∈ R , we have q i × ( q i × ¨ q i ) = ( q i · ¨ q i ) q i − ( q i · q i )¨ q i , = − ( ˙ q i · ˙ q i ) q i − ¨ q i , where we use the properties ddt ( q i · ˙ q i ) = q i · ¨ q i + ˙ q i · ˙ q i = 0 and q i · q i = 1. Substituting these into (2.6), weobtain an expression for ¨ q i , which is summarized as follows. Proposition 1
Consider a mechanical system on ( S ) n whose Lagrangian is expressed as (2.2). The con-tinuous equations of motion are given by M ii ¨ q i = q i × ( q i × n X j =1 j = i M ij ¨ q j ) − ( ˙ q i · ˙ q i ) M ii q i + q i × (cid:18) q i × ∂V∂q i (cid:19) (2.7) for i ∈ { , . . . , n } . Equivalently, this can be written in matrix form as M I × − M ˆ q ˆ q · · · − M n ˆ q ˆ q − M ˆ q ˆ q M I × · · · − M n ˆ q ˆ q ... ... ... − M n ˆ q n ˆ q n − M n ˆ q n ˆ q n · · · M nn I × ¨ q ¨ q ... ¨ q n = − ( ˙ q · ˙ q ) M q + ˆ q ∂V∂q − ( ˙ q · ˙ q ) M q + ˆ q ∂V∂q ... − ( ˙ q n · ˙ q n ) M nn q n + ˆ q n ∂V∂q n , (2.8) where the hat map ˆ · : R → R × is defined such that ˆ ab = a × b for any a, b ∈ R . Since ˙ q i = ω i × q i for the angular velocity ω i satisfying q i · ω i = 0, we have¨ q i = ˙ ω i × q i + ω i × ( ω i × q i ) , = ˙ ω i × q i − ( ω i · ω i ) q i . Substituting this into (2.5) and using the fact that q i · ˙ ω i = 0, we obtain continuous equations of motion interms of the angular velocity. Corollary 1
The continuous equations of motion given by (2.7) can be written in terms of the angularvelocity as M ii ˙ ω i = n X j =1 j = i ( M ij q i × ( q j × ˙ ω j ) + M ij ( ω j · ω j ) q i × q j ) − q i × ∂V∂q i , (2.9)˙ q i = ω i × q i (2.10) for i ∈ { , . . . , n } . Equivalently, this can be written in matrix form as M I × − M ˆ q ˆ q · · · − M n ˆ q ˆ q n − M ˆ q ˆ q M I × · · · − M n ˆ q ˆ q n ... ... ... − M n ˆ q n ˆ q − M n ˆ q n ˆ q · · · M nn I × ˙ ω ˙ ω ... ˙ ω n = P nj =2 M j ( ω j · ω j )ˆ q q j − ˆ q ∂V∂q P nj =1 ,j =2 M j ( ω j · ω j )ˆ q q j − ˆ q ∂V∂q ... P n − j =1 M nj ( ω j · ω j )ˆ q n q j − ˆ q n ∂V∂q n . (2.11) agrangian Mechanics and Variational Integrators on Two-Spheres S ) n . Theyavoid singularities completely, and they preserve the structure of T ( S ) n automatically, if an initial conditionis chosen properly. These equations are useful for understanding global characteristics of the dynamics. Inaddition, these expressions are dramatically more compact than the equations of motion written in terms ofany local parameterization.We need to check that the 3 n × n matrices given by the first terms of (2.8) and (2.11) are nonsingular.This is a property of the mechanical system itself, rather than a consequence of the particular form ofequations of motion. For example, when n = 2, it can be shown thatdet " M I × − M ˆ q ˆ q − M ˆ q ˆ q M I × = det " M I × − M ˆ q ˆ q − M ˆ q ˆ q M I × , = M M ( M M − M ( q · q ) )( M M − M ) . Since the inertia matrix is symmetric positive definite, M , M > M M > M , and from the Cauchy-Schwarz inequality, ( q · q ) ≤ ( q · q )( q · q ) = 1. Thus, the above matrices are non-singular. One mayshow a similar property for n >
2. Throughtout this paper, it is assumed that the 3 n × n matrices givenat the first terms of (2.8) and (2.11) are nonsingular. Under this assumption, the Legendre transformationgiven in the next subsection is a diffeomorphism; the Lagrangian is hyperregular. The Legendre transformation of the Lagrangian gives an equivalentHamiltonian form of equations of motion in terms of conjugate momenta if the Lagrangian is hyperregular.Here, we find expressions for the conjugate momenta, which are used in the following section for the discreteequations of motion. For q i ∈ S , the corresponding conjugate momentum p i lies in the dual space T ∗ q i S .We identify the tangent space T q i S and its dual space T ∗ q i S by using the usual dot product in R . TheLegendre transformation is given by p i · δq i = D ˙ q i L ( q , . . . , q n , ˙ q , . . . , ˙ q n ) · δq i , = n X j =1 M ij ˙ q j · δq i , which is satisfied for any δq i perpendicular to q i . Here D ˙ q i L denotes the derivative of the Lagrangian withrespect to ˙ q i . The momentum p i is an element of the dual space identified with the tangent space, and thecomponent parallel to q i has no effect since δq i · q i = 0. As such, the vector representing p i is perpendicularto q i , and p i is equal to the projection of P nj =1 M ij ˙ q j onto the orthogonal complement to q i , p i = n X j =1 ( M ij ˙ q j − ( q i · M ij ˙ q j ) q i ) = n X j =1 (( q i · q i ) M ij ˙ q j − ( q i · M ij ˙ q j ) q i ) , = M ii ˙ q i − q i × ( q i × n X j =1 j = i M ij ˙ q j ) . (2.12)
3. Variational integrators on ( S ) n . The dynamics of Lagrangian and Hamiltonian systems on ( S ) n have unique geometric properties; the Hamiltonian flow is symplectic, the total energy is conserved in theabsence of non-conservative forces, and the momentum map associated with a symmetry of the system ispreserved. The configuration space is a homogeneous manifold. These geometric features determine thequalitative dynamics of the system, and serve as a basis for theoretical study. T. Lee, M. Leok, and N. H. McClamroch
Conventional numerical integrators construct a discrete approximation of the flow using only informationabout the vector field. Other than the direction specified by the vector field, they completely ignore thephysical laws and the geometric properties inherent in the differential equations [17]. For example, if weintegrate (2.11) by using an explicit Runge-Kutta method, the unit length of the vector q i , and the totalenergy are not preserved numerically; we will see this later in this paper.Numerical integration methods that preserve the simplecticity of a Hamiltonian system have been stud-ied [20]. Coefficients of a Runge-Kutta method can be carefully chosen to satisfy a simplecticity criterionand order conditions to obtain a symplectic Runge-Kutta method. However, it can be difficult to constructsuch integrators, and it is not guaranteed that other invariants of the system, such as the momentum map,are preserved. Alternatively, variational integrators are constructed by discretizing Hamilton’s principle,rather than by discretizing the continuous Euler-Lagrange equation [18, 16]. The key feature of variationalintegrators is that they are derived by a discrete version of a physical principle, so the resulting integra-tors satisfy the physical properties automatically in a discrete sense; they are symplectic and momentumpreserving, and they exhibit good energy behavior for exponentially long times. Lie group methods arenumerical integrators that preserve the Lie group structure of the configuration space [7]. Recently, thesetwo approaches have been unified to obtain Lie group variational integrators that preserve the geometricproperties of the dynamics as well as the Lie group structure of the configuration manifold [10].The two-sphere is a homogeneous manifold. It does not have a Lie group structure by itself, but instead,the special orthogonal group, SO(3) = { F ∈ R × | F T F = I × , det F = 1 } , acts on S in a transitiveway; for any q , q ∈ S , there exist F ∈ SO(3) such that q = F q . If a group acts transitively on amanifold, a curve on the manifold can be represented as the action of a curve in the Lie group on an initialpoint on the manifold. As such, Lie group methods can be applied to obtain numerical integration schemesfor homogeneous manifolds [19, 13, 14]. However, it is not guaranteed that these methods preserve thegeometric properties of the dynamics. In this paper, we focus on a Lagrangian mechanical system evolvingon the homogeneous manifold, ( S ) n by extending the method of Lie group variational integrators [11, 10].The resulting integrator preserves the dynamic characteristics and the homogeneous manifold structureconcurrently. ( S ) n . The procedure to derive discrete Euler-Lagrangeequations follows the development of the continuous time case; the tangent bundle is replaced by a cartesianproduct of the configuration manifold, a discrete Lagrangian is chosen to approximate the integral of theLagrangian over a discrete time step, and the variation of the corresponding discrete action sum providesdiscrete Euler-Lagrange equations, referred to as a variational integrator. The discrete version of the Legendretransformation yields the discrete equations in Hamiltonian form.Let the number of timesteps be N , with constant timesteps h >
0. A variable with subscript k denotesthe value of variable at t = kh . Define a discrete Lagrangian L d : ( S ) n × ( S ) n → R such that it approximatesthe integral of the Lagrangian given by (2.2) over a discrete time step L d ( q k , . . . , q n k , q k +1 , . . . , q n k +1 ) = 12 h n X i,j =1 M ij ( q i k +1 − q i k ) · ( q j k +1 − q j k ) − h V k − h V k +1 , (3.1)where V k denotes the value of the potential at the k -th step, i.e. V k = V ( q k , . . . , q n k ). As given in (2.3),the infinitesimal variation of q i k is written as δq i k = ξ i k × q i k , (3.2) agrangian Mechanics and Variational Integrators on Two-Spheres ξ i k ∈ R is constrained to be orthogonal to q i k , i.e. ξ i k · q i k = 0. The variation of the discreteLagrangian can be written as δL d k = 1 h n X i,j =1 ( δq i k +1 − δq i k ) · M ij ( q j k +1 − q j k ) − h n X i =1 (cid:18) δq i k · ∂V k ∂q i k + δq i k +1 · ∂V k +1 ∂q i k +1 (cid:19) . (3.3)Substituting (3.2) into (3.3), and using the vector identity ( a × b ) · c = a · ( b × c ) for any a, b, c ∈ R , weobtain δL d k = 1 h n X i,j =1 (cid:0) ξ i k +1 · ( q i k +1 × M ij ( q j k +1 − q j k )) − ξ i k · ( q i k × M ij ( q j k +1 − q j k )) (cid:1) − h n X i =1 (cid:18) ξ i k · (cid:18) q i k × ∂V k ∂q i k (cid:19) + ξ i k +1 · (cid:18) q i k +1 × ∂V k +1 ∂q i k +1 (cid:19)(cid:19) . (3.4)Let G d be the discrete action sum defined as G d = P N − k =0 L d k , which approximates the action integral asthe discrete Lagrangian approximates a piece of the action integral over a discrete time step. The variationof the action sum is obtained by using (3.4). Using the fact that ξ i k vanish at k = 0 and k = N , we canreindex the summation, which is the discrete analog of integration by parts, to yield δ G d = N − X k =1 n X i =1 ξ i k · h ( q i k × n X j =1 M ij ( − q j k +1 + 2 q j k − q j k − )) − hq i k × ∂V k ∂q i k . From discrete Hamilton’s principle δ G d = 0 for any ξ i k perpendicular to q i k . Using the same argument givenin (2.5), the discrete equations of motion are given by1 h ( q i k × n X j =1 M ij ( − q j k +1 + 2 q j k − q j k − )) − hq i k × ∂V k ∂q i k = 0 (3.5)for i ∈ { , . . . n } . In addition, we require that the unit length of the vector q i k is preserved. This is achievedby viewing S as a homogeneous manifold. Since the special orthogonal group SO(3) acts on S transitively,we can define a discrete update map for q i k as q i k +1 = F i k q i k for F i k ∈ SO(3). Then, the unit length of the vector q i is preserved through the discrete equations of motion,since q i k +1 · q i k +1 = q Ti k F Ti k F i k q i k = 1. These results are summarized as follows. Proposition 2
Consider a mechanical system on ( S ) n whose Lagrangian is expressed as (2.2). The discreteequations of motion are given by M ii q i k × F i k q i k + q i k × n X j =1 j = i M ij ( F j k − I × ) q j k = q i k × n X j =1 M ij ( q j k − q j k − ) − h q i k × ∂V k ∂q i k = 0 , (3.6) q i k +1 = F i k q i k (3.7) for i ∈ { , . . . n } . For given ( q i k − , q i k ) , we solve (3.6) to obtain F i k ∈ SO(3) . Then, q i k +1 is computed by(3.7). This yields a discrete flow map ( q i k − , q i k ) ( q i k , q i k +1 ) , and this process is repeated. T. Lee, M. Leok, and N. H. McClamroch
We find discrete equations of motion in terms of the an-gular velocity. The discrete Legendre transformation is given as follows [16]. p i k · δq i k = − D q ik L d k · δq i k , = h n X j =1 M ij ( q j k +1 − q j k ) + h ∂V k ∂q i k · δq i k , which can be directly obtained from (3.3). This is satisfied for any δq i k perpendicular to q i k . Using the sameargument used to derive (2.12), the conjugate momenta p i k is the projection of the expression in bracketsonto the orthogonal complement of q i k . Thus, we obtain q i k = − h q i k × ( q i k × n X j =1 M ij ( q j k +1 − q j k )) − h q i k × (cid:18) q i k × ∂V k ∂q i k (cid:19) . Comparing this to (2.12), substituting ˙ q i k = ω i k × q i k , and rearranging, we obtain q j k × M ii ω i k + ( q i k × n X j =1 j = i M ij ( ω j k × q j k )) − h ( q i k × n X j =1 M ij ( q j k +1 − q j k )) − h q i k × ∂V k ∂q i k = 0 . Since the expression in the brackets is orthogonal to q i k , the left hand side is equal to zero if and only if theexpression in the brackets is zero. Thus, M ii ω i k + ( q i k × n X j =1 j = i M ij ( ω j k × q j k )) = 1 h ( q i k × n X j =1 M ij ( q j k +1 − q j k )) + h q i k × ∂V k ∂q i k . (3.8)This provides a relationship between ( q i k , ω i k ) and ( q i k , q i k +1 ). Comparing this with (3.5), we obtain M ii ω i k + ( q i k × n X j =1 j = i M ij ( ω j k × q j k )) = 1 h ( q i k × n X j =1 M ij ( q j k − q j k − )) − h q i k × ∂V k ∂q i k , (3.9)which provides a relationship between ( q i k , ω i k ) and ( q i k − , q i k ). Equations (3.8) and (3.9) give a discrete flowmap in terms of the angular velocity; for a given ( q i k , ω i k ), we find ( q i k , q i k +1 ) by using (3.8). Substitutingthis into (3.9) expressed at the k + 1th step, we obtain ( q i k +1 , ω i k +1 ). This procedure is summarized asfollows. Corollary 2
The discrete equations of motion given by (3.6) and (3.7) can be written in terms of the angularvelocity as M ii q i k × F i k q i k + q i k × n X j =1 j = i M ij ( F j k − I × ) q j k = M ii hω i k − ( q i k × n X j =1 j = i M ij ( q j k × hω j k )) − h q i k × ∂V k ∂q i k , (3.10) q i k +1 = F i k q i k , (3.11) M ii ω i k +1 − ( q i k +1 × n X j =1 j = i M ij ( q j k +1 × ω j k +1 )) = 1 h ( q i k +1 × n X j =1 M ij ( q j k +1 − q j k )) − h q i k +1 × ∂V k +1 ∂q i k +1 (3.12) agrangian Mechanics and Variational Integrators on Two-Spheres for i ∈ { , . . . , n } . Equivalently, (3.12) can be written in a matrix form as M I × − M ˆ q k +1 ˆ q k +1 · · · − M n ˆ q ˆ q n k +1 − M ˆ q k +1 ˆ q k +1 M I × · · · − M n ˆ q k +1 ˆ q n k +1 ... ... ... − M n ˆ q n k +1 ˆ q k +1 − M n ˆ q n k +1 ˆ q k +1 · · · M nn I × ω k +1 ω k +1 ... ω n k +1 = h ( q k +1 × P nj =1 M j ( q j k +1 − q j k )) − h q k +1 × ∂V k +1 ∂q k +1 h ( q k +1 × P nj =1 M j ( q j k +1 − q j k )) − h q k +1 × ∂V k +1 ∂q k +1 ... h ( q n k +1 × P nj =1 M nj ( q j k +1 − q j k )) − h q n k +1 × ∂V k +1 ∂q nk +1 . (3.13) For a given ( q i k , ω i k ) , we solve (3.10) to obtain F i k ∈ SO(3) . Then, q i k +1 and ω i k +1 are computed by(3.11) and (3.13), respectively. This yields a discrete flow map in terms of the angular velocity ( q i k , ω i k ) ( q i k +1 , ω i k +1 ) , and this process is repeated. For the discrete equations of motion, we need to solve (3.6) and(3.10) to obtain F i k ∈ SO(3). Here we present a computational approach. The implicit equations given by(3.6) and (3.10) have the following structure. M ii q i × F i q i + q i × n X j =1 j = i M ij ( F j − I × ) q j = d i (3.14)for i ∈ { , . . . , n } , where M ij ∈ R , q i ∈ S , d i ∈ R are known, and we need to find F i ∈ SO(3). We derive anequivalent equation in terms of local coordinates for F i . This is reasonable since F i represents the relativeupdate between two integration steps. Using the Cayley transformation [21], F i ∈ SO(3) can be expressedin terms of f i ∈ R as F i = ( I × + ˆ f i )( I × − ˆ f i ) − , = 11 + f i · f i ((1 − f i · f i ) I × + 2 f i f Ti + 2 ˆ f i ) . The operation F i q i can be considered as a rotation of the vector q i about the direction f i with rotation angle2 tan − k f i k . Since the rotation of the vector q i about the direction q i has no effect, we can assume that f i is orthogonal to q i , i.e. f i · q i = 0. Under this assumption, F i q i is given by F i q i = 11 + f i · f i ((1 − f i · f i ) q i + 2 ˆ f i q i ) . (3.15)Thus, we obtain q i × F i q i = 21 + f i · f i q i × ( f i × q i ) = 21 + f i · f i f i , ( F j − I × ) q j = −
21 + f j · f j ( q j f Tj + ˆ q j ) f j , T. Lee, M. Leok, and N. H. McClamroch where we use the property, ˆ q i f i = q i × f i = − ˆ f i q i . Substituting these into (3.14), we obtain M I × f · f − M ˆ q (ˆ q + q f T )1+ f · f · · · − M n ˆ q (ˆ q n + q n f Tn )1+ f n · f n − M ˆ q (ˆ q + q f T )1+ f · f M I × f · f · · · − M n ˆ q (ˆ q n + q n f Tn )1+ f n · f n ... ... ... − M n ˆ q n (ˆ q + q f T )1+ f · f − M n ˆ q n (ˆ q + q f T )1+ f · f · · · M nn I × f n · f n f f ... f n = d d ... d n , (3.16)which is an equation equivalent to (3.14), written in terms of local coordinates for F i using the Cayleytransformation. Any numerical method to solve nonlinear equations can be applied to find f i . Then, F i q i iscomputed by using (3.15). In particular, (3.16) is written in a form that can be readily applied to a fixedpoint iteration method [8].If there are no coupling terms in the kinetic energy, we can obtain an explicit solution of (3.14). When M ij = 0 for i = j , (3.16) reduces to 2 M ii f i · f i f i = d i . Using the identity, θ θ = sin 2 θ for any θ ∈ R , it can be shown that the solution of this equation is givenby f i = tan (cid:0) sin − ( k d i k /M ii ) (cid:1) d i k d i k . Substituting this into (3.15) and rearranging, we obtain F i q i = d i M ii × q i + − (cid:13)(cid:13)(cid:13)(cid:13) d i M ii (cid:13)(cid:13)(cid:13)(cid:13) ! / q i . Using this expression, we can rewrite the discrete equations of motion given in (3.10)–(3.13) in an explicitform.
Corollary 3
Consider a mechanical system on ( S ) n whose Lagrangian is expressed as (2.2) where M ij = 0 for i = j , i.e. the dynamics are coupled only though the potential energy. The explicit discrete equations ofmotion are given by q i k +1 = (cid:18) hω i k − h M ii q i k × ∂V k ∂q i k (cid:19) × q i k + − (cid:13)(cid:13)(cid:13)(cid:13) hω i k − h M ii q i k × ∂V k ∂q i k (cid:13)(cid:13)(cid:13)(cid:13) ! / q i k , (3.17) ω i k +1 = ω i k − h M ii q i k × ∂V k ∂q i k − h M ii q i k +1 × ∂V k +1 ∂q i k +1 (3.18) for i ∈ { , . . . , n } . ( S ) n . Since variational integrators are derived fromthe discrete Hamilton’s principle, they are symplectic, and momentum preserving. The discrete action sumcan be considered as a zero-form on ( S ) n × ( S ) n which maps the initial condition of a discrete flow satisfyingthe discrete Euler-Lagrange equation to the action sum for that trajectory. The simplecticity of the discreteflow follows from the fact the the iterated exterior derivative of any differential form is zero. If the discreteLagrangian exhibits a symmetry, the corresponding momentum map is preserved since by symmetry, thevariation of the discrete Lagrangian in the symmetry direction is zero, which in combination with the discreteEuler–Lagrange equations, implies a discrete version of Noether’s theorem. Detailed proofs for the symplecticproperty and the momentum preserving property can be found in [16]. The total energy oscillates around agrangian Mechanics and Variational Integrators on Two-Spheres S ) n without need oflocal parameterizations, explicit constraints or reprojection. Using the characteristics of the homogeneousmanifold, the discrete update map is represented by a group action of SO(3), and a proper subspace issearched to obtain a compact, possibly explicit, form for the numerical integrator. As a result, the followingnumerical problems are avoided: (i) local parameterizations yield singularities; (ii) numerical trajectoriesin the vicinity of a singularity experiences numerical ill-conditioning; (iii) unit length constraints lead toadditional computational complexity; (iv) reprojection corrupts the numerical accuracy of trajectories [5, 13].It can be shown that these variational integrators have second-order accuracy as the discrete action sumis a second-order approximation of the action integral. Higher-order integrators can be easily constructedby applying a symmetric composition method [23]. The computational properties of variational integrators on ( S ) n andexplicit Runge-Kutta methods are compared for several mechanical systems taken from variety of scientificareas, namely a double spherical pendulum, an n -body problem on a sphere, an interconnected system ofspherical pendula, pure bending of an elastic rod, a spatial array of magnetic dipoles, and molecular dynamicsthat evolves on a sphere. Example 1 (Double Spherical Pendulum)
A double spherical pendulum is defined by two mass par-ticles serially connected to frictionless two degree-of-freedom pivots by rigid massless links acting under auniform gravitational potential. The dynamics of a double spherical pendulum has been studied in [15], anda variational integrator is developed in [22] by explicitly using unit length constraints.Let the mass and the length of the pendulum be m , m , l , l ∈ R , respectively, and let e = [0 , , ∈ R be the direction of gravity. The vector q ∈ S represents the direction from the pivot to the first mass,and the vector q ∈ S represents the direction from the first mass to the second mass. The inertia matrixis given by M = ( m + m ) l , M = m l l , and M = m l . The gravitational potential is written as V ( q , q ) = − ( m + m ) gl e · q − m gl e · q for the gravitational acceleration g ∈ R . Substituting theseinto (2.10)–(2.11), the continuous equations of motion for the double spherical pendulum are given by˙ q = ω × q ˙ q = ω × q , (3.19) " ( m + m ) l I × − m l l ˆ q ˆ q − m l l ˆ q ˆ q m l I × ˙ ω ˙ ω = " m l l ( ω · ω )ˆ q q + ( m + m ) gl ˆ q e m l l ( ω · ω )ˆ q q + m gl ˆ q e , (3.20)which are more compact than existing equations written in terms of angles. Another nice property is thatthe same structure for the equations of motion is maintained for n >
2. Thus, it is easy to generalize theseequations of motion to a triple, or more generally, a multiple-link spherical pendulum.We compare the computational properties of the discrete equations of motion given by (3.10)–(3.13)with a 4(5)-th order variable step size Runge-Kutta method for (3.19)–(3.20). We choose m = m = 1 kg, l = l = 9 .
81 m. The initial conditions are q = [0 . , , . q = [0 , , ω = [ − . , , . ω = [0 , ,
0] rad / sec. The simulation time is 100 sec, and the step-size of the discrete equations of motion is h = 0 .
01. Figure 3.1 shows the computed total energy and the configuration manifold errors. The variationalintegrator preserves the total energy and the structure of ( S ) n well for this chaotic motion of the doublespherical pendulum. The mean total energy variation is 2 . × − Nm, and the mean unit length error is8 . × − . But, there is a notable increase of the computed total energy for the Runge-Kutta method,2 T. Lee, M. Leok, and N. H. McClamroch(a) Trajectory of pendulum (b) Computed total energy (c) Unit length error
Fig. 3.1 . Numerical simulation of a double spherical pendulum (RK45: blue, dotted, VI: red, solid) where the mean variation of the total energy is 7 . × − Nm. The Runge-Kutta method also fails topreserve the structure of ( S ) n . The mean unit length error is 6 . × − . Example 2 ( n -body Problem on Sphere) An n -body problem on the two-sphere deals with the motionof n mass particles constrained to lie on a two-sphere, acting under a mutual potential. Let m i ∈ R and q i ∈ S be the mass and the position vector of the i -th particle, respectively. The i, j -th element of the inertiamatrix is M ij = m i when i = j , and M ij = 0 otherwise. In [9], the following expression for the potential isintroduced as an analog of a gravitational potential, V ( q , . . . , q n ) = − γ n X i,j =1 i = j q i · q j p − ( q i · q j ) for a constant γ . Substituting these into (2.7), the continuous equations of motion for the n -body problemon a sphere are given by m i ¨ q i = − m i ( ˙ q i · ˙ q i ) q i − q i × (cid:0) q i × γ n X j =1 j = i q j (1 − ( q i · q j ) ) / (cid:1) (3.21)for i ∈ { , . . . , n } .A two-body problem on the two-sphere under this gravitational potential is studied in [6] by explicitlyusing unit length constraints. Here we study a three-body problem, n = 3. Since there are no coupling termsin the kinetic energy, we use the explicit form of the variational integrator. We compare the computationalproperties of the discrete equations of motion given by (3.17)–(3.18) with a 2-nd order fixed step size Runge-Kutta method for (3.21). We choose m = m = m = 1, and γ = 1. The initial conditions are q =[0 , − , q = [0 , , q = [ − , , ω = [0 , , − . ω = [1 , , ω = [0 , , . × − to1 . × − when the step size is reduced by 10 times from 10 − to 10 − , which verifies the second orderaccuracy numerically. agrangian Mechanics and Variational Integrators on Two-Spheres (a) Trajectory of particles −4 −3 −2 −1 −6 −5 −4 −3 −2 −1 Step size T o t a l E n e r gy E rr o r (b) Total energy error v.s. step size −4 −3 −2 −1 −15 −10 −5 Step size U n it l e ng t h E rr o r (c) Unit length error v.s. step size Fig. 3.2 . Numerical simulation of a 3-body problem on sphere (RK2: blue, square, VI: red, circle)
Example 3 (Interconnection of Spherical Pendula)
We study the dynamics of n spherical pendulaconnected by linear springs. Each pendulum is a mass particle connected to a frictionless two degree-of-freedom pivot by a rigid massless link acting under a uniform gravitational potential. It is assumed that allof the pivot points lie on a horizontal plane, and some pairs of pendulua are connected by linear springs atthe centers of links.Let the mass and the length of the i -th pendulum be m i , l i ∈ R , respectively. The vector q i ∈ S represents the direction from the i -th pivot to the i -th mass. The inertia matrix is given by M ij = m i l i when i = j , and M ij = 0 otherwise. Let Ξ be a set defined such that ( i, j ) ∈ Ξ if the i -th pendulum andthe j -th pendulum are connected. For a connected pair ( i, j ) ∈ Ξ, define κ ij ∈ R and r ij ∈ R as thecorresponding spring constant and the vector from the i -th pivot to the j -th pivot, respectively. The basesfor the inertial frame are chosen such that the direction along gravity is denoted by e = [0 , , ∈ R , andthe horizontal plane is spanned by e = [0 , , , e = [0 , , ∈ R . The potential energy is given by V ( q , . . . q n ) = − n X i =1 m i gl i q i · e + X ( i,j ) ∈ Ξ κ ij (cid:18)(cid:13)(cid:13)(cid:13)(cid:13) r ij + 12 l j q j − l i q i (cid:13)(cid:13)(cid:13)(cid:13) − k r ij k (cid:19) . Substituting these into (2.9)–(2.10), the continuous equations of motion for the interconnection of sphericalpendula are given by m i l i ˙ ω i = − q i × ∂V∂q i (3.22)˙ q i = ω i × q i (3.23)for i ∈ { , . . . , n } .We compare the computational properties of the discrete equations of motion given by (3.17)–(3.18)with a 2-nd order fixed step size explicit Runge-Kutta method for (3.22)–(3.23), and the same Runge-Kuttamethod with reprojection; at each time step, the vectors q i k are projected onto S by using normalization.We choose four interconnected pendula, n = 4, and we assume each pendulum has the same mass andlength; m i = 0 . l i = 0 . { (1 , , (2 , , (3 , , (4 , } , and thecorresponding spring constants and the relative vector between pivots are given by κ = 10, κ = 20, κ = 30, κ = 40 N / m, r = − r = l i e , and r = − r = − l i e . The initial conditions are chosen as q = q = q = e , q = [0 . , . , . ω = [ − , , ω = ω = ω = 0 rad / sec4 T. Lee, M. Leok, and N. H. McClamroch(a) Motion of pendula (b) Computed total energy (c) Unit length error
Fig. 3.3 . Numerical simulation of a system of 4 spherical pendula (RK2: blue, dotted, RK2 with projection: black, dashed,VI: red, solid)
Figure 3.3 shows the computed total energy and the unit length errors. The variational integratorpreserves the total energy and the structure of ( S ) n well. The mean total energy variation is 3 . × − Nm, and the mean unit length error is 4 . × − . For both Runge-Kutta methods, there is anotable increase of the computed total energy. It is interesting to see that the reprojection approach makesthe total energy error worse, even though it preserves the structure of ( S ) n accurately. This shows that astandard reprojection method can corrupt numerical trajectories [5, 13]. Example 4 (Pure Bending of Elastic Rod)
We study the dynamics of ( n + 1) rigid rod elements thatare serially connected by rotational springs, where the ‘zeroth’ rod is assumed to be fixed to a wall. Thus,the configuration space is ( S ) n . This can be considered as a simplified dynamics model for pure non-planarbending of a thin elastic rod that is clamped at one end and free at the other end. Notably, this approachis geometrically exact, and preserves the length of the elastic rod in the presence of large displacements.The mass and the length of the i -th rod element are denoted by m i , l i ∈ R , respectively. The inertiamatrix is given by M ii = 13 m i l i + n X k = i +1 m i l i , M ij = n X k =max { i,j } m k l k for i, j ∈ { , . . . n } and i = j . The potential energy is composed of gravitational terms and elastic bendingterms given by V ( q , . . . , q n ) = − n X i =1 m i g (cid:0) i − X j =1 l j q j + 12 l i q i (cid:1) · e + 12 κ i (1 − q i − · q i ) , where a constant vector q ∈ S denotes the direction of the zeroth rod element fixed to a wall, and κ i ∈ R denotes spring constants. The bases for the inertial frame are chosen such that the gravity direction isdenoted by e = [0 , , ∈ R , and the horizontal plane is spanned by e = [0 , , , e = [0 , , ∈ R .Suppose that the total mass and length of rod are given by m, l , and each rod element has the same massand length, i.e. m i = mn +1 , l i = ln +1 for i ∈ { , . . . , n } . Substituting these into (2.8), the continuous agrangian Mechanics and Variational Integrators on Two-Spheres (a) Deformation of rod t T o t a l E n e r gy (b) Computed total energy (c) Unit length error Fig. 3.4 . Numerical simulation of an elastic rod (RK45: blue, dotted, VI: red, solid) equations of motion for the pure bending of an elastic rod are given by n − / n +1) ml I × − n − n +1) ml ˆ q ˆ q · · · − n +1) ml ˆ q ˆ q − n − n +1) ml ˆ q ˆ q n − / n +1) ml I × · · · − n +1) ml ˆ q ˆ q ... ... ... − n +1) ml ˆ q n ˆ q n − n +1) ml ˆ q n ˆ q n · · · / n +1) ml I × ¨ q ¨ q ...¨ q n = − n − / n +1) ml ( ˙ q · ˙ q ) q + ˆ q ∂V∂q − n − / n +1) ml ( ˙ q · ˙ q ) q + ˆ q ∂V∂q ... − / n +1) ml ( ˙ q n · ˙ q n ) q n + ˆ q n ∂V∂q n , (3.24)We compare the computational properties of the discrete equations of motion given by (3.10)–(3.13)with a 4(5)-th order variable step size Runge-Kutta method for (3.24). We choose 10 rod elements, n = 10,and the total mass and the total length are m = 55 g, l = 1 . , m. The spring constants are chosen as κ i = 1000 Nm. Initially, the rod is aligned horizontally; q i = e for all i ∈ , . . . n . The initial angularvelocity for each rod element is zero except ω = [0 , ,
10] rad / sec. This represents the dynamics of the rodafter an initial impact. The simulation time is 3 sec, and the step-size of the discrete equations of motion is h = 0 . S ) n . The mean total energy variation is 1 . × − Nm,and the mean unit length error is 2 . × − . There is a notable dissipation of the computed totalenergy for the Runge-Kutta method, where the mean variation of the total energy is 3 . × − Nm. TheRunge-Kutta method also fail to preserve the structure of ( S ) n . The mean unit length error is 1 . × − . Example 5 (Spatial Array of Magnetic Dipoles)
We study dynamics of n magnetic dipoles uniformlydistributed on a plane. Each magnetic dipole is modeled as a spherical compass; a thin rod magnet supportedby a frictionless, two degree-of-freedom pivot acting under their mutual magnetic field. This can be consideredas a simplified model for the dynamics of micromagnetic particles [3].The mass and the length of the i -th magnet are denoted by m i , l i ∈ R , respectively. The magneticdipole moment of the i -th magnet is denoted by ν i q i , where ν i ∈ R is the constant magnitude of themagnetic moment measured in ampere square-meters, and q i ∈ S is the direction of the north pole fromthe pivot point. Thus, the configuration space is ( S ) n . The inertia matrix is given by M ij = m i l i when i = j , and M ij = 0 otherwise. Let r ij ∈ R be the vector from the i -pivot point to the j -th pivot point. The6 T. Lee, M. Leok, and N. H. McClamroch(a) Motion of magnetic dipoles −3 t T o t a l E n e r gy (b) Computed total energy (c) Unit length error Fig. 3.5 . Numerical simulation of an array of magnetic dipoles (RK45: blue, dotted, VI: red, solid) mutual potential energy of the array of magnetic dipoles are given by V ( q , . . . , q n ) = 12 n X i,j =1 j = i µ ν i ν j π k r ij k (cid:20) ( q i · q j ) − k r ij k ( q i · r ij )( q j · r ij ) (cid:21) , where µ = 4 π × − N · A − is the permeability constant. Substituting these into (2.9)–(2.10), the continuousequations of motion for the spatial array of magnetic dipoles are given by112 m i l i ˙ ω i = − q i × n X j =1 j = i µ ν i ν j π k r ij k (cid:20) q j − k r ij k r ij ( q j · r ij ) (cid:21) , (3.25)˙ q i = ω i × q i (3.26)for i ∈ { , . . . , n } .We compare the computational properties of the discrete equations of motion given by (3.17)–(3.18) witha 4(5)-th order variable step size Runge-Kutta method for (3.25)–(3.26). We choose 16 magnetic dipoles, n = 16, and we assume each magnetic dipole has the same mass, length, and magnitude of magnetic moment; m i = 0 .
05 kg, l i = 0 .
02 m, ν i = 0 . · m . The magnetic dipoles are located at vertices of a 4 × . l i . The initial conditions are chosen as q i = [1 , , ω i = [0 , ,
0] for all i ∈ { , . . . , } except q = [0 . , . , − . ω = [0 , . ,
0] rad / sec.Figure 3.5 shows the computed total energy and the unit length errors. The variational integratorpreserves the total energy and the structure of ( S ) n well. The mean total energy variation is 8 . × − Nm, and the mean unit length error is 1 . × − . There is a notable dissipation of the computedtotal energy for the Runge-Kutta method, where the mean variation of the total energy is 2 . × − Nm.The Runge-Kutta method also fail to preserve the structure of ( S ) n . The mean unit length error is 1 . × − . Example 6 (Molecular Dynamics on a Sphere)
We study molecular dynamics on S . Each moleculeis modeled as a particle moving on S . Molecules are subject to two distinct forces: an attractive force atlong range and a repulsive force at short range. Let m i ∈ R and q i ∈ S be the mass and the positionvector of the i -th molecule, respectively. The i, j -th element of the inertia matrix is M ij = m i when i = j ,and M ij = 0 otherwise. The Lennard-Jones potential is a simple mathematical model that represents the agrangian Mechanics and Variational Integrators on Two-Spheres (a) Initial trajectories t T o t a l E n e r gy (b) Computed total energy Fig. 3.6 . Numerical simulation of molecular dynamics on a sphere (a) t = 0 (b) t = 0 .
25 (c) t = 0 . t = 0 .
75 (e) t = 5 Fig. 3.7 . Kinetic energy distributions over time behavior of molecules [12] V ( q , . . . , q n ) = 12 n X i,j =1 j = i ǫ "(cid:18) σ k q i − q j k (cid:19) − (cid:18) σ k q i − q j k (cid:19) , where the first term models repulsion between molecules at short distance according to the Pauli principle,and the second term models attraction at long distance generated by van der Walls forces. The constant ǫ and σ are molecular constants; ǫ is proportional to the strength of the mutual potential, and σ characterizeinter-molecular force. Substituting these into (2.7), the continuous equations of motion for the moleculardynamics on a sphere are given by m i ¨ q i = − m i ( ˙ q i · ˙ q i ) q i − q i × (cid:0) q i × n X j =1 j = i ǫ q i − q j k q i − q j k " σ k q i − q j k − σ k q i − q j k (3.27)for i ∈ { , . . . , n } .We choose 642 molecules, n = 642, and we assume each molecule has the same mass, m i = 1. Initially,molecules are uniformly distributed on a sphere. The strength of the potential is chosen as ǫ = 0 .
01, andthe constant σ is chosen such that the inter-molecular force between neighboring molecules is close to zero.The initial velocities are modeled as two vortices separated by 30 degrees. The simulation time is 5 sec, andthe step size is h = 0 . . × − , and the mean unit length error is 5 . × − . In molecular dynamicssimulations, macroscopic quantities such as temperature and pressure are more useful than trajectories8 T. Lee, M. Leok, and N. H. McClamroch of molecules. Figure 3.7 shows the change of kinetic energy distributions over time, which measures thetemperature [1]; the sphere is discretized by an icosahedron with 5120 triangular faces, and the color ofa face is determined by the average kinetic energy for molecules that lie within the face and within itsneighboring faces. The local kinetic energy is represented by color shading of blue, green, yellow, and redcolors in ascending order.
4. Conclusions.
Euler-Lagrange equations and variational integrators are developed for Lagrangianmechanical systems evolving on ( S ) n where the Lagrangian is written in a particular form given by (2.2).The structure of S is carefully considered to obtain global equations of motion on ( S ) n without localparameterizations or explicit constraints.In the continuous time setting, this provides a remarkably compact form of the equations of motioncompared to the popular angular description. For example, it is not practical to study a triple sphericalpendulum by using angles due to the complexity of the trigonometric expressions involved. The globalEuler-Lagrange equations on ( S ) n maintain the same compact structure for arbitrary n . In particular, it ispossible to use them as a finite element model for a continuum problem as shown in Example 4. They arealso useful for the theoretical study of global dynamic characteristics.The variational integrators on ( S ) n preserve the geometric properties of the dynamics as well as thestructure of the configuration manifold concurrently. They are symplectic, momentum preserving, and theyexhibit good energy behavior for exponentially long time as they are derived from discrete Hamilton’s prin-ciple. Using the characteristics of the homogeneous manifold ( S ) n , the discrete update map is representedby a group action of SO(3) to obtain compactly represented, and possibly explicit, numerical integrators. Inparticular, variational integrators on ( S ) n completely avoid the singularities and complexity introduced bylocal parameterizations and explicit constraints. REFERENCES[1]
M. P. Allen and D. J. Tildesley , Computer Simulation of Liquids , Clarendon Press, 1987.[2]
S. Bendersky and B. Sandler , Investigation of spatial double pendulum: an engineering approach , Discrete Dynamicsin Nature and Society, 2006 (2006), pp. 1–12.[3]
X. Z. Cheng, M. B. A. Jalil, and H. K. Lee , Time-quantified monte carlo algorithm for interacting spin array micro-magnetic dynamics , Physical Review B, 73 (2006), p. 224438.[4]
E. Hairer , Backward analysis of numerical integrators and symplectic methods , Annals of Numerical Mathematics, 1(1994), pp. 107–132. Scientific computation and differential equations (Auckland, 1993).[5]
E. Hairer, C. Lubich, and G. Wanner , Geometric Numerical Integration , Springer, 2000.[6] ,
Geometric numerical integration illustrated by the St¨ormer-Verlet method , Acta Numer., 12 (2003), pp. 399–450.[7]
A. Iserles, H. Munthe-Kaas, S. P. Nørsett, and A. Zanna , Lie-group methods , in Acta Numerica, vol. 9, CambridgeUniversity Press, 2000, pp. 215–365.[8]
C. T. Kelley , Iterative Methods for Linear and Nonlinear Equations , SIAM, 1995.[9]
V. V. Kozlov and A. O. Harin , Kepler’s problem in constant curvature spaces , Celestial Mechanics and DynamicalAstronomy, 54 (1992), pp. 393–399.[10]
T. Lee, M. Leok, and N. H. McClamroch , Lie group variational integrators for the full body problem , ComputerMethods in Applied Mechanics and Engineering, 196 (2007), pp. 2907–2924.[11] ,
Lie group variational integrators for the full body problem in orbital mechanics , Celestial Mechanics and DynamicalAstronomy, 98 (2007), pp. 121–144.[12]
J. E. Lennard-Jones , Cohesion , The Proceedings of the Physical Society, 43 (1931), pp. 461–482.[13]
D. Lewis and N. Nigam , Geometric integration on spheres and some interesting applications , Journal of Computationaland Applied Mathematics, 151 (2003), pp. 141–170.[14]
D. Lewis and P. J. Olver , Geometric integration algorithms on homogeneous manifolds , Foundations of ComputationalMathematics, 2 (2001), pp. 363–392.agrangian Mechanics and Variational Integrators on Two-Spheres [15] J. E. Marsden, E. Scheurle, and J. Wendlandt , Visualization of orbits and pattern evocation for the double sphericalpendulum , in International Congress on Industrial and Applied Mathematics, vol. 87, 1995.[16]
J. E. Marsden and M. West , Discrete mechanics and variational integrators , in Acta Numerica, vol. 10, CambridgeUniversity Press, 2001, pp. 317–514.[17]
R. McLachlan and R. Quispel , Six Lectures on The Geometric Integration of ODEs, in Foundations of ComputationalMathematics , London Mathematical Society Lecture Note, 284, Cambridge University Press, 2001, pp. 155–210.[18]
J. Moser and A. P. Veselov , Discrete versions of some classical integrable systems and factorization of matrix polyno-mials , Communications in Mathematical Physics, 139 (1991), pp. 217–243.[19]
H. Munthe-Kaas and A. Zanna , Numerical Integration of Differential Equations on Homogeneous Manifolds, in Foun-dations of Computational Mathematics , Springer, 1997, pp. 305–315.[20]
J. M. Sanz-Serna , Symplectic integrators for hamiltonian problems: an overview , Acta Numerica, (1992), pp. 243–286.[21]
M. D. Shuster , Survey of attitude representations , Journal of the Astronautical Sciences, 41 (1993), pp. 439–517.[22]
J. M. Wendlandt and J. E. Marsden , Mechanical integrators derived from a discrete variational principle , Physica D,106 (1997), pp. 223–246.[23]