Lie Group integrators for mechanical systems
Elena Celledoni, Ergys ?okaj, Andrea Leone, Davide Murari, Brynjulf Owren
LLie Group integrators for mechanical systems
Elena Celledoni Ergys C¸ okaj Andrea LeoneDavide Murari Brynjulf OwrenNovember 2020, last compiled: February 26, 2021
Abstract
Since they were introduced in the 1990s, Lie group integrators havebecome a method of choice in many application areas. These includemultibody dynamics, shape analysis, data science, image registration andbiophysical simulations. Two important classes of intrinsic Lie group inte-grators are the Runge–Kutta–Munthe–Kaas methods and the commutatorfree Lie group integrators. We give a short introduction to these classes ofmethods. The Hamiltonian framework is attractive for many mechanicalproblems, and in particular we shall consider Lie group integrators forproblems on cotangent bundles of Lie groups where a number of differ-ent formulations are possible. There is a natural symplectic structure onsuch manifolds and through variational principles one may derive sym-plectic Lie group integrators. We also consider the practical aspects ofthe implementation of Lie group integrators, such as adaptive time step-ping. The theory is illustrated by applying the methods to two nontrivialapplications in mechanics. One is the N-fold spherical pendulum wherewe introduce the restriction of the adjoint action of the group SE (3) to T S , the tangent bundle of the two-dimensional sphere. Finally, we showhow Lie group integrators can be applied to model the controlled pathof a payload being transported by two rotors. This problem is modeledon R × ( SO (3) × so (3)) × ( T S ) and put in a format where Lie groupintegrators can be applied. In many physical problems, including multi-body dynamics, the configurationspace is not a linear space, but rather consists of a collection of rotations andtranslations. A simple example is the free rigid body whose configuration spaceconsists of rotations in 3D. A more advanced example is the simplified modelof the human body, where the skeleton at a given time is described as a systemof interacting rods and joints. Mathematically, the structure of such problemsis usually best described as a manifold. Since manifolds by definition can beequipped with local coordinates, one can always describe and simulate suchsystems locally as if they were linear spaces. There are of course many choicesof local coordinates, for rotations some famous ones are: Euler angles, the1 a r X i v : . [ m a t h . NA ] F e b ait-Bryan angles commonly used in aerospace applications, the unit lengthquaternions, and the exponentiated skew-symmetric 3 × α -method for linear spaces proposedby Hilber, Hughes and Taylor [21] has been popular for solving problems inmultibody dynamics, and in [1, 2] this method is generalised to a Lie groupintegrator.The literature on Lie group integrators is rich and diverse, the interestedreader may consult the surveys [25, 9, 6, 43] and Chapter 4 of the monograph[17] for further details.In this paper we discuss different ways of applying Lie group integratorsto simulating the dynamics of mechanical multi-body systems. Our point ofdeparture is the formulation of the models as differential equations on manifolds.Assuming to be given either a Lie group acting transitively on the manifold M or a set of frame vector fields on M , we use them to describe the mechanicalsystem and further to build the numerical integrator. We shall here mostlyconsider schemes of the types commonly known as Crouch–Grossman methods[10], Runge–Kutta–Munthe–Kaas methods [38, 39] and Commutator-free Liegroup methods [5]. 2he choice of Lie group action is often not unique and thus the same me-chanical system can be described in different equivalent ways. Under numericaldiscretization the different formulations can lead to the conservation of differentgeometric properties of the mechanical system. In particular, we explore theeffect of these different formulations on a selection of examples in multi-bodydynamics. Lie group integrators have been succesfully applied for the simula-tion of mechanical systems, and in problems of control, bio-mechnics and otherengineering applications, see for example [45], [26] [8], [24]. The present work ismotivated by applications in modeling and simulation of slender structures likeCosserat rods and beams [48], and one of the examples presented here is theapplication to a chain of pendula. Another example considers an applicationfor the controlled dynamics of a multibody system.In section 2 we give a review of the methods using only the essential in-trinsic tools of Lie group integrators. The algorithms are simple and amenablefor a coordinate-free description suited to object oriented implementations. Insection 3, we discuss Hamiltonian systems on Lie groups, and we present threedifferent Lie group formulations of the heavy top equations. These systems (andtheir Lagrangian counterpart) often arise in applications as building blocks ofmore realistic systems which comprise also damping and control forces. In sec-tion 4, we discuss some ways of adapting the integration step size in time. Insection 5 we consider the application to a chain of pendula. And in section 6we consider the application of a multi-body system of interest in the simulationand control of drone dynamics. Lie group integrators solve differential equations whose solution evolve on amanifold M . This means that we seek a curve y ( t ) ∈ M whose tangent at anypoint coincides with a vector field F ∈ X ( M ) and passing through a designatedinitial value y at t = t ˙ y ( t ) = F | y ( t ) , y ( t ) = y . (1)Before addressing numerical methods for solving (1) it is necessary to introducea convenient way of representing the vector field F . There are different waysof doing this. One is to furnish M with a transitive action ψ : G × M → M by some Lie group G of dimension d ≥ dim M . We denote the action of g on m as g · m , i.e. g · m = ψ ( g, m ). Let g be the Lie algebra of G , and denoteby exp : g → G the exponential map. We define ψ ∗ : g → X ( M ) to be theinfinitesimal generator of the action, i.e. F ξ | m = ψ ∗ ( ξ ) | m = ddt (cid:12)(cid:12)(cid:12)(cid:12) t =0 ψ (exp( tξ ) , m ) (2)The transitivity of the action now ensures that ψ ∗ ( g ) | m = T m M for any m ∈ M ,such that any tangent vector v m ∈ T m M can be represented as v m = ψ ∗ ( ξ v ) | m ξ v ∈ g ( ξ v may not be unique). Consequently, for any vector field F ∈ X ( M ) there exists a map f : M → g such that F | m = ψ ∗ ( f ( m )) | m , for all m ∈ M (3)This is the original tool [39] for representing a vector field on a manifold witha group action. Another approach was used in [10] where a set of frame vectorfields E , . . . , E d in X ( M ) was introduced assuming that for every m ∈ M ,span { E | m , . . . , E d | m } = T m M . Then, for any vector field F ∈ X ( M ) there are, in general non-unique, functions f i : M → R such that F | m = d (cid:88) i =1 f i ( m ) E i | m . A fixed vector ξ ∈ R d will define a vector field F ξ on M similar to (2) F ξ | m = d (cid:88) i =1 ξ i E i | m (4)If ξ i = f i ( p ) for some p ∈ M , the corresponding F ξ will be a vector field in thelinear span of the frame which coincides with F at the point p . Such a vectorfield was named by [10] as a the vector field frozen at p .The two formulations just presented are in many cases connected, and canthen be used in an equivalent manner. Suppose that e , . . . , e d is a basis of theLie algebra g , then we can simply define frame vector fields as E i = ψ ∗ ( e i ) andthe vector field we aim to describe is, F | m = ψ ∗ ( f ( m )) | m = ψ ∗ ( (cid:88) i f i ( m ) e i ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) m = (cid:88) i f i E i | m . As mentioned above there is a non-uniqueness issue when defining a vector fieldby means of a group action or a frame. A more fundamental description canbe obtained using the machinery of connections. The assumption is that thesimply connected manifold M is equipped with a connection which is flat andhas constant torsion. Then F p , the frozen vector field of F at p defined above,can be defined as the unique element F p ∈ X ( M ) satisfying1. F p | p = F | p ∇ X F p = 0 for any X ∈ X ( M ).So F p is the vector field that coincides with F at p and is parallel transportedto any other point on M by the connection ∇ . Since the connection is flat, theparallel transport from the point p to another point m ∈ M does not dependon the chosen path between the two points. For further details, see e.g. [31].4 xample 1. For mechanical systems on Lie groups, two important construc-tions are the adjoint and coadjoint representatons. For every g ∈ G there is anautomorphism Ad g : g → g defined as Ad g ( ξ ) = T L g ◦ T R g − ( ξ ) where L g and R g are the left and right multiplications respectively, L g ( h ) = gh and R g ( h ) = hg . Since Ad is a representation, i.e. Ad gh = Ad g ◦ Ad h it alsodefines a left Lie group action by G on g . From this definition and a dualitypairing (cid:104)· , ·(cid:105) between g and g ∗ , we can also derive a representation on g ∗ denoted Ad ∗ g , simply by (cid:104) Ad ∗ g ( µ ) , ξ (cid:105) = (cid:104) µ, Ad g ( ξ ) (cid:105) , ξ ∈ g , µ ∈ g ∗ . The action g · µ = Ad ∗ g − ( µ ) has infinitesimal generator given as ψ ∗ ( ξ ) | µ = − ad ∗ ξ µ Following [33], for a Hamiltonian H : T ∗ G → R , define H − to be its restrictionto g ∗ . Then the Lie-Poisson reduction of the dynamical system is defined on g ∗ as ˙ µ = − ad ∗ δH − δµ µ and this vector field is precisely of the form (3) with f ( µ ) = δH − δµ ( µ ) . A sideeffect of this is that the integral curves of these Lie-Poisson systems preservecoadjoint orbits, making the coadjoint action an attractive choice for Lie groupintegrators.Let us now detail the situation for the very simple case where G = SO (3) .The Lie algebra so (3) can be modeled as × skew-symmetric matrices, andvia the standard basis we identify each such matrix ˆ ξ by a vector ξ ∈ R , thisidentification is known as the hat map ˆ ξ = − ξ ξ ξ − ξ − ξ ξ (5) Now, we also write the elements of so (3) ∗ as vectors in R with duality pairing (cid:104) µ, ξ (cid:105) = µ T ξ . With these representations, we find that the coadjoint action canbe expressed as g · µ = ψ ( g, µ ) = Ad ∗ g − µ = gµ the rightmost expression being a simple matrix-vector multiplication. Since g is orthogonal, it follows that the coadjoint orbits foliate 3-space into sphericalshells, and the coadjoint action is transitive on each of these orbits. The freerigid body can be cast as a problem on T SO (3) ∗ with a left invariant Hamiltonianwhich reduces to the function H − ( µ ) = 12 (cid:104) µ, I − µ (cid:105) n so (3) ∗ where I : so (3) → so (3) ∗ is the inertia tensor. From this, we can nowset f ( µ ) = δH − /δµ = I − µ . We the recover the Euler free rigid body equationas ˙ µ = ψ ∗ ( f ( µ ) | µ = − ad ∗ I − µ µ = − I − µ × µ where the last expression involves the cross product of vectors in R . The simplest numerical integrator for linear spaces is the explicit Euler method.Given an initial value problem ˙ y = F ( y ), y (0) = y the method is defined as y n +1 = y n + hF ( y n ) for some stepsize h . In the spirit of the previous section,one could think of the Euler method as the h -flow of the constant vector field F y n ( y ) = F ( y n ), that is y n +1 = exp( hF y n ) y n This definition of the Euler method makes sense also when F is replaced bya vector field on some manifold. In this general situation it is known as theLie–Euler method.We shall here consider the two classes of methods known as Runge–Kutta–Munthe–Kaas (RKMK) methods and Commutator-free Lie group methods.For RKMK methods the underlying idea is to transform the problem fromthe manifold M to the Lie algebra g , take a time step, and map the result backto M . The transformation we use is y ( t ) = exp( σ ( t )) · y , σ (0) = 0 . The transformed differential equation for σ ( t ) makes use of the derivative ofthe exponential mapping, the reader should consult [39] for details about thederivation, we give the final result˙ σ ( t ) = dexp − σ ( t ) ( f (exp( σ ( t )) · y )) (6)The map v (cid:55)→ dexp u ( v ) is linear and invertible when u belongs to some suffi-ciently small neighborhood of 0 ∈ g . It has an expansion in nested commutators[20]. Using the operator ad u ( v ) = [ u, v ] and its powers ad u v = [ u, [ u, v ]] etc, onecan writedexp u ( v ) = e z − z (cid:12)(cid:12)(cid:12)(cid:12) z =ad u ( v ) = v + 12 [ u, v ] + 16 [ u, [ u, v ]] + · · · (7)and the inverse isdexp − u ( v ) = ze z − (cid:12)(cid:12)(cid:12)(cid:12) z =ad u ( v ) = v −
12 [ u, v ] + 112 [ u, [ u, v ]] + · · · (8)The RKMK methods are now obtained simply by applying some standardRunge–Kutta method to the transformed equation (6) with a time step h , usinginitial value σ (0) = 0. This leads to an output σ ∈ g and one simply sets6 = exp( σ ) · y . Then one repeats the procedure replacing y by y in thenext step etc. While solving (6) one needs to evaluate dexp − u ( v ) as a part ofthe process. This can be done by truncating the series (8) since σ (0) = 0 im-plies that we always evaluate dexp − u with u = O ( h ), and thus, the k th iteratedcommutator ad ku = O ( h k ). For a given Runge–Kutta method, there are someclever tricks that can be done to minimise the total number of commutatorsto be included from the expansion of dexp − u v , see [4, 40]. We give here oneconcrete example of an RKMK method proposed in [4] f n, = hf ( y n ) ,f n, = hf (exp( f n, ) · y n ) ,f n, = hf (exp( f n, − [ f n, , f n, ]) · y n ) ,f n, = hf (exp( f n, ) · y n ) ,y n +1 = exp( ( f n, + 2 f n, + 2 f n, + f n, − [ f n, , f n, ])) · y n . The other option is to compute the exact expression for dexp − u ( v ) for theparticular Lie algebra we use. For instance, it was shown in [7] that for the Liealgebra so (3) one hasdexp − u ( v ) = v − u × v + α − (1 − α cot α ) u × ( u × v )We will present the corresponding formula for se (3) in Section 2.3.The second class of Lie group integrators to be considered here are thecommutator-free methods, named this way in [5] to emphasize the contrastto RKMK schemes which usually include commutators in the method format.These schemes include the Crouch-Grossman methods [10] and they have theformat Y n,r = exp (cid:32) h (cid:88) k α kr,J f n,k (cid:33) · · · exp (cid:32) h (cid:88) k α kr, f n,k (cid:33) · y n f n,r = f ( Y n,r ) y n +1 = exp (cid:32) h (cid:88) k β kJ f n,k (cid:33) · · · exp (cid:32) h (cid:88) k β k f n,k (cid:33) · y n Here the Runge–Kutta coefficients α kr,j , β rj are related to a classical Runge–Kutta scheme with coefficients a kr , b r in that a kr = (cid:80) j α kr,j and b r = (cid:80) j β rj .The α kr,j , β rj are usually chosen to obtain computationally inexpensive schemeswith the highest possible order of convergence. The computational complexityof the above schemes depends on the cost of computing an exponential as wellas of evaluating the vector field. Therefore it makes sense to keep the number ofexponentials J in each stage as low as possible, and possibly also the number ofstages s . A trick proposed in [5] was to select coefficients that make it possiblereuse exponential from one stage to another. This is perhaps best illustrated7hrough the following example from [5], a generalisation of the classical 4thorder Runge–Kutta method. Y n, = y n Y n, = exp( hf n, ) · y n Y n, = exp( hf n, ) · y n Y n, = exp( hf n, − hf n, ) · Y n, y n + = exp( h (3 f n, + 2 f n, + 2 f n, − f n, )) · y n y n +1 = exp( h ( − f n, + 2 f n, + 2 f n, + 3 f n, )) · y n + (9)where f n,i = f ( Y n,i ). Here, we see that one exponential is saved in computing Y n, by making use of Y n, . dexp − u ( v ) in se (3) As an alternative to using a truncated version of the infinite series for dexp − u (8), one can consider exact expressions obtained for certain Lie algebras. Since se (3) is particularly important in applications to mechanics, we give here itsexact expression. For this, we represent elements of se (3) as a pair ( A, a ) ∈ R × R ∼ = R , the first component corresponding to a skew-symmetric matrixˆ A via (5) and a is the translational part. Now, let ϕ ( z ) be a real analyticfunction at z = 0. We define ϕ + ( z ) = ϕ ( iz ) + ϕ ( − iz )2 , ϕ − ( z ) = ϕ ( iz ) − ϕ ( − iz )2 i We next define the four functions g ( z ) = ϕ − ( z ) z , ˜ g ( z ) = g (cid:48) ( z ) z , g ( z ) = ϕ (0) − ϕ + ( z ) z , ˜ g ( z ) = g (cid:48) ( z ) z and the two scalars ρ = A T a , α = (cid:107) A (cid:107) . One can show that for any ( A, a ) and(
B, b ) in se (3), it holds that ϕ (ad ( A,a ) )( B, b ) = (
C, c )where C = ϕ (0) B + g ( α ) A × B + g ( α ) A × ( A × B ) c = ϕ (0) b + g ( α ) ( a × B + A × b ) + ρ ˜ g ( α ) A × B + ρ ˜ g ( α ) A × ( A × B )+ g ( α ) ( a × ( A × B ) + A × ( a × B ) + A × ( A × b ))Considering for instance (8), we may now use ϕ ( z ) = ze z − to calculate g ( z ) = − , ˜ g ( z ) = 0 , g ( z ) = 1 − z cot z z , ˜ g ( z ) = 1 z ddz g ( z ) , ϕ (0) = 1 . and thereby obtain an expression for dexp − A,a ) ( B, b ) with the formula above.8imilar types of formulas are known for computing the matrix exponentialas well as functions of the ad-operator for several other Lie groups of small andmedium dimension. For instance in [37] a variety of coordinate mappings forrigid body motions are discussed. For Lie algebras of larger dimension, boththe exponential mapping and dexp − u may become computationally infeasible.For these cases, one may benefit from replacing the exponential by some othercoordinate map for the Lie group φ : g → G . One option is to use canonicalcoordinates of the second kind [44]. Then for some Lie groups such as theorthogonal, unitary and symplectic groups, there exist other maps that can beused and which are computationally less expensive. A popular choice is theCayley transformation [12]. In this section we consider Hamiltonian systems on Lie groups. These sys-tems (and their Lagrangian counterpart) often appear in mechanics applicationsas building blocks for more realistic systems with additional damping and con-trol forces. We consider canonical systems on the cotangent bundle of a Liegroup and Lie-Poisson systems which can arise by symmetry reduction or oth-erwise. We illustrate the various cases with different formulations of the heavytop system.
The coadjoint action by G on g ∗ is denoted Ad ∗ g defined for any g ∈ G as (cid:104) Ad ∗ g µ, ξ (cid:105) = (cid:104) µ, Ad g ξ (cid:105) , ∀ ξ ∈ g , (10)where Ad : g → g is the adjoint representation and for a duality pairing (cid:104)· , ·(cid:105) between g ∗ and g .We consider the cotangent bundle of a Lie group G , T ∗ G and identify itwith G × g ∗ using the right multiplication R g : G → G and its tangent map-ping R g ∗ := T R g . The cartesian product G × g ∗ can be given a semi-directproduct structure that turns it into a Lie group G := G × g ∗ where the groupmultiplication is ( g , µ ) · ( g , µ ) = ( g · g , µ + Ad ∗ g − µ ) . (11)Acting by left multiplication any vector field F ∈ X ( G ) is expressed bymeans of a map f : G → T e G , F ( g, µ ) = T e R ( g,µ ) f ( g, µ ) = ( R g ∗ f , f − ad ∗ f µ ) , (12)where f = f ( g, µ ) ∈ g , f = f ( g, µ ) ∈ g ∗ are the two components of f .9 .2 Symplectic form and Hamiltonian vector fields The right trivialised symplectic form pulled back to G reads ω ( g,µ ) (( R g ∗ ξ , δν ) , ( R g ∗ ξ , δν )) = (cid:104) δν , ξ (cid:105) − (cid:104) δν , ξ (cid:105) − (cid:104) µ, [ ξ , ξ ] (cid:105) , ξ , ξ ∈ g . (13) See [30] for more details, proofs and for a the left trivialized symplectic form.The vector field F is a Hamiltonian vector field if it satisfiesi F ω = dH, for some Hamiltonian function H : T ∗ G → R , where i F is defined as i F ( X ) := ω ( F, X ) for any vector field X . This implies that the map f for such a Hamil-tonian vector field gets the form f ( g, µ ) = (cid:18) ∂H∂µ ( g, µ ) , − R ∗ g ∂H∂g ( g, µ ) (cid:19) . (14)The following is a one-parameter family of symplectic Lie group integratorson T ∗ G : ( ξ, ¯ n ) = hf (cid:0) exp( θξ ) · g , dexp ∗− ξ µ + (1 − θ )dexp ∗− (1 − θ ) ξ ¯ n (cid:1) , (15)( g , µ ) = (exp( ξ ) , Ad ∗ exp( − (1 − θ ) ξ ) ¯ n ) · ( g , µ ) . (16)For higher order integrators of this type and a complete treatment see [3]. A mechanical system formulated on the cotangent bundle T ∗ G with a left orright invariant Hamiltonian can be reduced to a system on g ∗ [32]. In fact fora Hamiltonian H right invariant under the left action of G , ∂H∂g = 0, and from(12) and (14) we get for the second equation˙ µ = ∓ ad ∗ ∂H∂µ µ, (17)where the positive sign is used in case of left invariance (see e.g. section 13.4 in[34]). The solution to this system preserves coadjoint orbits, thus using the Liegroup action g · µ = Ad ∗ g − µ, to build a Lie group integrator results in preservation of such coadjoint orbits.Lie group integrators for this interesting case were studied in [14].The Lagrangian counterpart to these Hamiltonian equations are the theEuler–Poincar´e equations , [23]. ω ( g,µ ) is obtained from the natural symplectic form on T ∗ G (which is a differential two-form), defined as Ω ( g,p g ) (( δv , δπ ) , ( δv , δπ )) = (cid:104) δπ , δv (cid:105) − (cid:104) δπ , δv (cid:105) , by right trivialization. The Euler–Poincar´e equations are Euler–Lagrange equations with respect to a Lagrange–d’Alembert principle obtained taking constraint variations. .4 Three different formulations of the heavy top equa-tions The heavy top is a simple test example for illustrating the behaviour of Liegroup methods. We will consider three different formulations for this mechanicalsystem. The first formulation is on T ∗ SO (3) where the equations are canonicalHamiltonian, a second point of view is that the system is a Lie–Poisson systemon se (3) ∗ , and finally it is canonical Hamiltonian on a larger group with aquadratic Hamiltonian function. The three different formulations suggest theuse of different Lie group integrators. T ∗ SO (3) . The heavy top is a rigid body with a fixed point in a gravitational field. Thephase space of this mechanical system is T ∗ SO (3) where the equations of theheavy top are in canonical Hamiltonian form. Assuming ( Q, p ) are coordinatesfor T ∗ SO (3), Π = ( T e L Q ) ∗ ( p ) is the left trivialized or body momentum. TheHamiltonian of the heavy top is given in terms of ( Q, Π) as H : SO (3) (cid:110) so (3) ∗ → R , H ( Q, Π) = 12 (cid:104) Π , I − Π (cid:105) + M g(cid:96) Γ · X , Γ = Q − Γ , where I : so (3) → so (3) ∗ is the inertia tensor, here represented as a diagonal 3 × Γ = Q − Γ , where Γ ∈ R is the axis of the spatial coordinate systemparallel to the direction of gravity but pointing upwards, M is the mass of thebody, g is the gravitational acceleration, X is the body fixed unit vector of theoriented line segment pointing from the fixed point to the center of mass of thebody, (cid:96) is the length of this segment. The equations of motion on SO (3) (cid:110) so (3) ∗ are ˙Π = Π × I − Π +
M g(cid:96) Γ × X , (18)˙ Q = Q (cid:91)I − Π . (19)The identification of T ∗ SO (3) with SO (3) (cid:110) so (3) ∗ via right trivializationleads to the spatial momentim variable π = ( T e R Q ) ∗ ( p ) = Q Π. The equationswritten in the space variables (
Q, π ) get the form˙ π = M g(cid:96) Γ × Q X , (20)˙ Q = ˆ ωQ ω = Q I − Q T π. (21)where, the first equation states that the component of π parallel to Γ is constantin time. These equations can be obtained from (12) and (14) on the righttrivialized T ∗ SO (3), SO (3) (cid:110) so (3) ∗ , with the heavy top Hamiltonian and thesymplectic Lie group integrators (15)-(16) can be applied in this case. Similarmethods where proposed in [30] and [47].11 .4.2 Heavy top equations on se ∗ (3)The Hamiltonian of the heavy top is not invariant under the action of SO (3),so the equations (18)-(19)given in section (3.4.1) cannot be reduced to so ∗ (3),nevertheless the heavy top equations are Lie–Poisson on se ∗ (3), [51, 16, 46].Observe that the equations of the heavy top on T ∗ SO (3) (18)-(19) can beeasily modified eliminating the variable Q ∈ SO (3) and replacing it with Γ ∈ R Γ = Q − Γ to obtain ˙Π = Π × I − Π +
M g(cid:96) Γ × X , (22)˙ Γ = Γ × ( I − Π) . (23)We will see that the solutions of these equations evolve on se ∗ (3). In whatfollows, we consider elements of se ∗ (3) to be pairs of vectors in R , e.g. (Π , Γ ).Correspondingly the elements of SE (3) are represented as pairs ( g, u ) with g ∈ SO (3) and u ∈ R . The group multiplication in SE (3) is then( g , u ) · ( g , u ) = ( g g , g u ) , where g g is the product in SO (3) and g u is the product of a 3 × R . The coadjoint representation and its infinitesimalgenerator on se ∗ (3) take the formAd ∗ ( g, u ) (Π , Γ ) = ( g − (Π − u × Γ ) , g − Γ ) , ad ∗ ( ξ, u ) (Π , Γ ) = ( − ξ × Π − u × Γ , − ξ × Γ ) . Using this expression for ad ∗ ( ξ, u ) with ( ξ = ∂H∂ Π , u = ∂H∂ Γ ), it can be easily seenthat the equations (17) in this setting reproduce the heavy top equations (22)-(23). Therefore the equations are Lie–Poisson equations on se ∗ (3). Howeversince the heavy top is a rigid body with a fixed point and there are no transla-tions, these equations do not arise from a reduction of T ∗ SE (3). Moreover theHamiltonian on se (3) ∗ is not quadratic and the equations are not geodesic equa-tions. Implicit and explicit Lie group integrators applicable to this formulationof the heavy top equations and preserving coadjoint orbits were discussed in[14], for a variable step-size integrator applied to this formulation of the heavytop see [11]. We rewrite the heavy top equations one more time considering the constantvector p = − M g(cid:96) X as a momentum variable conjugate to the position q ∈ R and where p = Q − Γ + ˙ q , and the Hamiltonian is a quadratic function of Π, Q , p and q : H : T ∗ SO (3) × R ∗ × R → R ,H ((Π , Q ) , ( p , q )) = 12 (cid:104) Π , I − Π (cid:105) + 12 (cid:107) p − Q − Γ (cid:107) − (cid:107) Q − Γ (cid:107) , see [22, section 8.5]. This Hamiltonian is invariant under the left action of SO (3). The corresponding equations are canonical on T ∗ S ≡ S (cid:110) s ∗ where12 = SO (3) × R with Lie algebra s := so (3) × R and T ∗ S can be identifiedwith T ∗ SO (3) × R ∗ × R . The equations are˙Π = Π × I − Π − ( Q − Γ ) × p , (24)˙ Q = Q (cid:91)I − Π , (25)˙ p = 0 , (26)˙ q = p − Q − Γ . (27)and in the spatial momentum variables˙ π = − Γ × Q p , (28)˙ Q = ˆ ωQ, ω = Q I − Q T π, (29)˙ p = 0 , (30)˙ q = p − Q − Γ . (31)Similar formulations were considered in [29] for the stability analysis of an un-derwater vehicle. One approach for varying the step size is based on the use of an embeddedRunge–Kutta pair. This principle can be carried from standard Runge–Kuttamethods in vector spaces to the present situation with RKMK and commutator-free schemes via minor modifications. We briefly summarise the main principleof embedded pairs before giving more specific details for the case of Lie groupintegrators. This approach is very well documented in the literature and goesback to Merson [35] and a detailed treatment can be found in [18, p. 165–168].An embedded pair consists of a main method used to propagate the numer-ical solution, together with some auxiliary method that is only used to obtainan estimate of the local error. This local error estimate is in turn used to derivea step size adjustment formula that attempts to keep the local error estimateapproximately equal to some user defined tolerance tol in every step. Supposethe main method is of order p and the auxiliary method is of order ˆ p (cid:54) = p. Both methods are applied to the input value y n and yields approximations y n +1 and ˆ y n +1 respectively, using the same step size h n +1 . Now, some distancemeasure between y n +1 and ˆ y n +1 provides an estimate e n +1 for the size of thelocal truncation error. Thus, e n +1 = Ch ˆ p +1 n +1 + O ( h ˆ p +2 ). Aiming at e n +1 ≈ tolin every step, one may use a formula of the type h n +1 = θ (cid:18) tol e n (cid:19) p +1 h n (32) In this paper we will assume ˆ p < p in which case the local error estimate is relevant forthe approximation ˆ y n +1 θ is a ‘safety factor’, typically chosen between 0 . .
9. In case thestep is rejected because e n > tol we can redo the step with a step size obtainedby the same formula. We summarise the approach in the following algorithmGiven y n , h n , tolLet h := h n repeat Compute y n +1 , ˆ y n +1 , e n +1 from y n , h Update stepsize h := θ (cid:16) tol e n +1 (cid:17) α h accepted := e n +1 < tol if acceptedupdate step index: n := n + 1 h n := h until acceptedHere we have used again the safety factor θ , and the parameter α is generallychosen as α = p, ˆ p ) . We need to specify how to calculate the quantity e n +1 in each step. For RKMKmethod the situation is simplified by the fact that we are solving the local prob-lem (6) in the linear space g , where the known theory can be applied directly. Soany standard embedded pair of Runge–Kutta methods described by coefficients( a ij , b i , ˆ a ij , ˆ b i ) of orders ( p, ˆ p ) can be applied to the full dexpinv-equation (6) toobtain local Lie algebra approximations σ , ˆ σ and one use e.g. e n = (cid:107) σ − ˆ σ (cid:107) (note that the equation itself depends on y n ). For methods which use a trun-cated version of the series for dexp − u one may also try to optimise performanceby including commutators that are shared between the main method and theauxiliary scheme. For the commutator-free methods of section 2.2 the situation is different sincesuch methods do not have a natural local representation in a linear space. Onecan still derive embedded pairs, and this can be achieved by studying orderconditions [42] as was done in [11]. Now one obtains after each step two ap-proximations y n +1 and ˆ y n +1 on M both by using the same initial value y n and step size h n +1 . One must also have access to some metric d to calculate e n +1 = d ( y n +1 , ˆ y n +1 ) We give a few examples of embedded pairs. ( p, ˆ p ) = (3 , Y n, = y n Y n, = exp( hf n, ) · y n Y n, = exp( hf n, ) · y n y n +1 = exp( h ( − f n, + f n, )) · Y n, ˆ y n +1 = exp( h ( f n, + f n, )) · y n One could also have reused the third stage Y n, in the update, rather than Y n, . Y n, = y n Y n, = exp( hf n, ) · y n Y n, = exp( h ( f n, + f n, ) · y n y n +1 = exp( h ( − f n, − f n, + f n, )) · Y n, ˆ y n +1 = exp( h ( f n, + 3 f n, )) · y n It is always understood that the frozen vector fields are f n,i := f Y n,i . (4 , ,
3) was derived, in the sensethat one attempted to minimise the number of stages and exponential in theembedded pair as a whole. This came, however, at the expense of a relativelylarge error constant. So rather than presenting the method from that paper,we suggest a simpler procedure at the cost of some more computational workper step, we simply furnish the commutator-free method of section 2 by a thirdorder auxiliary scheme. It can be described as follows:1. Compute Y n,i , i = 1 . . . , y n +1 from (9)2. Compute an additional stage ¯ Y n, and then ˆ y n +1 as¯ Y n, = exp( hf n, ) · y n ˆ y n +1 = exp( h ( − f n, + 3 f n, + 4 ¯ f n, )) · exp( h f n, ) · y n (33) N -fold 3D pendulum In this section, we present a model for a system of N connected 3-dimensionalpendulums. The modelling part comes from [27], and here we study the vectorfield describing the dynamics, in order to re-frame it into the Lie group inte-grators setting described in the previous sections. The model we use is not15ompletely realistic since, for example, it neglects possible interactions betweenpendulums, and it assumes ideal spherical joints between them. However, thisis still a relevant example from the point of view of geometric numerical inte-gration. More precisely, we show a possible way to work with a configurationmanifold which is not a Lie group, applying the theoretical instruments intro-duced before. The Lagrangian we consider is a function from ( T S ) N to R . Figure 1: 3 − fold pendulum at a fixed time instant, with fixed point placed at the origin. Instead of the coordinates ( q , ..., q N , ˙ q , ..., ˙ q N ), where ˙ q i ∈ T q i S , we choose towork with the angular velocities. Precisely, T q i S = { v ∈ R : v T q i = 0 } = (cid:104) q i (cid:105) ⊥ ⊂ R , and hence for any ˙ q i ∈ T q i S there exist ω i ∈ R such that ˙ q i = ω i × q i , whichcan be interpreted as the angular velocity of q i . So we can assume withoutloss of generality that ω Ti q i = 0 (i.e. ω i ∈ T q i S ) and pass to the coordinates( q , ω , q , ω , ..., q N , ω N ) ∈ ( T S ) N to describe the dynamics. In this sectionwe denote with m , ..., m N the masses of the pendulums and with L , ..., L N their lengths. Figure 1 shows the case N = 3. We organize the section intothree parts:1. We define the transitive Lie group action used to integrate this modelnumerically,2. We show a possible way to express the dynamics in terms of the infinites-imal generator of this action, for the general case of N joint pendulums,3. We focus on the case N = 2, as a particular example. For this setting,we present some numerical experiment comparing various Lie group inte-grators and some classical numerical integrator. Then we conclude withnumerical experiments on variable step size.16 .1 Transitive group action on ( T S ) N We characterize a transitive action for (
T S ) N , starting with the case N = 1 andgeneralizing it to N > se (3), the Lie algebra of SE (3), and R . We start from the Ad-actionof SE (3) on se (3) (see [22]), which writesAd : SE (3) × se (3) → se (3) , Ad((
R, r ) , ( u, v )) = ( Ru, Rv + ˆ rRu ) . Since se (3) (cid:39) R , the Ad-action allows us to define the following Lie groupaction on R ψ : SE (3) × R → R , ψ (( R, r ) , ( u, v )) = ( Ru, Rv + ˆ rRu ) . We can think of ψ as a Lie group action on T S since, for any q ∈ R , it mapspoints of T S | q | := { (˜ q, ˜ ω ) ∈ R × R : ˜ ω T ˜ q = 0 , | ˜ q | = | q |} ⊂ R into other points of T S | q | . Moreover, with standard arguments (see [41]), it ispossible to prove that the orbit of a generic point m = ( q, ω ) ∈ R with ω T q = 0coincides with O ( m ) = T S | q | . In particular, when q ∈ R is a unit vector (i.e. q ∈ S ), ψ allows us to define atransitive Lie group action on T S = T S | q | =1 which writes ψ : SE (3) × T S → T S ψ (( A, a ) , ( q, ω )) := ψ ( A,a ) ( q, ω ) = ( Aq, Aω + ˆ aAq ) = (¯ q, ¯ ω ) . To conclude the description of the action, we report here its infinitesimal gen-erator which is fundamental in the Lie group integrators setting ψ ∗ (( u, v )) | ( q,ω ) = (ˆ uq, ˆ uω + ˆ vq ) . We can extend this construction to the case
N > SE (3)and equipped with the direct product structure. More precisely, we considerthe group G = ( SE (3)) N and by direct product structure we mean that for anypair of elements δ (1) = ( δ (1)1 , ..., δ (1) N ) , δ (2) = ( δ (2)1 , ..., δ (2) N ) ∈ G, denoted with ∗ the semidirect product of SE (3), we define the product ◦ on G as δ (1) ◦ δ (2) := ( δ (1)1 ∗ δ (2)1 , ..., δ (1) N ∗ δ (2) N ) ∈ G. N = 1 to larger N s as follows ψ : ( SE (3)) N × ( T S ) N → ( T S ) N ,ψ (( A , a , ..., A N , a n ) , ( q , ω , ..., q N , ω N )) == ( A q , A ω + ˆ a A q , ..., A N q N , A N ω N + ˆ a N A N q N ) , whose infinitesimal generator writes ψ ∗ ( ξ ) | m = (ˆ u q , ˆ u ω + ˆ v q , ..., ˆ u N q N , ˆ u N ω N + ˆ v N q N ) , where ξ = [ u , v , ..., u N , v N ] ∈ se (3) N and m = ( q , ω , ..., q N , ω N ) ∈ ( T S ) N .We have now the only group action we need to deal with the N − fold sphericalpendulum. In the following part of this section we work on the vector fielddescribing the dynamics and adapt it to the Lie group integrators setting. We consider the vector field F ∈ X (( T S ) N ), describing the dynamics of the N -fold 3D pendulum, and we express it in terms of the infinitesimal generator of theaction defined above. More precisely, we find a function F : ( T S ) N → se (3) N such that ψ ∗ ( f ( m )) | m = F | m , ∀ m ∈ ( T S ) N . We omit the derivation of F starting from the Lagrangian of the system, whichcan be found in the section devoted to mechanical systems on ( S ) N of [27]. Theconfiguration manifold of the system is ( S ) N , while the Lagrangian, expressedin terms of the variables ( q , ω , ..., q N , ω N ) ∈ ( T S ) N , writes L ( q, ω ) = T ( q, ω ) − U ( q ) = 12 N (cid:88) i,j =1 (cid:16) M ij ω Ti ˆ q Ti ˆ q j ω j (cid:17) − N (cid:88) i =1 (cid:16) N (cid:88) j = i m j (cid:17) gL i e T q i , where M ij = (cid:16) N (cid:88) k = max { i,j } m k (cid:17) L i L j I ∈ R × is the inertia matrix of the system and e = [0 , , T . Noticing that when i = j we get ω Ti ˆ q Ti ˆ q i ω i = ω Ti ( I − q i q Ti ) ω i = ω Ti ω i , we simplify the notation writing T ( q, ω ) = 12 N (cid:88) i,j =1 (cid:16) ω Ti R ( q ) ij ω j (cid:17) R ( q ) ∈ R N × N is a symmetric block matrix defined as R ( q ) ii = (cid:16) N (cid:88) j = i m j (cid:17) L i I ∈ R × ,R ( q ) ij = (cid:16) N (cid:88) k = j m k (cid:17) L i L j ˆ q Ti ˆ q j ∈ R × = R ( q ) Tji , i < j.
The vector field on which we need to work defines the following first-order ODE˙ q i = ω i × q i , i = 1 , ..., N,R ( q ) ˙ ω = N (cid:88) j =1 j (cid:54) = i M ij | ω j | ˆ q i q j − (cid:16) N (cid:88) j = i m j (cid:17) gL i ˆ q i e i =1 ,...,N ∈ R N By direct computation it is possible to see that, for any q = ( q , ..., q N ) ∈ ( S ) N and ω ∈ T q S × ... × T q N S , we have( R ( q ) ω ) i ∈ T q i S . Therefore, there is a well-defined linear map A q : T q S × ... × T q N S → T q S × ... × T q N S , A q ( ω ) := R ( q ) ω. We can even notice that R ( q ) defines a positive-definite bilinear form on thislinear space, since ω T R ( q ) ω = N (cid:88) i,j =1 ω Ti ˆ q Ti M ij ˆ q j ω j = N (cid:88) i,j =1 (ˆ q i ω i ) T M ij (ˆ q j ω j ) = v T M v > . The last inequality holds because M is the inertia matrix of the system and henceit defines a symmetric positive-definite bilinear form on T q S × ... × T q N S , seee.g. [15] . This implies the map A q is invertible and hence we are ready toexpress the vector field in terms of the infinitesimal generator. We can rewritethe ODEs for the angular velocities as follows:˙ ω = A − q (cid:16) [ g , ..., g N ] T (cid:17) = h ( q, ω ) ...h N ( q, ω ) = a ( q, ω ) × q ...a N ( q, ω ) × q N It follows from the definition of the inertia tensor, i.e.0 ≤ ˜ T ( q, ˙ q ) = 12 N (cid:88) i =1 (cid:16) (cid:88) j ≥ i m j (cid:17) L i L j ˙ q Ti ˙ q j := 12 ˙ q T M ˙ q. Moreover, in this situation it is even possible to explicitly find the Cholesky factorization thematrix M with an iterative algorithm. g i = g i ( q, ω ) = N (cid:88) j =1 j (cid:54) = i M ( q ) ij | ω j | ˆ q i q j − (cid:16) N (cid:88) j = i m j (cid:17) gL i ˆ q i e , i = 1 , ..., N and a , ..., a N : ( T S ) N → R are N functions whose existence is guaranteedby the analysis done above. Indeed, we can set a i ( q, ω ) := q i × h i ( q, ω ) andconclude that a mapping f from ( T S ) N to ( se (3)) N such that ψ ∗ ( f ( q, ω )) | ( q,ω ) = F | ( q,ω ) is the following one f ( q, ω ) = ω q × h ......ω N q N × h N ∈ se (3) N (cid:39) R N . We will not go into the Hamiltonian formulation of this problem; however, weremark that a similar approach works even in that situation. Indeed, followingthe derivation presented in [27], we see that for a mechanical system on ( S ) N the conjugate momentum writes T ∗ q S × ...T ∗ q N S (cid:51) π = ( π , ..., π N ) , where π i = − ˆ q i ∂L∂ω i and its components are still orthogonal to the respective base points q i ∈ S .Moreover, Hamilton’s equations take the form˙ q i = ∂H ( q, π ) ∂π i × q i , ˙ π i = ∂H ( q, π ) ∂q i × q i + ∂H ( q, π ) ∂π i × π i , which implies that setting f ( q, π ) = (cid:2) ∂ q H ( q, π ) , ∂ π H ( q, π ) , . . . , ∂ q N H ( q, π ) , ∂ π N H ( q, π ) (cid:3) we can represent even the Hamiltonian vector field of the N − fold 3D pendulumin terms of this group action. N = 2We have seen how it is possible to turn the equations of motion of a N − chain ofpendulums into the Lie group integrators setting. Now we focus on the examplewith N = 2 pendulums. The equations of motion write˙ q = ˆ ω q , ˙ q = ˆ ω q , ( q ) (cid:20) ˙ ω ˙ ω (cid:21) = (cid:20) ( − m L L | ω | ˆ q + ( m + m ) gL ˆ e ) q ( − m L L | ω | ˆ q + m gL ˆ e ) q (cid:21) , (34)where R ( q ) = (cid:20) ( m + m ) L I m L L ˆ q T ˆ q m L L ˆ q T ˆ q m L I (cid:21) . As presented above, the matrix R ( q ) defines a linear invertible map of the space T q S × T q S onto itself: A ( q ,q ) : T q S × T q S → T q S × T q S , [ ω , ω ] T → R ( q )[ ω , ω ] T . We can easily see that it is well defined since R ( q ) (cid:20) ω ω (cid:21) = (cid:20) ( m + m ) L I m L L ˆ q T ˆ q m L L ˆ q T ˆ q m L I (cid:21) (cid:20) ˆ v q ˆ v q (cid:21) = (cid:20) ˆ r q ˆ r q (cid:21) ∈ ( T S ) with r ( q, ω ) := ( m + m ) L v + m L L ˆ q ˆ v q ,r ( q, ω ) := m L L ˆ q ˆ v q + m L v . This map guarantees that if we rewrite the pair of equations for the angularvelocities in (34) as˙ ω = R − ( q ) (cid:20) ( − m L L | ω | ˆ q + ( m + m ) gL ˆ e ) q ( − m L L | ω | ˆ q + m gL ˆ e ) q (cid:21) = R − ( q ) b == A − q ,q ) ( b ) = (cid:20) h h (cid:21) ∈ T q S × T q S , then we are assured that there exist a pair of functions a , a : T S × T S → R such that ˙ ω = (cid:20) a ( q, ω ) × q a ( q, ω ) × q (cid:21) = (cid:20) h ( q ) h ( q ) (cid:21) . Since we want a i × q i = h i , we just impose a i = q i × h i and hence the wholevector field can be rewritten as ˙ q ˙ ω ˙ q ˙ ω = ω × q ( q × h ) × q ω × q ( q × h ) × q = F | ( q,ω ) , with h i = h i ( q, ω ) and (cid:20) h ( q, ω ) h ( q, ω ) (cid:21) = R − ( q ) (cid:20) ( − m L L | ω | ˆ q + ( m + m ) gL ˆ e ) q ( − m L L | ω | ˆ q + m gL ˆ e ) q (cid:21) . Therefore, we can express the whole vector field in terms of the infinitesimalgenerator of the action of SE (3) × SE (3) as ψ ∗ ( f ( q, ω )) | ( q,ω ) = F | ( q,ω ) through the function f : T S × T S → se (3) × se (3) (cid:39) R , ( q, ω ) → ( ω , q × h , ω , q × h ) . .3 Numerical experiments In this section, we present some numerical experiment for the N − chain of pen-dulums. We start by comparing the various Lie group integrators that we havetested (with the choice N = 2), and conclude by analyzing an implementationof variable step size. Lie group integrators allow to keep the evolution of thesolution in the correct manifold, which is T S × T S when N = 2. Hence, webriefly report two sets of numerical experiments. In the first one, we show theconvergence rate of all the Lie group integrators tested on this model. In thesecond one, we check how they behave in terms of preserving the two followingrelations: • q i ( t ) T q i ( t ) = 1 , i.e. q i ( t ) ∈ S , i = 1 , , • q i ( t ) T ω i ( t ) = 0 , i.e. ω i ( t ) ∈ T q i ( t ) S , i = 1 , , completing the analysis with a comparison with the classical Runge–Kutta 4and with ODE45 of MATLAB. The Lie group integrators used to obtain thefollowing experiments are Lie Euler, Lie Euler Heun, three versions of Runge–Kutta–Munthe–Kaas methods of order four and one of order three. The RKMK4with two commutators mentioned in the plots, is the one presented in Section2, while the other schemes can be found for example in [6].Figure 2 presents the plots of the errors, in logarithmic scale, obtained con-sidering as a reference solution the one given by the ODE45 method, with stricttolerance. Here, we used an exact expression for the dexp − σ function. How-ever, we could obtain the same results with a truncated version of this function,keeping a sufficiently high number of commutators. The schemes show the rightconvergence rates, so we can move to the analysis of the time evolution on T S × T S .In Figure 3 we can see the comparison of the time evolution of the 2 − norms of q ( t ) and q ( t ), for 0 ≤ t ≤ T = 5. As highlighted above, unlike classical numer-ical integrators like the one implemented in ODE45 or the Runge–Kutta 4, theLie group methods preserve the norm of the base components of the solutions,i.e. | q ( t ) | = | q ( t ) | = 1 ∀ t ∈ [0 , T ]. Therefore, as expected, these integratorspreserve the configuration manifold. However, to complete this analysis, weshow the plots making a similar comparison but with the tangentiality condi-tions. Indeed, in Figure 4 we compare the time evolutions of the inner products q ( t ) T ω ( t ) and q ( t ) T ω ( t ) for t ∈ [0 , T S × T S . As we can see, while for Liegroup methods these inner products are of the order of 10 − and 10 − , theones obtained with classical integrators show that the tangentiality conditionsare not preserved with the same accuracy.We now move to some experiments on variable stepsize. In this last part wefocus on the RKMK pair coming from Dormand–Prince method (DOPRI 5(4)[13]), which we denote with RKMK(5,4). The aim of the plots we show is tocompare the same schemes, both with constant and variable stepsize. We startby setting a tolerance and solving the system with the RKMK(5,4) scheme.22 igure 2: Convergence rate of the implemented Lie group integrators, based on global errorconsidering as a reference solution the one of ODE45, with strict tolerance.Figure 3: Visualization of the quantity 1 − q i ( t ) T q i ( t ), i = 1 ,
2, for time t ∈ [0 , S . igure 4: Visualization of the inner product q i ( t ) T ω i ( t ), i = 1 ,
2, for t ∈ [0 , T q i ( t ) S . Using the same number of time steps, we solve it again with RKMK of order 5.These experiments show that, for some tolerance and some initial conditions,the step size’s adaptivity improves the numerical approximation accuracy. Sincewe do not have an available analytical solution to quantify these two schemes’accuracy, we compare them with the solution obtained with a strict toleranceand ODE45.
Figure 5: Comparison of accuracy at final time (on the left) and step adaptation for the case N = 20 (on the right). In Figure 5, we compare the performance of the constant and variable step-size methods, where the structure of the initial condition is always the same,but what changes is the number of connected pendulums. The considered ini-24ial condition is ( q i , ω i ) = (cid:0) √ / , , √ / , , , (cid:1) , ∀ i = 1 , ..., N , and all themasses and lengths are set to 1. From these experiments we can notice situa-tions where the variable step size beats the constant one in terms of accuracy atthe final time, like the case N = 2 which we discuss in more detail afterwards.Usually these are situations with time intervals where the step can be enlargedconsiderably and others where it has to be restricted significantly. On the otherhand, as we can see in the plot on the right of Figure 5, for the case N = 20 theconstant stepsize does better than the variable one, and this seems to happenbecause there is no region where the behaviour of the solution is sufficientlysmooth. (a) ( q ( t ) , ω ( t )) (b) ( q ( t ) , ω ( t )) Figure 6: In these plots we represent the six components of the solution describing the dynam-ics of the first mass (on the left) and of the second mass (on the right), for the case N = 2. Wecompare the behaviour of the solution obtained with constant stepsize RKMK5, the variablestepsize RKMK(5,4) and ODE45. For the case N = 2, we notice a relevant improvement passing to variablestepsize. In Figures 6 and 7 we can see that, for this choice of the parameters,the solution behaves smoothly in most of the time interval, but then there is apeak in the second component of the angular velocities of both the masses, at t ≈ .
2. We can observe this behaviour both in the plots of Figure 6, where weproject the solution on the twelve components and even in Figure 7c. In thelatter, we plot two of the vector field components, i.e. the second componentsof the angular accelerations ˙ ω i ( t ), i = 1 ,
2. They show an abrupt change inthe vector field in correspondence to t ≈ .
2, where the step is considerablyrestricted. To conclude, we notice that the variable stepsize method betterfollows the reference solution in correspondence with this quite abrupt changein the solution. Figure 7b, where we focus on the last time instants, highlightsthis behaviour.At the end, we can summarize that there are situations in which variablestepsize outperforms the constant one. In these cases, the solution obtainedwith variable stepsize follows better the reference solution and adapts the step25easonably in terms of the solution’s behaviour. (a) Step adaptation (b) Zoom at final times (c) Values of ˙ ω (2) i ( t ) Figure 7: On the left, we compare the adaptation of the stepsize of RKMK(5,4) with the one ofODE45 and with the constant stepsize of RKMK5. In the center we plot the second componentof the angular velocities ω (2) i , i = 1 ,
2, and we zoom in the last time interval t ∈ [2 . ,
3] to seethat the variable stepsize version of the method better reproduces the reference solution. Onthe right, we visualize the speed of variation of second component of the angular velocities.
In this section we consider a multibody system made of two cooperatingquadrotor unmanned aerial vehicles (UAV) connected to a point mass (sus-pended load) via rigid links. This model is described in [27, 49].We consider an inertial frame whose third axis goes in the direction of grav-ity, but opposite orientation, and we denote with y ∈ R the mass point andwith y , y ∈ R the two quadrotors. We assume that the links between the twoquadrotors and the mass point are of a fixed length L , L ∈ R + . The configu-ration variables of the system are: the position of the mass point in the inertialframe, y ∈ R , the attitude matrices of the two quadrotors, ( R , R ) ∈ ( SO (3)) and the directions of the links which connect the center of mass of each quadrotorrespectively with the mass point, ( q , q ) ∈ ( S ) . The configuration manifoldof the system is Q = R × ( SO (3)) × ( S ) . In order to present the equationsof motion of the system we start by identifying T SO (3) (cid:39) SO (3) × so (3) vialeft-trivialization. This choice allows us to write the kinematic equations of thesystem as ˙ R i = R i ˆΩ i , ˙ q i = ˆ ω i q i i = 1 , , (35)where Ω , Ω ∈ R represent the angular velocities of each quadrotor, respec-tively, and ω , ω express the time derivatives of the orientations q , q ∈ S ,respectively, in terms of angular velocities, expressed with respect to the body-fixed frames. From these equations we define the trivialized Lagrangian L ( y, ˙ y, R , Ω , R , Ω , q , ω , q , ω ) : R × ( SO (3) × so (3)) × ( T S ) → R , as the difference of the total kinetic energy of the system and the total potential(gravitational) energy, L = T − U , with:26 igure 8: Two quadrotors connected to the mass point m y via massless links of lengths L i . T = 12 m y (cid:107) ˙ y (cid:107) + 12 (cid:88) i =1 ( m i (cid:107) ˙ y − L i ˆ ω i q i (cid:107) + Ω Ti J i Ω i ) , and U = − m y ge T y − (cid:88) i =1 m i ge T ( y − L i q i ) , where J , J ∈ R × are the inertia matrices of the two quadrotors and m , m ∈ R + are their respective total masses. In this system each of the two quadrotorsgenerates a thrust force, which we denote with u i = − T i R i e ∈ R , where T i is the magnitude, while e is the direction of this vector in the i − th body-fixedframe, i = 1 ,
2. The presence of these forces make it a non conservative system.Moreover, the rotors of the two quadrotors generate a moment vector, and wedenote with M , M ∈ R the cumulative moment vector of each of the twoquadrotors. To derive the Euler–Lagrange equations, a possible approach isthrough Lagrange–d’Alambert’s principle, as presented in [27]. We write themin matrix form as A ( z ) ˙ z = h ( z ) (36)where z = [ y, v, Ω , Ω , ω , ω ] T ∈ R ,A ( z ) = I M q J J − L ˆ q I − L ˆ q I , ( z ) = h ( z ) h ( z ) h ( z ) h ( z ) h ( z ) h ( z ) = v − (cid:80) i =1 m i L i (cid:107) ω i (cid:107) q i + M q ge + (cid:80) i =1 u (cid:107) i − Ω × J Ω + M − Ω × J Ω + M − L g ˆ q e − m L q × u ⊥ − L g ˆ q e − m L q × u ⊥ , where M q = m y I + (cid:80) i =1 m i q i q Ti , and u (cid:107) i , u ⊥ i are respectively the orthogonalprojection of u i along q i and to the plane T q i S , i = 1 ,
2, i.e. u (cid:107) i = q i q Ti u i , u ⊥ i = ( I − q i q Ti ) u i . These equations, coupled with the kinematic equations in(35), describe the dynamics of a point P = [ y, v, R , Ω , R , Ω , q , ω , q , ω ] ∈ M = T Q.
Since the matrix A ( z ) is invertible, we pass to the following set of equations˙ z = A − ( z ) h ( z ) := ˜ h ( z ) := ¯ h ( P ) = [¯ h ( P ) , ..., ¯ h ( P )] T . (37) We identify the phase space M with M (cid:39) T R × ( T SO (3)) × ( T S ) . Thegroup we consider is ¯ G = R × ( T SO (3)) × ( SE (3)) , where the groups are combined with a direct-product structure and R is theadditive group. For a group element g = (( a , a ) , (( B , b ) , ( B , b )) , (( C , c ) , ( C , c ))) ∈ ¯ G and a point P ∈ M in the manifold, we consider the following left action ψ g ( P ) = [ y + a , v + a , B R , Ω + b , B R , Ω + b ,C q , C ω + c × C q , C q , C ω + c × C q ] . The well-definiteness and transitivity of this action come from standard argu-ments, see for example [41]. The infinitesimal generator associated to ξ = [ ξ , ξ , η , η , η , η , µ , µ , µ , µ ] ∈ ¯ g , where ¯ g = T e ¯ G , writes ψ ∗ ( ξ ) | P = [ ξ , ξ , ˆ η R , η , ˆ η R , η , ˆ µ q , ˆ µ ω + ˆ µ q , ˆ µ q , ˆ µ ω + ˆ µ q ] . We can now focus on the construction of the function f : M → ¯ g such that ψ ∗ ( f ( P )) | P = F | P , where 28 | P = [¯ h ( P ) , ¯ h ( P ) , R ˆΩ , ¯ h ( P ) , R ˆΩ , ¯ h ( P ) , ˆ ω q , ¯ h ( P ) , ˆ ω q , ¯ h ( P )] ∈ T P M is the vector field obtained combining the equations (35) and (37). We have f ( P ) = [¯ h ( P ) , ¯ h ( P ) , R Ω , ¯ h ( P ) , R Ω , ¯ h ( P ) ,ω , q × ¯ h ( P ) , ω , q × ¯ h ( P )] ∈ ¯ g . We have obtained the local representation of the vector field F ∈ X ( M ) interms of the infinitesimal generator of the transitive group action ψ , hence wecan solve for one time step ∆ t the IVP (cid:40) ˙ σ ( t ) = dexp − σ ( t ) (cid:16) f (cid:0) ψ (exp( σ ( t )) , P ( t )) (cid:1)(cid:17) σ (0) = 0 ∈ ¯ g and then update the solution P ( t + ∆ t ) = ψ (exp( σ (∆ t )) , P ( t )).The above construction is completely independent of the control functions { u (cid:107) i , u ⊥ i , M i } i =1 , and hence it is compatible with any choice of these parameters. We tested Lie group numerical integrators for a load transportation problempresented in [49]. The control inputs { u (cid:107) i , u ⊥ i , M i } i =1 , are constructed suchthat the point mass asymptotically follows a given desired trajectory y d ∈ R ,given by a smooth function of time, and the quadrotors maintain a prescribedformation relative to the point mass. In particular, the parallel components u (cid:107) i are designed such that the payload follows the desired trajectory y d (loadtransportation problem), while the normal components u ⊥ i are designed suchthat q i converge to desired directions q id (tracking problem in S ). Finally, M i are designed to control the attitude of the quadrotors.In this experiment we focus on a simplified dynamics model, i.e. we neglectthe construction of the controllers M i for the attitude dynamics of the quadro-tors. However, the full dynamics model can also be easily integrated, once theexpressions for the attitude controllers are available.In Figure 9 we show the convergence rate of four different RKMK methods.In Figures 10-14 we show results in the tracking of a parabolic trajectory,obtained by integrating the system (36) with a RKMK method of order 4.In this paper we have considered Lie group integrators with a particular focuson problems from mechanics. In mathematical terms this means that the Liegroups and manifolds of particular interest are SO ( n ) , n = 2 , SE ( n ) , n = 2 , S and T S . Many important mechanical problems canbe modelled by using cartesian or semidirect products of these groups withoutexplicitly introducing constraints. It will be the subject of future work to exploremore examples and to aim for a more systematic approach to applying Lie groupintegrators to mechanical problems. 29 igure 9: Convergence rate of the numerical schemes compared with ODE45Figure 10: Snapshots at 0 ≤ t ≤
5. Figure 11: Components of the load position(in blue) and the desired trajectory (in red) asa function time. igure 12: Deviation of the load position fromthe target trajectory. Figure 13: Direction error of the links.Figure 14: Preservation of the norms of q , q ∈ S . eferences [1] M. Arnold and O. Br¨uls. Convergence of the generalized- α scheme for con-strained mechanical systems. Multibody Syst. Dyn. , 18(2):185–202, 2007.[2] M. Arnold, O. Br¨uls, and A. Cardona. Error analysis of generalized- α Lie group time integration methods for constrained mechanical systems.
Numer. Math. , 129(1):149–179, 2015.[3] G. Bogfjellmo and H. Marthinsen. High-order symplectic partitioned liegroup methods.
Foundations of Computational Mathematics , pages 1–38,2015.[4] F. Casas and B. Owren. Cost Efficient Lie Group Integrators in the RKMKClass.
BIT Numerical Mathematics , 43(4):723–742, 2003.[5] E. Celledoni, A. Marthinsen, and B. Owren. Commutator-free Lie groupmethods.
Future Generation Computer Systems , 19:341–352, 2003.[6] E. Celledoni, H. Marthinsen, and B. Owren. An introduction to Lie groupintegrators—basics, new developments and applications.
J. Comput. Phys. ,257(part B):1040–1061, 2014.[7] E. Celledoni and B. Owren. Lie group methods for rigid body dynamicsand time integration on manifolds.
Comput. Methods Appl. Mech. Engrg. ,192(3-4):421–438, 2003.[8] J. ´Cesi´c, V. Jukov, I. Petrovic, and D. Kuli´c. Full body human motionestimation on lie groups using 3D marker position measurements.
In Pro-ceedings of the IEEE-RAS International Conference on HumanoidRobotics,Cancun, Mexico, 15–17 November 2016 , 2016.[9] S.H. Christiansen, H.Z. Munthe-Kaas, and B. Owren. Topics in structure-preserving discretization.
Acta Numerica , 20:1–119, 2011.[10] P. E. Crouch and R. Grossman. Numerical integration of ordinary differ-ential equations on manifolds.
J. Nonlinear Sci. , 3:1–33, 1993.[11] C. Curry and B. Owren. Variable step size commutator free Lie groupintegrators.
Numer. Algorithms , 82(4):1359–1376, 2019.[12] F. Diele, L. Lopez, and R. Peluso. The Cayley transform in the numericalsolution of unitary differential systems.
Adv. Comput. Math. , 8(4):317–334,1998.[13] J. R Dormand and P. J Prince. A family of embedded runge-kutta formulae.
Journal of computational and applied mathematics , 6(1):19–26, 1980.[14] K. Engø and S. Faltinsen. Numerical integration of Lie–Poisson systemswhile preserving coadjoint orbits and energy.
SIAM J. Numer. Anal. ,39(1):128–145, 2001. 3215] H. Goldstein, C.P. Poole, and J. Safko.
Classical Mechanics . Pearson, 2013.[16] V. Guillemin and S. Sternberg. he moment map and collective motion.
Ann. Physics , 127:220–253, 1980.[17] E. Hairer, Ch. Lubich, and G. Wanner.
Geometric numerical integration ,volume 31 of
Springer Series in Computational Mathematics . Springer,Heidelberg, 2010. Structure-preserving algorithms for ordinary differentialequations, Reprint of the second (2006) edition.[18] E. Hairer, S. P. Nørsett, and G. Wanner.
Solving Ordinary DifferentialEquations I, Nonstiff Problems . Springer-Verlag, Second revised edition,1993.[19] J. Hall and M. Leok. Lie group spectral variational integrators.
Found.Comput. Math. , 17(1):199–257, 2017.[20] F. Hausdorff. Die symbolische Exponential Formel in der Gruppen theorie.
Leipziger Ber. , 58:19–48, 1906.[21] H. M. Hilber, T. J. R. Hughes, and R. L. Taylor. Improved numerical dissi-pation for time integration algorithms in structural dynamics.
EarthquakeEngineering & Structural Dynamics , 5(3):283–292, 1977. cited By 1683.[22] D. Holm.
Geometric Mechanics: Part II: Rotating, Translating and Rolling .World Scientific Publishing Company, 2008.[23] D. Holm, J. Marsden, and T. Ratiu. The Euler–poincar´e equations andsemidirect products with applications to continuum theories.
Adv. in Math. ,137:1–81, 1998.[24] S. Holzinger and J. Gerstmayr. Time integration of rigid bodies modelledwith three rotation parameters.
Multibody Sys Dyn , 2021.[25] A. Iserles, H. Z. Munthe-Kaas, S. P. Nørsett, and A. Zanna. Lie-groupmethods.
Acta Numerica , 9:215–365, 2000.[26] T. Lee, M. Leok, and N. H. McClamroch. Lie group variational integratorsfor the full body problem.
Comput. Methods Appl. Mech. Engrg. , 196(29-30):2907–2924, 2007.[27] T. Lee, M. Leok, and N. H. McClamroch.
Global formulations of Lagrangianand Hamiltonian dynamics on manifolds . Interaction of Mechanics andMathematics. Springer, Cham, 2018. A geometric approach to modelingand analysis.[28] T. Leitz and S. Leyendecker. Galerkin Lie-group variational integratorsbased on unit quaternion interpolation.
Comput. Methods Appl. Mech.Engrg. , 338:333–361, 2018. 3329] N. E. Leonard and J. E. Marsden. Stability and drift of underwater vehicledynamics: Mechanical systems with rigid motion symmetry.
Physica D ,105(1–3):130–162, 1997.[30] D. Lewis and J. C. Simo. Conserving algorithms for the dynamics of Hamil-tonian systems of Lie groups.
J. Nonlinear Sci. , 4:253–299, 1994.[31] A. Lundervold and H. Z. Munthe-Kaas. On algebraic structures of nu-merical integration on vector spaces and manifolds. In
Fa`a di Bruno Hopfalgebras, Dyson-Schwinger equations, and Lie-Butcher series , volume 21 of
IRMA Lect. Math. Theor. Phys. , pages 219–263. Eur. Math. Soc., Z¨urich,2015.[32] J. E. Marsden, T. Ratiu, and A. Weinstein. Semi-direct products and re-duction in mechanics.
Transactions of the American Mathematical Society ,281:147–77, 1884.[33] J. E. Marsden and T. S. Ratiu.
Introduction to Mechanics and Symmetry .Springer-Verlag, 1994.[34] J. E. Marsden and T. S. Ratiu.
Introduction to Mechanics and Symme-try . Number 17 in Texts in Applied Mathematics. Springer-Verlag, secondedition, 1999.[35] R. H. Merson. An operational method for the study of integration processes.In
Proc. Symp. Data Processing , 1957.[36] J. Moser and A. P. Veselov. Discrete versions of some classical integrablesystems and factorization of matrix polynomials.
Comm. Math. Phys. ,139(2):217–243, 1991.[37] A. M¨uller. Coordinate mappings for rigid body motions.
ASME Journalof Computational and Nonlinear Dynamics , 12, 2017.[38] H. Munthe-Kaas. Runge–Kutta methods on Lie groups.
BIT , 38(1):92–111,1998.[39] H. Munthe-Kaas. High order Runge–Kutta methods on manifolds.
Appl.Numer. Math. , 29:115–127, 1999.[40] H. Munthe-Kaas and B. Owren. Computations in a free Lie algebra.
Phil.Trans. Royal Soc. A , 357:957–981, 1999.[41] P. J. Olver.
Applications of Lie groups to differential equations , volume107. Springer Science & Business Media, 2000.[42] B. Owren. Order conditions for commutator-free Lie group methods.
J.Phys. A , 39(19):5585–5599, 2006. 3443] B. Owren. Lie group integrators. In
Discrete mechanics, geometric inte-gration and Lie-Butcher series , volume 267 of
Springer Proc. Math. Stat. ,pages 29–69. Springer, Cham, 2018.[44] B. Owren and A. Marthinsen. Integration methods based on canonicalcoordinates of the second kind.
Numer. Math. , 87(4):763–790, 2001.[45] J. Park and W. Chung. Geometric integration on Euclidean group withapplication to articulated multibody systems.
J. CAM , 21, 2005.[46] T. Ratiu. Euler-poisson equations on lie algebras and the n-dimensionalheavy rigid body.
Proc. Nat. Acad. Sci. USA , pages 1327–1328, 1981.[47] A. Saccon. Midpoint rule for variational integrators on Lie groups.
Internat.J. Numer. Methods Engrg. , 78(11):1345–1364, 2009.[48] J. C. Simo and L. Vu-Quoc. On the dynamics of finite-strain rods under-going large motions — a geometrically exact approach.
Comput. MethodsAppl. Mech. Engrg. , 66:125–161, 1988.[49] Lee T., K. Sreenath, and V. Kumar. Geometric control of cooperatingmultiple quadrotor uavs with a suspended payload. , pages 5510–5515, 2013.[50] A. P. Veselov. Integrable systems with discrete time, and difference opera-tors.
Funktsional. Anal. i Prilozhen. , 22(2):1–13, 96, 1988.[51] A. M. Vinogradov and B. Kupershmidt. The structure of hamiltonianmechanics.
Funktsional. Anal. i Prilozhen. , 32:177–243, 1977.[52] V. Wieloch and M. Arnold. BDF integrators for constrained mechanicalsystems on Lie groups.