Lie group integrators for mechanical systems

Abstract

Since they were introduced in the 1990s, Lie group integrators have become a method of choice in many application areas. These include multibody dynamics, shape analysis, data science, image registration and biophysical simulations. Two important classes of intrinsic Lie group integrators are the Runge–Kutta–Munthe–Kaas methods and the commutator free Lie group integrators. We give a short introduction to these classes of methods. The Hamiltonian framework is attractive for many mechanical problems, and in particular we shall consider Lie group integrators for problems on cotangent bundles of Lie groups where a number of different formulations are possible. There is a natural symplectic structure on such manifolds and through variational principles one may derive symplectic Lie group integrators. We also consider the practical aspects of the implementation of Lie group integrators, such as adaptive time stepping. The theory is illustrated by applying the methods to two nontrivial applications in mechanics. One is the N-fold spherical pendulum where we introduce the restriction of the adjoint action of the group $SE\left(3\right)$ to $T{S}^2$, the tangent bundle of the two-dimensional sphere. Finally, we show how Lie group integrators can be applied to model the controlled path of a payload being transported by two rotors. This problem is modelled on ${\mathbb{R}}^6\times {\left(SO\left(3\right)\times \mathfrak{s}\mathfrak{o}\left(3\right)\right)}^2\times (T{S}^2{)}^2$ and put in a format where Lie group integrators can be applied.

Keywords:

• Lie group integrators
• mechanical systems
• adaptive time stepping
• hamiltonian systems
• multibody systems

2010 AMS Subject Classifications:

• 65L05
• 70E55

1. Introduction

In many physical problems, including multi-body dynamics, the configuration space is not a linear space, but rather consists of a collection of rotations and translations. A simple example is the free rigid body whose configuration space consists of rotations in 3D. A more advanced example is the simplified model of the human body, where the skeleton at a given time is described as a system of interacting rods and joints. Mathematically, the structure of such problems is usually best described as a manifold. Since manifolds by definition can be equipped with local coordinates, one can always describe and simulate such systems locally as if they were linear spaces. There are of course many choices of local coordinates, for rotations some famous ones are: Euler angles, the Tait–Bryan angles commonly used in aerospace applications, the unit length quaternions and the exponentiated skew-symmetric $3\times 3$-matrices. Lie group integrators represent a somewhat different strategy. Rather than specifying a choice of local coordinates from the outset, in this approach the model and the numerical integrator are expressed entirely in terms of a Lie group and its action on the phase space. This often leads to a more abstract and simpler formulation of the mechanical system and of the numerical schemes, deferring further details to the implementation phase.

In the literature, one can find many different types and formats of Lie group integrators. Some of these are completely general and intrinsic, meaning that they only make use of inherent properties of Lie groups and manifolds as was suggested in [6,11,41]. But many numerical methods have been suggested that add structure or utilize properties which are specific to a particular Lie group or manifold. Notable examples of this are the methods based on canonical coordinates of the second kind [45], and the methods based on the Cayley transformation [13,32], applicable, e.g. to the rotation groups and Euclidean groups. In some applications, e.g. in multi-body systems, it may be useful to formulate the problem as a mix between Lie groups and kinematic constraints, introducing for instance Lagrange multipliers. Sometimes this may lead to more practical implementations where a basic general setup involving Lie groups can be further equipped with different choices of constraints depending on the particular application. Such constrained formulations are outside the scope of the present paper. It should also be noted that the Lie group integrators devised here do not make any a-priori assumptions about how the manifold is represented.

The applications of Lie group integrators for mechanical problems also have a long history, two of the early important contributions were the Newmark methods of Simo and Vu–Quoc [50] and the symplectic and energy-momentum methods by Lewis and Simo [32]. Mechanical systems are often described as Euler–Lagrange equations or as Hamiltonian systems on manifolds, with or without external forces [28]. Important ideas for the discretization of mechanical systems originated also from the work of Moser and Veselov [39,51] on discrete integrable systems. This work served as motivation for further developments in the field of geometric mechanics and for the theory of (Lie group) discrete variational integrators [20,27,30]. The majority of Lie group methods found in the literature are one-step type generalizations for classical methods, such as Runge–Kutta type formulas. In mechanical engineering, the classical BDF methods have played an important role and were recently generalized [54] to Lie groups. Similarly, the celebrated α-method for linear spaces proposed by Hilber, Hughes and Taylor [22] has been popular for solving problems in multibody dynamics, and in [1,2,4] this method is generalized to a Lie group integrator.

The literature on Lie group integrators is rich and diverse, the interested reader may consult the surveys [7,10,26,46] and Chapter 4 of the monograph [18] for further details.

In this paper, we discuss different ways of applying Lie group integrators to simulating the dynamics of mechanical multi-body systems. Our point of departure is the formulation of the models as differential equations on manifolds. Assuming to be given either a Lie group acting transitively on the manifold $\mathcal{M}$ or a set of frame vector fields on $\mathcal{M}$, we use them to describe the mechanical system and further to build the numerical integrator. We shall here mostly consider schemes of the types commonly known as Crouch–Grossman methods [11], Runge–Kutta–Munthe–Kaas methods [40,41] and Commutator-free Lie group methods [6].

The choice of Lie group action is often not unique and thus the same mechanical system can be described in different equivalent ways. Under numerical discretization, different formulations can lead to the conservation of different geometric properties of the mechanical system. In particular, we explore the effect of these different formulations on a selection of examples in multi-body dynamics. Lie group integrators have been successfully applied for the simulation of mechanical systems, and in problems of control, bio-mechanics and other engineering applications, see, e.g. [9,25,27,47]. The present work is motivated by applications in modelling and simulation of slender structures like Cosserat rods and beams [50], and one of the examples presented here is the application to a chain of pendula. Another example considers an application for the controlled dynamics of a multibody system.

In Section 2, we give a review of the methods using only the essential intrinsic tools of Lie group integrators. The algorithms are simple and amenable for a coordinate-free description suited to object oriented implementations. In Section 3, we discuss Hamiltonian systems on Lie groups, and we present three different Lie group formulations of the heavy top equations. These systems (and their Lagrangian counterpart) often arise in applications as building blocks of more realistic systems which comprise also damping and control forces. In Section 4, we discuss some ways of adapting the integration step size in time. In Section 5, we consider the application to a chain of pendula. And in Section 6, we consider the application of a multi-body system of interest in the simulation and control of drone dynamics.

2. Lie group integrators

2.1. The formulation of differential equations on manifolds

Lie group integrators solve differential equations whose solution evolve on a manifold $\mathcal{M}$. For ease of notation, we restrict the discussion to the case of autonomous vector fields, although allowing for explicit t-dependence could easily have been included. This means that we seek a curve $y\left(t\right)\in \mathcal{M}$ whose tangent at any point coincides with a vector field $F\in \mathcal{X}\left(\mathcal{M}\right)$ and passing through a designated initial value $y_0$ at $t=t_0$ (1) $\dot{y}\left(t\right)=F{|}_{y\left(t\right)},y\left(t_0\right)=y_0.$(1) Before addressing numerical methods for solving (1(1) $\dot{y}\left(t\right)=F{|}_{y\left(t\right)},y\left(t_0\right)=y_0.$(1) ) it is necessary to introduce a convenient way of representing the vector field F. There are different ways of doing this. One is to furnish $\mathcal{M}$ with a transitive action $\psi :G\times \mathcal{M}\to \mathcal{M}$ by some Lie group G of dimension $d\ge \dim \mathcal{M}$. We denote the action of g on m as $g\cdot m$, i.e. $g\cdot m=\psi (g,m)$. Let $\mathfrak{g}$ be the Lie algebra of G, and denote by $\exp :\mathfrak{g}\to G$ the exponential map. We define ${\psi }_{\ast }:\mathfrak{g}\to \mathcal{X}\left(\mathcal{M}\right)$ to be the infinitesimal generator of the action, i.e. (2) ${\left.F_{\xi }\right|}_m={\left.{\psi }_{\ast }\left(\xi \right)\right|}_m={\left.\frac{\mathrm{d}}{\mathrm{d}t}\right|}_{t=0}\psi (\exp (t\xi ),m)$(2) The transitivity of the action now ensures that ${\psi }_{\ast }\left(\mathfrak{g}\right){|}_m=T_m\mathcal{M}$ for any $m\in \mathcal{M}$, such that any tangent vector $v_m\in T_m\mathcal{M}$ can be represented as $v_m={\psi }_{\ast }\left({\xi }_v\right){|}_m$ for some ${\xi }_v\in \mathfrak{g}$ (${\xi }_v$ may not be unique). Consequently, for any vector field $F\in \mathcal{X}\left(\mathcal{M}\right)$ there exists a map $f:\mathcal{M}\to \mathfrak{g}$¹ such that (3) $F{|}_m={\left.{\psi }_{\ast }\left(f\right(m\left)\right)\right|}_m,\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}m\in \mathcal{M}$(3) This is the original tool [41] for representing a vector field on a manifold with a group action. Another approach was used in [11] where a set of frame vector fields $E_1,\dots ,E_d$ in $\mathcal{X}\left(\mathcal{M}\right)$ was introduced assuming that for every $m\in \mathcal{M}$, $\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}\{{\left.E_1\right|}_m,\dots ,{\left.E_d\right|}_m\}=T_m\mathcal{M}.$ Then, for any vector field $F\in \mathcal{X}\left(\mathcal{M}\right)$ there are, in general non-unique, functions $f_i:\mathcal{M}\to \mathbb{R}$, which can be chosen with the same regularity as F, such that $F{|}_m=\sum _{i=1}^d{f}_i\left(m\right){\left.E_i\right|}_m.$ A fixed vector $\xi \in {\mathbb{R}}^d$ will define a vector field $F_{\xi }$ on $\mathcal{M}$ similar to (2(2) ${\left.F_{\xi }\right|}_m={\left.{\psi }_{\ast }\left(\xi \right)\right|}_m={\left.\frac{\mathrm{d}}{\mathrm{d}t}\right|}_{t=0}\psi (\exp (t\xi ),m)$(2) ) (4) ${\left.F_{\xi }\right|}_m=\sum _{i=1}^d{\xi }_i{E}_i{|}_m$(4) If ${\xi }_i=f_i\left(p\right)$ for some $p\in \mathcal{M}$, the corresponding $F_{\xi }$ will be a vector field in the linear span of the frame which coincides with F at the point p. Such a vector field was named by Crouch and Grossman [11] as a the vector field frozen at p.

The two formulations just presented are in many cases connected and can then be used in an equivalent manner. Suppose that $e_1,\dots ,e_d$ is a basis of the Lie algebra $\mathfrak{g}$, then we can simply define frame vector fields as $E_i={\psi }_{\ast }\left(e_i\right)$ and the vector field we aim to describe is $F{|}_m={\left.{\psi }_{\ast }\left(f\right(m\left)\right)\right|}_m={\left.{\psi }_{\ast }\left(\sum _i{f}_i\left(m\right)e_i\right)\right|}_m=\sum _i{f}_i{\left.E_i\right|}_m.$ As mentioned above, there is a non-uniqueness issue when defining a vector field by means of a group action or a frame. A more fundamental description can be obtained using the machinery of connections. The assumption is that the simply connected manifold $\mathcal{M}$ is equipped with a connection which is flat and has constant torsion. Then $F_p$, the frozen vector field of F at p defined above, can be defined as the unique element $F_p\in \mathcal{X}\left(\mathcal{M}\right)$ satisfying

1. $F_p{|}_p=F{|}_p$

2. ${\mathrm{\nabla }}_X{F}_p=0$ for any $X\in \mathcal{X}\left(M\right)$.

So $F_p$ is the vector field that coincides with F at p and is parallel transported to any other point on $\mathcal{M}$ by the connection ∇. Since the connection is flat, the parallel transport from the point p to another point $m\in \mathcal{M}$ does not depend on the chosen path between the two points. For further details, see, e.g. Lundervold and Munthe-Kaas [33].

Example 2.1

For mechanical systems on Lie groups, two important constructions are the adjoint and coadjoint representations. For every $g\in G$, there is an automorphism ${\mathrm{A}\mathrm{d}}_g:\mathfrak{g}\to \mathfrak{g}$ defined as ${\mathrm{A}\mathrm{d}}_g\left(\xi \right)=T{L}_g\circ T{R}_{g^{-1}}\left(\xi \right)$ where $L_g$ and $R_g$ are the left and right multiplications respectively, $L_g\left(h\right)=gh$ and $R_g\left(h\right)=hg$. Since $\mathrm{A}\mathrm{d}$ is a representation, i.e. ${\mathrm{A}\mathrm{d}}_{gh}={\mathrm{A}\mathrm{d}}_g\circ {\mathrm{A}\mathrm{d}}_h$ it also defines a left Lie group action by G on $\mathfrak{g}$. From this definition and a duality pairing $〈\cdot ,\cdot 〉$ between $\mathfrak{g}$ and ${\mathfrak{g}}^{\ast }$, we can also derive a representation on ${\mathfrak{g}}^{\ast }$ denoted ${\mathrm{A}\mathrm{d}}_g^{\ast }$, simply by $〈{\mathrm{A}\mathrm{d}}_g^{\ast }\left(\mu \right),\xi 〉=〈\mu ,{\mathrm{A}\mathrm{d}}_g\left(\xi \right)〉,\xi \in \mathfrak{g},\mu \in {\mathfrak{g}}^{\ast }.$ The action $g\cdot \mu ={\mathrm{A}\mathrm{d}}_{g^{-1}}^{\ast }\left(\mu \right)$ has infinitesimal generator given as ${\left.{\psi }_{\ast }\left(\xi \right)\right|}_{\mu }=-{\mathrm{a}\mathrm{d}}_{\xi }^{\ast }\mu$ Following [35], for a Hamiltonian $H:T^{\ast }G\to \mathbb{R}$, define $H^{-}$ to be its restriction to ${\mathfrak{g}}^{\ast }$. Then the Lie–Poisson reduction of the dynamical system is defined on ${\mathfrak{g}}^{\ast }$ as $\dot{\mu }=-{\mathrm{a}\mathrm{d}}_{\frac{\mathrm{\partial }{H}^{-}}{\mathrm{\partial }\mu }}^{\ast }\mu$ and this vector field is precisely of the form (3(3) $F{|}_m={\left.{\psi }_{\ast }\left(f\right(m\left)\right)\right|}_m,\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{l}m\in \mathcal{M}$(3) ) with $f\left(\mu \right)=\frac{\mathrm{\partial }{H}^{-}}{\mathrm{\partial }\mu }\left(\mu \right)$. A side effect of this is that the integral curves of these Lie–Poisson systems preserve coadjoint orbits, making the coadjoint action an attractive choice for Lie group integrators.

Let us now detail the situation for the very simple case where $G=SO\left(3\right)$. The Lie algebra $\mathfrak{s}\mathfrak{o}\left(3\right)$ can be modelled as $3\times 3$ skew-symmetric matrices, and via the standard basis we identify each such matrix $\hat{\xi }$ by a vector $\xi \in {\mathbb{R}}^3$, this identification is known as the hat map (5) $\hat{\xi }=\left[\begin{array}{ccc}0& -{\xi }_3& {\xi }_2\\ {\xi }_3& 0& -{\xi }_1\\ -{\xi }_2& {\xi }_1& 0\end{array}\right]$(5) Now, we also write the elements of $\mathfrak{s}\mathfrak{o}(3{)}^{\ast }$ as vectors in ${\mathbb{R}}^3$ with duality pairing $〈\mu ,\xi 〉={\mu }^T\xi$. With these representations, we find that the coadjoint action can be expressed as $g\cdot \mu =\psi (g,\mu )={\mathrm{A}\mathrm{d}}_{g^{-1}}^{\ast }\mu =g\mu$ the rightmost expression being a simple matrix-vector multiplication. Since g is orthogonal, it follows that the coadjoint orbits foliate 3-space into spherical shells, and the coadjoint action is transitive on each of these orbits. The free rigid body can be cast as a problem on $TSO(3{)}^{\ast }$ with a left invariant Hamiltonian which reduces to the function $H^{-}\left(\mu \right)=\frac{1}{2}〈\mu ,{\mathbb{I}}^{-1}\mu 〉$ on $\mathfrak{s}\mathfrak{o}(3{)}^{\ast }$ where $\mathbb{I}:\mathfrak{s}\mathfrak{o}\left(3\right)\to \mathfrak{s}\mathfrak{o}(3{)}^{\ast }$ is the inertia tensor. From this, we can now set $f\left(\mu \right)=\mathrm{\partial }{H}^{-}/\mathrm{\partial }\mu ={\mathbb{I}}^{-1}\mu$. We then recover the Euler free rigid body equation as $\dot{\mu }={\left.{\psi }_{\ast }\left(f\right(\mu )\right|}_{\mu }=-{\mathrm{a}\mathrm{d}}_{{\mathbb{I}}^{-1}\mu }^{\ast }\mu =-{\mathbb{I}}^{-1}\mu \times \mu$ where the last expression involves the cross product of vectors in ${\mathbb{R}}^3$.

2.2. Two classes of Lie group integrators

The simplest numerical integrator for linear spaces is the explicit Euler method. Given an initial value problem $\dot{y}=F\left(y\right)$, $y\left(0\right)=y_0$ the method is defined as $y_{n+1}=y_n+hF\left(y_n\right)$ 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}\left(y\right)=F\left(y_n\right)$, that is $y_{n+1}=\exp \left(h{F}_{y_n}\right)y_n$ This definition of the Euler method makes sense also when F is replaced by a vector field on some manifold. In this general situation, it is known as the Lie–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 from the manifold $\mathcal{M}$ to the Lie algebra $\mathfrak{g}$, take a time step, and map the result back to $\mathcal{M}$. The transformation we use is $y\left(t\right)=\exp \left(\sigma \right(t\left)\right)\cdot y_0,\sigma \left(0\right)=0.$ The transformed differential equation for $\sigma \left(t\right)$ makes use of the derivative of the exponential mapping, the reader should consult [41] for details about the derivation, we give the final result (6) $\dot{\sigma }\left(t\right)={\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_{\sigma \left(t\right)}^{-1}\left(f\right(\exp \left(\sigma \right(t\left)\right)\cdot y_0\left)\right)$(6) The map $v\mapsto {\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_u\left(v\right)$ is linear and invertible when u belongs to some sufficiently small neighbourhood of $0\in \mathfrak{g}$. It has an expansion in nested Lie brackets [21]. Using the operator ${\mathrm{a}\mathrm{d}}_u\left(v\right)=[u,v]$ and its powers ${\mathrm{a}\mathrm{d}}_u^{2}v=[u,[u,v\left]\right]$ , etc., one can write (7) ${\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_u\left(v\right)={\left.\frac{e^z-1}{z}\right|}_{z={\mathrm{a}\mathrm{d}}_u}\left(v\right)=v+\frac{1}{2}[u,v]+\frac{1}{6}[u,[u,v\left]\right]+\cdots$(7) and the inverse is (8) ${\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_u^{-1}\left(v\right)={\left.\frac{z}{e^z-1}\right|}_{z={\mathrm{a}\mathrm{d}}_u}\left(v\right)=v-\frac{1}{2}[u,v]+\frac{1}{12}[u,[u,v\left]\right]+\cdots$(8) The RKMK methods are now obtained simply by applying some standard Runge–Kutta method to the transformed Equation (6(6) $\dot{\sigma }\left(t\right)={\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_{\sigma \left(t\right)}^{-1}\left(f\right(\exp \left(\sigma \right(t\left)\right)\cdot y_0\left)\right)$(6) ) with a time step h, using initial value $\sigma \left(0\right)=0$. This leads to an output ${\sigma }_1\in \mathfrak{g}$ and one simply sets $y_1=\exp \left({\sigma }_1\right)\cdot y_0$. Then one repeats the procedure replacing $y_0$ by $y_1$ in the next step, etc. While solving (6(6) $\dot{\sigma }\left(t\right)={\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_{\sigma \left(t\right)}^{-1}\left(f\right(\exp \left(\sigma \right(t\left)\right)\cdot y_0\left)\right)$(6) ), one needs to evaluate ${\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_u^{-1}\left(v\right)$ as a part of the process. This can be done by truncating the series (8(8) ${\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_u^{-1}\left(v\right)={\left.\frac{z}{e^z-1}\right|}_{z={\mathrm{a}\mathrm{d}}_u}\left(v\right)=v-\frac{1}{2}[u,v]+\frac{1}{12}[u,[u,v\left]\right]+\cdots$(8) ) since $\sigma \left(0\right)=0$ implies that we always evaluate ${\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_u^{-1}$ with $u=\mathcal{O}\left(h\right)$, and thus, the kth iterated commutator ${\mathrm{a}\mathrm{d}}_u^k=\mathcal{O}\left(h^k\right)$. For a given Runge–Kutta method, there are some clever tricks that can be done to minimize the total number of commutators to be included from the expansion of ${\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_u^{-1}v$ [5,42]. We give here one concrete example of an RKMK method proposed in [5] $\begin{array}{rl}{f}_{n,1}& =hf\left(y_n\right),\\ f_{n,2}& =hf(\exp (\frac{1}{2}{f}_{n,1})\cdot y_n),\\ f_{n,3}& =hf(\exp (\frac{1}{2}{f}_{n,2}-\frac{1}{8}[f_{n,1},f_{n,2}])\cdot y_n),\\ f_{n,4}& =hf(\exp (f_{n,3})\cdot y_n),\\ y_{n+1}& =\exp \left(\frac{1}{6}\right(f_{n,1}+2f_{n,2}+2f_{n,3}+f_{n,4}-\frac{1}{2}[f_{n,1},f_{n,4}]\left)\right)\cdot y_n.\end{array}$ The other option is to compute the exact expression for ${\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_u^{-1}\left(v\right)$ for the particular Lie algebra we use. For instance, it was shown in [8] that for the Lie algebra $\mathfrak{s}\mathfrak{o}\left(3\right)$ one has ${\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_u^{-1}\left(v\right)=v-\frac{1}{2}u\times v+{\alpha }^{-2}\left(1-\frac{\alpha }{2}\cot \frac{\alpha }{2}\right)u\times (u\times v)$ We will present the corresponding formula for $\mathfrak{s}\mathfrak{e}\left(3\right)$ in Section 2.3.

The second class of Lie group integrators to be considered here are the commutator-free methods, named this way in [6] to emphasize the contrast to RKMK schemes which usually include commutators in the method format. These schemes include the Crouch–Grossman methods [11] and they have the format $\begin{array}{rl}{Y}_{n,r}& =\exp \left(h\sum _k{\alpha }_{r,J}^k{f}_{n,k}\right)\cdots \exp \left(h\sum _k{\alpha }_{r,1}^k{f}_{n,k}\right)\cdot y_n\\ f_{n,r}& =f\left(Y_{n,r}\right)\\ y_{n+1}& =\exp \left(h\sum _k{\text{β}}_J^k{f}_{n,k}\right)\cdots \exp \left(h\sum _k{\text{β}}_1^k{f}_{n,k}\right)\cdot y_n\end{array}$ Here the Runge–Kutta coefficients ${\alpha }_{r,j}^k$, ${\text{β}}_j^r$ are related to a classical Runge–Kutta scheme with coefficients $a_r^k$, $b_r$ in that $a_r^k=\sum _j{\alpha }_{r,j}^k$ and $b_r=\sum _j{\text{β}}_j^r$. The ${\alpha }_{r,j}^k$, ${\text{β}}_j^r$ are usually chosen to obtain computationally inexpensive schemes with the highest possible order of convergence. The computational complexity of the above schemes depends on the cost of computing an exponential as well as of evaluating the vector field. Therefore it makes sense to keep the number of exponentials J in each stage as low as possible, and possibly also the number of stages s. A trick proposed in [6] was to select coefficients that make it possible to reuse exponentials from one stage to another. This is perhaps best illustrated through the following example from [6], a generalization of the classical fourth-order Runge–Kutta method. (9) $\begin{array}{r}\begin{array}{rl}{Y}_{n,1}& =y_n\\ Y_{n,2}& =\exp \left(\frac{1}{2}h{f}_{n,1}\right)\cdot y_n\\ Y_{n,3}& =\exp \left(\frac{1}{2}h{f}_{n,2}\right)\cdot y_n\\ Y_{n,4}& =\exp (h{f}_{n,3}-\frac{1}{2}h{f}_{n,1})\cdot Y_{n,2}\\ y_{n+\frac{1}{2}}& =\exp \left(\frac{1}{12}h\right(3f_{n,1}+2f_{n,2}+2f_{n,3}-f_{n,4}\left)\right)\cdot y_n\\ y_{n+1}& =\exp \left(\frac{1}{12}h\right(-f_{n,1}+2f_{n,2}+2f_{n,3}+3f_{n,4}\left)\right)\cdot y_{n+\frac{1}{2}}\end{array}\end{array}$(9) where $f_{n,i}=f\left(Y_{n,i}\right)$. Here, we see that one exponential is saved in computing $Y_{n,4}$ by making use of $Y_{n,2}$.

2.3. An exact expression for ${\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_u^{-1}\left(v\right)$ in $\mathfrak{s}\mathfrak{e}\left(3\right)$

As an alternative to using a truncated version of the infinite series for ${\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_u^{-1}$ (8(8) ${\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_u^{-1}\left(v\right)={\left.\frac{z}{e^z-1}\right|}_{z={\mathrm{a}\mathrm{d}}_u}\left(v\right)=v-\frac{1}{2}[u,v]+\frac{1}{12}[u,[u,v\left]\right]+\cdots$(8) ), one can consider exact expressions obtained for certain Lie algebras. Since $\mathfrak{s}\mathfrak{e}\left(3\right)$ is particularly important in applications to mechanics, we give here its exact expression. For this, we represent elements of $\mathfrak{s}\mathfrak{e}\left(3\right)$ as a pair $(A,a)\in {\mathbb{R}}^3\times {\mathbb{R}}^3\cong {\mathbb{R}}^6$, the first component corresponding to a skew-symmetric matrix $\hat{A}$ via (5(5) $\hat{\xi }=\left[\begin{array}{ccc}0& -{\xi }_3& {\xi }_2\\ {\xi }_3& 0& -{\xi }_1\\ -{\xi }_2& {\xi }_1& 0\end{array}\right]$(5) ) and a is the translational part. Now, let $\phi \left(z\right)$ be a real analytic function at z = 0. We define ${\phi }_{+}\left(z\right)=\frac{\phi \left(iz\right)+\phi (-iz)}{2},{\phi }_{-}\left(z\right)=\frac{\phi \left(iz\right)-\phi (-iz)}{2i}$ We next define the four functions $g_1\left(z\right)=\frac{{\phi }_{-}\left(z\right)}{z},{\tilde{g}}_1\left(z\right)=\frac{g_1^{\prime }\left(z\right)}{z},g_2\left(z\right)=\frac{\phi \left(0\right)-{\phi }_{+}\left(z\right)}{z^2},{\tilde{g}}_2\left(z\right)=\frac{g_2^{\prime }\left(z\right)}{z}$ and the two scalars $\rho =A^{T}a$, $\alpha =\parallel A{\parallel }_2$. One can show that for any $(A,a)$ and $(B,b)$ in $\mathfrak{s}\mathfrak{e}\left(3\right)$, it holds that $\phi \left({\mathrm{a}\mathrm{d}}_{(A,a)}\right)(B,b)=(C,c)$ where $\begin{array}{rl}C& =\phi \left(0\right)B+g_1\left(\alpha \right)A\times B+g_2\left(\alpha \right)A\times (A\times B)\\ c& =\phi \left(0\right)b+g_1\left(\alpha \right)(a\times B+A\times b)+\rho {\tilde{g}}_1\left(\alpha \right)A\times B+\rho {\tilde{g}}_2\left(\alpha \right)A\times (A\times B)\\ & +g_2\left(\alpha \right)(a\times (A\times B)+A\times (a\times B)+A\times (A\times b\left)\right)\end{array}$ Considering for instance (8(8) ${\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_u^{-1}\left(v\right)={\left.\frac{z}{e^z-1}\right|}_{z={\mathrm{a}\mathrm{d}}_u}\left(v\right)=v-\frac{1}{2}[u,v]+\frac{1}{12}[u,[u,v\left]\right]+\cdots$(8) ), we may now use $\phi \left(z\right)=\frac{z}{e^z-1}$ to calculate $g_1\left(z\right)=-\frac{1}{2},{\tilde{g}}_1\left(z\right)=0,g_2\left(z\right)=\frac{1-\frac{z}{2}\cot \frac{z}{2}}{z^2},{\tilde{g}}_2\left(z\right)=\frac{1}{z}\frac{\mathrm{d}}{\mathrm{d}z}{g}_2\left(z\right),\phi \left(0\right)=1.$ and thereby obtain an expression for ${\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_{(A,a)}^{-1}(B,b)$ with the formula above.

Similar types of formulas are known for computing the matrix exponential as well as functions of the $\mathrm{a}\mathrm{d}$-operator for several other Lie groups of small and medium dimension. For instance in [34], a variety of coordinate mappings for rigid body motions are discussed. For Lie algebras of larger dimension, both the exponential mapping and ${\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_u^{-1}$ may become computationally infeasible. For these cases, one may benefit from replacing the exponential by some other coordinate map for the Lie group $\varphi :\mathfrak{g}\to G$. One option is to use canonical coordinates of the second kind [45]. Then for some Lie groups such as the orthogonal, unitary and symplectic groups, there exist other maps that can be used and which are computationally less expensive. A popular choice is the Cayley transformation [13].

3. Hamiltonian systems on Lie groups

In this section, we consider Hamiltonian systems on Lie groups. These systems (and their Lagrangian counterpart) often appear in mechanics applications as building blocks for more realistic systems with additional damping and control forces. We consider canonical systems on the cotangent bundle of a Lie group and Lie–Poisson systems which can arise by symmetry reduction or otherwise. We illustrate various cases with different formulations of the heavy top system.

3.1. Semi-direct products

The coadjoint action by G on ${\mathfrak{g}}^{\ast }$ is denoted ${\mathrm{A}\mathrm{d}}_g^{\ast }$ defined for any $g\in G$ as (10) $〈{\mathrm{A}\mathrm{d}}_g^{\ast }\mu ,\xi 〉=〈\mu ,{\mathrm{A}\mathrm{d}}_g\xi 〉,\mathrm{\forall }\xi \in \mathfrak{g},$(10) where $\mathrm{A}\mathrm{d}:\mathfrak{g}\to \mathfrak{g}$ is the adjoint representation and for a duality pairing $〈\cdot ,\cdot 〉$ between ${\mathfrak{g}}^{\ast }$ and $\mathfrak{g}$.

We consider the cotangent bundle of a Lie group G, $T^{\ast }G$ and identify it with $G\times {\mathfrak{g}}^{\ast }$ using the right multiplication $R_g:G\to G$ and its tangent mapping $R_{g\ast }:=T{R}_g$. The cartesian product $G\times {\mathfrak{g}}^{\ast }$ can be given a semi-direct product structure that turns it into a Lie group $\mathbf{G}:=G⋉{\mathfrak{g}}^{\ast }$ where the group multiplication is (11) $(g_1,{\mu }_1)\cdot (g_2,{\mu }_2)=(g_1\cdot g_2,{\mu }_1+{\mathrm{A}\mathrm{d}}_{g_1^{-1}}^{\ast }{\mu }_2).$(11) Acting by left multiplication any vector field $F\in \mathcal{X}\left(\mathbf{G}\right)$ is expressed by means of a map $f:\mathbf{G}\to T_e\mathbf{G}$, (12) $F(g,\mu )=T_e{R}_{(g,\mu )}f(g,\mu )=(R_{g\ast }{f}_1,f_2-{\mathrm{a}\mathrm{d}}_{f_1}^{\ast }\mu ),$(12) where $f_1=f_1(g,\mu )\in \mathfrak{g}$, $f_2=f_2(g,\mu )\in {\mathfrak{g}}^{\ast }$ are the two components of f.

3.2. Symplectic form and Hamiltonian vector fields

The right trivialized² symplectic form pulled back to $\mathbf{G}$ reads (13) $\begin{array}{rl}{\omega }_{(g,\mu )}\left(\right(R_{g\ast }{\xi }_1,\delta {\nu }_1),(R_{g\ast }{\xi }_2,\delta {\nu }_2\left)\right)& =〈\delta {\nu }_2,{\xi }_1〉+\\ & -〈\delta {\nu }_1,{\xi }_2〉-〈\mu ,[{\xi }_1,{\xi }_2]〉,{\xi }_1,{\xi }_2\in \mathfrak{g}.\end{array}$(13) See [32] for more details, proofs and for a the left trivialized symplectic form. The vector field F is a Hamiltonian vector field if it satisfies ${\mathrm{i}}_F\omega =dH,$ for some Hamiltonian function $H:T^{\ast }G\to \mathbb{R}$, where ${\mathrm{i}}_F$ is defined as ${\mathrm{i}}_F\left(X\right):=\omega (F,X)$ for any vector field X. This implies that the map f for such a Hamiltonian vector field gets the form (14) $f(g,\mu )=\left(\frac{\mathrm{\partial }H}{\mathrm{\partial }\mu }(g,\mu ),-R_g^{\ast }\frac{\mathrm{\partial }H}{\mathrm{\partial }g}(g,\mu )\right).$(14) The following is a one-parameter family of symplectic Lie group integrators on $T^{\ast }G$: (15) $\begin{array}{rl}{M}_{\theta }& ={\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_{-\xi }^{\ast }({\mu }_0+{\mathrm{A}\mathrm{d}}_{\exp \left(\theta \xi \right)}^{\ast }(\overline{n}\left)\right)-\theta {\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_{-\theta \xi }^{\ast }{\mathrm{A}\mathrm{d}}_{\exp \left(\theta \xi \right)}^{\ast }\left(\overline{n}\right),\end{array}$(15) (16) $\begin{array}{rl}(\xi ,\overline{n})& =hf(\exp \left(\theta \xi \right)\cdot g_0,M_{\theta }),\end{array}$(16) (17) $\begin{array}{rl}(g_1,{\mu }_1)& =(\exp (\xi ),{\mathrm{A}\mathrm{d}}_{\exp \left(\right(\theta -1\left)\xi \right)}^{\ast }\overline{n})\cdot (g_0,{\mu }_0).\end{array}$(17) For higher order integrators of this type and a complete treatment, see [3].

3.3. Reduced equations Lie Poisson systems

A mechanical system formulated on the cotangent bundle $T^{\ast }G$ with a left or right invariant Hamiltonian can be reduced to a system on ${\mathfrak{g}}^{\ast }$ [36]. In fact for a Hamiltonian H right invariant under the left action of G, $\frac{\mathrm{\partial }H}{\mathrm{\partial }g}=0$, and from (12(12) $F(g,\mu )=T_e{R}_{(g,\mu )}f(g,\mu )=(R_{g\ast }{f}_1,f_2-{\mathrm{a}\mathrm{d}}_{f_1}^{\ast }\mu ),$(12) ) and (14(14) $f(g,\mu )=\left(\frac{\mathrm{\partial }H}{\mathrm{\partial }\mu }(g,\mu ),-R_g^{\ast }\frac{\mathrm{\partial }H}{\mathrm{\partial }g}(g,\mu )\right).$(14) ) we get for the second equation (18) $\dot{\mu }=\mp {\mathrm{a}\mathrm{d}}_{\frac{\mathrm{\partial }H}{\mathrm{\partial }\mu }}^{\ast }\mu ,$(18) where the positive sign is used in case of left invariance (see, e.g. Section 13.4 in [37]). The solution to this system preserves coadjoint orbits, thus using the Lie group action $g\cdot \mu ={\mathrm{A}\mathrm{d}}_{g^{-1}}^{\ast }\mu ,$ to build a Lie group integrator results in preservation of such coadjoint orbits. Lie group integrators for this interesting case were studied in [15].

The Lagrangian counterpart to these Hamiltonian equations are the Euler–Poincaré equations³ [24].

3.4. Three different formulations of the heavy top equations

The heavy top is a simple test example for illustrating the behaviour of Lie group methods. We will consider three different formulations for this mechanical system. The first formulation is on $T^{\ast }SO\left(3\right)$ where the equations are canonical Hamiltonian, a second point of view is that the system is a Lie–Poisson system on $\mathfrak{s}\mathfrak{e}(3{)}^{\ast }$, and finally it is canonical Hamiltonian on a larger group with a quadratic Hamiltonian function. The three different formulations suggest the use of different Lie group integrators.

3.4.1. Heavy top equations on $T^{\ast }SO\left(3\right)$

The heavy top is a rigid body with a fixed point in a gravitational field. The phase space of this mechanical system is $T^{\ast }SO\left(3\right)$ where the equations of the heavy top are in canonical Hamiltonian form. Assuming $(Q,p)$ are coordinates for $T^{\ast }SO\left(3\right)$, $\mathrm{\Pi }=\left(T_e{L}_Q{)}^{\ast }\right(p)$ is the left trivialized or body momentum. The Hamiltonian of the heavy top is given in terms of $(Q,\mathrm{\Pi })$ as $H:SO\left(3\right)⋉\mathfrak{s}\mathfrak{o}(3{)}^{\ast }\to \mathbb{R},H(Q,\mathrm{\Pi })=\frac{1}{2}〈\mathrm{\Pi },{\mathbb{I}}^{-1}\mathrm{\Pi }〉+Mg\ell \mathbf{\Gamma }\cdot \mathcal{X},\mathbf{\Gamma }=Q^{-1}{\mathbf{\Gamma }}_0,$ where $\mathbb{I}:\mathfrak{s}\mathfrak{o}\left(3\right)\to \mathfrak{s}\mathfrak{o}(3{)}^{\ast }$ is the inertia tensor, here represented as a diagonal$3\times 3$ matrix, $\mathbf{\Gamma }=Q^{-1}{\mathbf{\Gamma }}_0$, where ${\mathbf{\Gamma }}_0\in {\mathbb{R}}^3$ is the axis of the spatial coordinate system parallel to the direction of gravity but pointing upwards, M is the mass of the body, g is the gravitational acceleration, $\mathcal{X}$ is the body fixed unit vector of the oriented line segment pointing from the fixed point to the centre of mass of the body, ℓ is the length of this segment. The equations of motion on $SO\left(3\right)⋉\mathfrak{s}\mathfrak{o}(3{)}^{\ast }$ are (19) $\begin{array}{rl}\dot{\mathrm{\Pi }}& =\mathrm{\Pi }\times {\mathbb{I}}^{-1}\mathrm{\Pi }+Mg\ell \mathbf{\Gamma }\times \mathcal{X},\end{array}$(19) (20) $\begin{array}{rl}\dot{Q}& =Q\widehat{{\mathbb{I}}^{-1}\mathrm{\Pi }}.\end{array}$(20) The identification of $T^{\ast }SO\left(3\right)$ with $SO\left(3\right)⋉\mathfrak{s}\mathfrak{o}(3{)}^{\ast }$ via right trivialization leads to the spatial momentum variable $\pi =\left(T_e{R}_Q{)}^{\ast }\right(p)=Q\mathrm{\Pi }$. The equations written in the space variables $(Q,\pi )$ get the form (21) $\begin{array}{rl}\dot{\pi }& =Mg\ell {\mathbf{\Gamma }}_0\times Q\mathcal{X},\end{array}$(21) (22) $\begin{array}{rl}\dot{Q}& =\hat{\omega }Q\omega =Q{\mathbb{I}}^{-1}{Q}^T\pi .\end{array}$(22) where the first equation states that the component of π parallel to ${\mathbf{\Gamma }}_0$ is constant in time. These equations can be obtained from (12(12) $F(g,\mu )=T_e{R}_{(g,\mu )}f(g,\mu )=(R_{g\ast }{f}_1,f_2-{\mathrm{a}\mathrm{d}}_{f_1}^{\ast }\mu ),$(12) ) and (14(14) $f(g,\mu )=\left(\frac{\mathrm{\partial }H}{\mathrm{\partial }\mu }(g,\mu ),-R_g^{\ast }\frac{\mathrm{\partial }H}{\mathrm{\partial }g}(g,\mu )\right).$(14) ) on the right trivialized $T^{\ast }SO\left(3\right)$, $SO\left(3\right)⋉\mathfrak{s}\mathfrak{o}(3{)}^{\ast }$, with the heavy top Hamiltonian and the symplectic Lie group integrators (16(16) $\begin{array}{rl}(\xi ,\overline{n})& =hf(\exp \left(\theta \xi \right)\cdot g_0,M_{\theta }),\end{array}$(16) )–(17(17) $\begin{array}{rl}(g_1,{\mu }_1)& =(\exp (\xi ),{\mathrm{A}\mathrm{d}}_{\exp \left(\right(\theta -1\left)\xi \right)}^{\ast }\overline{n})\cdot (g_0,{\mu }_0).\end{array}$(17) ) can be applied in this case. Similar methods were proposed in [32] and [49].

3.4.2. Heavy top equations on ${\mathfrak{s}\mathfrak{e}}^{\ast }\left(3\right)$

The Hamiltonian of the heavy top is not invariant under the action of $SO\left(3\right)$, so Equations (19(19) $\begin{array}{rl}\dot{\mathrm{\Pi }}& =\mathrm{\Pi }\times {\mathbb{I}}^{-1}\mathrm{\Pi }+Mg\ell \mathbf{\Gamma }\times \mathcal{X},\end{array}$(19) )–(20(20) $\begin{array}{rl}\dot{Q}& =Q\widehat{{\mathbb{I}}^{-1}\mathrm{\Pi }}.\end{array}$(20) ) given in Section 3.4.1 cannot be reduced to ${\mathfrak{s}\mathfrak{o}}^{\ast }\left(3\right)$, nevertheless the heavy top equations are Lie–Poisson on ${\mathfrak{s}\mathfrak{e}}^{\ast }\left(3\right)$, [17,48,52].

Observe that the equations of the heavy top on $T^{\ast }SO\left(3\right)$ (19(19) $\begin{array}{rl}\dot{\mathrm{\Pi }}& =\mathrm{\Pi }\times {\mathbb{I}}^{-1}\mathrm{\Pi }+Mg\ell \mathbf{\Gamma }\times \mathcal{X},\end{array}$(19) )–(20(20) $\begin{array}{rl}\dot{Q}& =Q\widehat{{\mathbb{I}}^{-1}\mathrm{\Pi }}.\end{array}$(20) ) can be easily modified eliminating the variable $Q\in SO\left(3\right)$ and replacing it with $\mathbf{\Gamma }\in {\mathbb{R}}^3$ $\mathbf{\Gamma }=Q^{-1}{\mathbf{\Gamma }}_0$ to obtain (23) $\begin{array}{rl}\dot{\mathrm{\Pi }}& =\mathrm{\Pi }\times {\mathbb{I}}^{-1}\mathrm{\Pi }+Mg\ell \mathbf{\Gamma }\times \mathcal{X},\end{array}$(23) (24) $\begin{array}{rl}\dot{\mathbf{\Gamma }}& =\mathbf{\Gamma }\times \left({\mathbb{I}}^{-1}\mathrm{\Pi }\right).\end{array}$(24) We will see that the solutions of these equations evolve on ${\mathfrak{s}\mathfrak{e}}^{\ast }\left(3\right)$. In what follows, we consider elements of ${\mathfrak{s}\mathfrak{e}}^{\ast }\left(3\right)$ to be pairs of vectors in ${\mathbb{R}}^3$, e.g. $(\mathrm{\Pi },\mathbf{\Gamma })$. Correspondingly the elements of $SE\left(3\right)$ are represented as pairs $(g,\mathbf{u})$ with $g\in SO\left(3\right)$ and $\mathbf{u}\in {\mathbb{R}}^3$. The group multiplication in $SE\left(3\right)$ is then $(g_1,{\mathbf{u}}_1)\cdot (g_2,{\mathbf{u}}_2)=(g_1{g}_2,g_1{\mathbf{u}}_2+{\mathbf{u}}_1),$ where $g_1{g}_2$ is the product in $SO\left(3\right)$ and $g_1\mathbf{u}$ is the product of a $3\times 3$ orthogonal matrix with a vector in ${\mathbb{R}}^3$. The coadjoint representation and its infinitesimal generator on ${\mathfrak{s}\mathfrak{e}}^{\ast }\left(3\right)$ take the form ${\mathrm{A}\mathrm{d}}_{(g,\mathbf{u})}^{\ast }(\mathrm{\Pi },\mathbf{\Gamma })=\left(g^{-1}\right(\mathrm{\Pi }-\mathbf{u}\times \mathbf{\Gamma }),g^{-1}\mathbf{\Gamma }),{\mathrm{a}\mathrm{d}}_{(\xi ,\mathbf{u})}^{\ast }(\mathrm{\Pi },\mathbf{\Gamma })=(-\xi \times \mathrm{\Pi }-\mathbf{u}\times \mathbf{\Gamma },-\xi \times \mathbf{\Gamma }).$ Using this expression for ${\mathrm{a}\mathrm{d}}_{(\xi ,\mathbf{u})}^{\ast }$ with $(\xi =\frac{\mathrm{\partial }H}{\mathrm{\partial }\mathrm{\Pi }},\mathbf{u}=\frac{\mathrm{\partial }H}{\mathrm{\partial }\mathbf{\Gamma }})$, it can be easily seen that Equation (18(18) $\dot{\mu }=\mp {\mathrm{a}\mathrm{d}}_{\frac{\mathrm{\partial }H}{\mathrm{\partial }\mu }}^{\ast }\mu ,$(18) ) in this setting reproduce the heavy top Equations (23(23) $\begin{array}{rl}\dot{\mathrm{\Pi }}& =\mathrm{\Pi }\times {\mathbb{I}}^{-1}\mathrm{\Pi }+Mg\ell \mathbf{\Gamma }\times \mathcal{X},\end{array}$(23) )–(24(24) $\begin{array}{rl}\dot{\mathbf{\Gamma }}& =\mathbf{\Gamma }\times \left({\mathbb{I}}^{-1}\mathrm{\Pi }\right).\end{array}$(24) ). Therefore the equations are Lie–Poisson equations on ${\mathfrak{s}\mathfrak{e}}^{\ast }\left(3\right)$. However, since the heavy top is a rigid body with a fixed point and there are no translations, these equations do not arise from a reduction of $T^{\ast }SE\left(3\right)$. Moreover, the Hamiltonian on $\mathfrak{s}\mathfrak{e}(3{)}^{\ast }$ is not quadratic and the equations are not geodesic equations. Implicit and explicit Lie group integrators applicable to this formulation of the heavy top equations and preserving coadjoint orbits were discussed in [15], for a variable stepsize integrator applied to this formulation of the heavy top, see [12].

3.4.3. Heavy top equations with quadratic Hamiltonian

We rewrite the heavy top equations one more time considering the constant vector $\mathbf{p}=-Mg\ell \mathcal{X}$ as a momentum variable conjugate to the position $\mathbf{q}\in {\mathbb{R}}^3$ and where $\mathbf{p}=Q^{-1}{\mathbf{\Gamma }}_0+\dot{\mathbf{q}}$, and the Hamiltonian is a quadratic function of Π, Q, $\mathbf{p}$ and $\mathbf{q}$: $\begin{array}{rl}& H:T^{\ast }SO\left(3\right)\times {{\mathbb{R}}^3}^{\ast }\times {\mathbb{R}}^3\to \mathbb{R},\\ & H\left(\right(\mathrm{\Pi },Q),(\mathbf{p},\mathbf{q}\left)\right)=\frac{1}{2}〈\mathrm{\Pi },{\mathbb{I}}^{-1}\mathrm{\Pi }〉+\frac{1}{2}\parallel \mathbf{p}-Q^{-1}{\mathbf{\Gamma }}_0{\parallel }^2-\frac{1}{2}\parallel Q^{-1}{\mathbf{\Gamma }}_0{\parallel }^2,\end{array}$ see [23, section 8.5]. This Hamiltonian is invariant under the left action of $SO\left(3\right)$. The corresponding equations are canonical on $T^{\ast }S\equiv S⋉{\mathfrak{s}}^{\ast }$ where $S=SO\left(3\right)\times {\mathbb{R}}^3$ with Lie algebra $\mathfrak{s}:=\mathfrak{s}\mathfrak{o}\left(3\right)\times {\mathbb{R}}^3$ and $T^{\ast }S$ can be identified with $T^{\ast }SO\left(3\right)\times {{\mathbb{R}}^3}^{\ast }\times {\mathbb{R}}^3$. The equations are (25) $\begin{array}{rl}\dot{\mathrm{\Pi }}& =\mathrm{\Pi }\times {\mathbb{I}}^{-1}\mathrm{\Pi }-\left(Q^{-1}{\mathbf{\Gamma }}_0\right)\times \mathbf{p},\end{array}$(25) (26) $\begin{array}{rl}\dot{Q}& =Q\widehat{{\mathbb{I}}^{-1}\mathrm{\Pi }},\end{array}$(26) (27) $\begin{array}{rl}\dot{\mathbf{p}}& =\mathbf{0},\end{array}$(27) (28) $\begin{array}{rl}\dot{\mathbf{q}}& =\mathbf{p}-Q^{-1}{\mathbf{\Gamma }}_0.\end{array}$(28) and in the spatial momentum variables (29) $\begin{array}{rl}\dot{\pi }& =-{\mathbf{\Gamma }}_0\times Q\mathbf{p},\end{array}$(29) (30) $\begin{array}{rl}\dot{Q}& =\hat{\omega }Q,\omega =Q{\mathbb{I}}^{-1}{Q}^T\pi ,\end{array}$(30) (31) $\begin{array}{rl}\dot{\mathbf{p}}& =\mathbf{0},\end{array}$(31) (32) $\begin{array}{rl}\dot{\mathbf{q}}& =\mathbf{p}-Q^{-1}{\mathbf{\Gamma }}_0.\end{array}$(32) Similar formulations were considered in [31] for the stability analysis of an underwater vehicle. A similar but different formulation of the heavy top was considered in [4].

3.4.4. Numerical experiments

We apply various implicit Lie group integrators to the heavy top system. The test problem we consider is the same as in [4], where $Q\left(0\right)=I$, $\ell =2$, M = 15 $\mathbb{I}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(0.234375,0.46875,0.234375)$, $\pi \left(0\right)=\mathbb{I}(0,150,-4.61538)$, $\mathcal{X}=(0,1,0)$ ${\mathrm{\Gamma }}_0=(0,0,-9.81)$.

In Figure 2, we report the performance of the symplectic Lie group integrators (15(15) $\begin{array}{rl}{M}_{\theta }& ={\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_{-\xi }^{\ast }({\mu }_0+{\mathrm{A}\mathrm{d}}_{\exp \left(\theta \xi \right)}^{\ast }(\overline{n}\left)\right)-\theta {\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_{-\theta \xi }^{\ast }{\mathrm{A}\mathrm{d}}_{\exp \left(\theta \xi \right)}^{\ast }\left(\overline{n}\right),\end{array}$(15) )–(17(17) $\begin{array}{rl}(g_1,{\mu }_1)& =(\exp (\xi ),{\mathrm{A}\mathrm{d}}_{\exp \left(\right(\theta -1\left)\xi \right)}^{\ast }\overline{n})\cdot (g_0,{\mu }_0).\end{array}$(17) ) applied both on Equations (21(21) $\begin{array}{rl}\dot{\pi }& =Mg\ell {\mathbf{\Gamma }}_0\times Q\mathcal{X},\end{array}$(21) )–(22(22) $\begin{array}{rl}\dot{Q}& =\hat{\omega }Q\omega =Q{\mathbb{I}}^{-1}{Q}^T\pi .\end{array}$(22) ) with $\theta =0$, $\theta =\frac{1}{2}$ and $\theta =1$ (SLGI), and to Equations (29(29) $\begin{array}{rl}\dot{\pi }& =-{\mathbf{\Gamma }}_0\times Q\mathbf{p},\end{array}$(29) )–(32(32) $\begin{array}{rl}\dot{\mathbf{q}}& =\mathbf{p}-Q^{-1}{\mathbf{\Gamma }}_0.\end{array}$(32) ) with $\theta =\frac{1}{2}$ (SLGIKK). The methods with $\theta =\frac{1}{2}$ attain order 2. In Figure 3, we show the energy error for the symplectic Lie group integrators with $\theta =\frac{1}{2}$ and $\theta =0$ integrating with stepsize h = 0.01 for 6000 steps.

Figure 1. Illustration of the heavy top, where CM is the centre of mass of the body, O is the fixed point, $\overrightarrow{g}$ is the gravitational acceleration vector, and $\ell ,Q,\overrightarrow{\chi }$ follow the notation introduced in Section 3.4.1.

Figure 2. Symplectic Lie group integrators integration on the time interval $[0,1]$. Left: 3D plot of $M\ell Q^{-1}{\mathrm{\Gamma }}_0$. Centre: components of $Q\mathcal{X}$. The left and centre plots are computed with the same step-size. Right: verification of the order of the methods.

Figure 3. Symplectic Lie group integrators, long time integration, h = 0.01, 6000 steps.. Top: energy error, bottom 3D plot of $M\ell Q^{-1}{\mathrm{\Gamma }}_0$.

4. Variable step size

One approach for varying the step size is based on the use of an embedded Runge–Kutta pair. This principle can be carried from standard Runge–Kutta methods in vector spaces to the present situation with RKMK and commutator-free schemes via minor modifications. We briefly summarize the main principle of embedded pairs before giving more specific details for the case of Lie group integrators. This approach is very well documented in the literature and goes back to Merson [38] and a detailed treatment can be found in [19, pp. 165–168].

An embedded pair consists of a main method used to propagate the numerical solution, together with some auxiliary method that is only used to obtain an estimate of the local error. This local error estimate is in turn used to derive a step size adjustment formula that attempts to keep the local error estimate approximately equal to some user defined tolerance $\mathrm{t}\mathrm{o}\mathrm{l}$ in every step. Suppose the main method is of order p and the auxiliary method is of order $\tilde{p}\ne p.$⁴ Both methods are applied to the input value $y_n$ and yields approximations $y_{n+1}$ and ${\tilde{y}}_{n+1}$ respectively, using the same step size $h_{n+1}$. Now, some distance measure⁵ between $y_{n+1}$ and ${\tilde{y}}_{n+1}$ provides an estimate $e_{n+1}$ for the size of the local truncation error. Thus $e_{n+1}=C{h}_{n+1}^{\tilde{p}+1}+\mathcal{O}\left(h^{\tilde{p}+2}\right)$. Aiming at $e_{n+1}\approx \mathrm{t}\mathrm{o}\mathrm{l}$ in every step, one may use a formula of the type (33) $h_{n+1}=\theta {\left(\frac{\mathrm{t}\mathrm{o}\mathrm{l}}{e_{n+1}}\right)}^{\frac{1}{\tilde{p}+1}}{h}_n$(33) where θ is a ‘safety factor’, typically chosen between 0.8 and 0.9. In case the step is rejected because $e_n>\mathrm{t}\mathrm{o}\mathrm{l}$ we can redo the step with a step size obtained by the same formula. We summarize the approach in the following algorithm

Here we have used again the safety factor θ, and the parameter α is generally chosen as $\alpha =\frac{1}{1+min(p,\tilde{p})}$.

4.1. RKMK methods with variable stepsize

We need to specify how to calculate the quantity $e_{n+1}$ in each step. For RKMK methods, the situation is simplified by the fact that we are solving the local problem (6(6) $\dot{\sigma }\left(t\right)={\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_{\sigma \left(t\right)}^{-1}\left(f\right(\exp \left(\sigma \right(t\left)\right)\cdot y_0\left)\right)$(6) ) in the linear space $\mathfrak{g}$, where the known theory can be applied directly. So any standard embedded pair of Runge–Kutta methods described by coefficients $(a_{ij},b_i,{\tilde{a}}_{ij},{\tilde{b}}_i)$ of orders $(p,\tilde{p})$ can be applied to the full dexpinv-equation (6(6) $\dot{\sigma }\left(t\right)={\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_{\sigma \left(t\right)}^{-1}\left(f\right(\exp \left(\sigma \right(t\left)\right)\cdot y_0\left)\right)$(6) ) to obtain local Lie algebra approximations ${\sigma }_1$, ${\tilde{\sigma }}_1$ and one uses, e.g. $e_{n+1}=\parallel {\sigma }_1-{\tilde{\sigma }}_1\parallel$ (note that the equation itself depends on $y_n$). For methods which use a truncated version of the series for ${\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_u^{-1}$ one may also try to optimize performance by including commutators that are shared between the main method and the auxiliary scheme.

4.2. Commutator-free methods with variable stepsize

For the commutator-free methods of Section 2.2, the situation is different since such methods do not have a natural local representation in a linear space. One can still derive embedded pairs, and this can be achieved by studying order conditions [44] as was done in [12]. Now one obtains after each step two approximations $y_{n+1}$ and ${\tilde{y}}_{n+1}$ on $\mathcal{M}$ both by using the same initial value $y_n$ and step size $h_n$. One must also have access to some metric d to calculate $e_{n+1}=d(y_{n+1},{\tilde{y}}_{n+1})$ We give a few examples of embedded pairs.

4.2.1. Pairs of order $(p,\tilde{p})=(3,2)$

It is possible to obtain embedded pairs of order 3(2) which satisfy the requirements above. We present two examples from [12]. The first one reuses the second stage exponential in the update $\begin{array}{rl}{Y}_{n,1}& =y_n\\ Y_{n,2}& =\exp \left(\frac{1}{3}h{f}_{n,1}\right)\cdot y_n\\ Y_{n,3}& =\exp \left(\frac{2}{3}h{f}_{n,2}\right)\cdot y_n\\ y_{n+1}& =\exp \left(h\right(-\frac{1}{12}{f}_{n,1}+\frac{3}{4}{f}_{n,3}\left)\right)\cdot Y_{n,2}\\ {\tilde{y}}_{n+1}& =\exp \left(\frac{1}{2}h\right(f_{n,2}+f_{n,3}\left)\right)\cdot y_n\end{array}$ One could also have reused the third stage $Y_{n,3}$ in the update, rather than $Y_{n,2}$. $\begin{array}{rl}{Y}_{n,1}& =y_n\\ Y_{n,2}& =\exp \left(\frac{2}{3}h{f}_{n,1}\right)\cdot y_n\\ Y_{n,3}& =\exp \left(h\right(\frac{5}{12}{f}_{n,1}+\frac{1}{4}{f}_{n,2})\cdot y_n\\ y_{n+1}& =\exp \left(h\right(-\frac{1}{6}{f}_{n,1}-\frac{1}{2}{f}_{n,2}+f_{n,3}\left)\right)\cdot Y_{n,3}\\ {\tilde{y}}_{n+1}& =\exp \left(\frac{1}{4}h\right(f_{n,1}+3f_{n,3}\left)\right)\cdot y_n\end{array}$ It is always understood that the frozen vector fields are $f_{n,i}:=f_{Y_{n,i}}$.

4.2.2. Order $(4,3)$

The procedure of deriving efficient pairs becomes more complicated as the order increases. In [12], a low cost pair of order $(4,3)$ was derived, in the sense that one attempted to minimize the number of stages and exponentials in the embedded pair as a whole. This came, however, at the expense of a relatively large error constant. So rather than presenting the method from that paper, we suggest a simpler procedure at the cost of some more computational work per step, we simply furnish the commutator-free method of Section 2 by a third-order auxiliary scheme. It can be described as follows:

1. Compute $Y_{n,i},i=1\dots ,4$ and $y_{n+1}$ from (9(9) $\begin{array}{r}\begin{array}{rl}{Y}_{n,1}& =y_n\\ Y_{n,2}& =\exp \left(\frac{1}{2}h{f}_{n,1}\right)\cdot y_n\\ Y_{n,3}& =\exp \left(\frac{1}{2}h{f}_{n,2}\right)\cdot y_n\\ Y_{n,4}& =\exp (h{f}_{n,3}-\frac{1}{2}h{f}_{n,1})\cdot Y_{n,2}\\ y_{n+\frac{1}{2}}& =\exp \left(\frac{1}{12}h\right(3f_{n,1}+2f_{n,2}+2f_{n,3}-f_{n,4}\left)\right)\cdot y_n\\ y_{n+1}& =\exp \left(\frac{1}{12}h\right(-f_{n,1}+2f_{n,2}+2f_{n,3}+3f_{n,4}\left)\right)\cdot y_{n+\frac{1}{2}}\end{array}\end{array}$(9) )

2. Compute an additional stage ${\overline{Y}}_{n,3}$ and then ${\tilde{y}}_{n+1}$ as (34) $\begin{array}{r}\begin{array}{rl}{\overline{Y}}_{n,3}& =\exp \left(\frac{3}{4}h{f}_{n,2}\right)\cdot y_n\\ {\tilde{y}}_{n+1}& =\exp \left(\frac{h}{9}\right(-f_{n,1}+3f_{n,2}+4{\overline{f}}_{n,3}\left)\right)\cdot \exp \left(\frac{h}{3}{f}_{n,1}\right)\cdot y_n\end{array}\end{array}$(34)

5. The N-fold 3D pendulum

In this section, we present a model for a system of N connected three-dimensional pendulums. The modelling part comes from [28], and here we study the vector field describing the dynamics, in order to re-frame it into the Lie group integrators setting described in the previous sections. The model we use is not completely realistic since, for example, it neglects possible interactions between pendulums, and it assumes ideal spherical joints between them. However, this is still a relevant example from the point of view of geometric numerical integration. More precisely, we show a possible way to work with a configuration manifold which is not a Lie group, applying the theoretical instruments introduced before.

The Lagrangian we consider is a function from $(T{S}^2{)}^N$ to $\mathbb{R}$. Instead of the coordinates $(q_1,\dots ,q_N,{\dot{q}}_1,\dots ,{\dot{q}}_N)$, where ${\dot{q}}_i\in T_{q_i}{S}^2$, we choose to work with the angular velocities. Precisely, $T_{q_i}{S}^2=\{v\in {\mathbb{R}}^3:v^T{q}_i=0\}=〈q_i{〉}^{\perp }\subset {\mathbb{R}}^3,$ and hence for any ${\dot{q}}_i\in T_{q_i}{S}^2$ there exist ${\omega }_i\in {\mathbb{R}}^3$ such that ${\dot{q}}_i={\omega }_i\times q_i$, which can be interpreted as the angular velocity of $q_i$. So we can assume without loss of generality that ${\omega }_i^T{q}_i=0$ (i.e. ${\omega }_i\in T_{q_i}{S}^2$) and pass to the coordinates $(q_1,{\omega }_1,q_2,{\omega }_2,\dots ,q_N,{\omega }_N)\in (T{S}^2{)}^N$ to describe the dynamics. In this section, we denote with $m_1,\dots ,m_N$ the masses of the pendulums and with $L_1,\dots ,L_N$ their lengths. Figure 4 shows the case N = 3. We organize the section into three parts:

1. We define the transitive Lie group action used to integrate this model numerically,

2. We show a possible way to express the dynamics in terms of the infinitesimal 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 integrators and some classical numerical integrator. Then we conclude with numerical experiments on variable step size.

5.1. Transitive group action on $(T{S}^2{)}^N$

We characterize a transitive action for $(T{S}^2{)}^N$, starting with the case N = 1 and generalizing it to N>1. The action we consider is based on the identification between $\mathfrak{s}\mathfrak{e}\left(3\right)$, the Lie algebra of $SE\left(3\right)$ and ${\mathbb{R}}^6$. We start from the Ad-action of $SE\left(3\right)$ on $\mathfrak{s}\mathfrak{e}\left(3\right)$ (see [23]), which writes $\begin{array}{rl}& \mathrm{A}\mathrm{d}:SE\left(3\right)\times \mathfrak{s}\mathfrak{e}\left(3\right)\to \mathfrak{s}\mathfrak{e}\left(3\right),\\ & \mathrm{A}\mathrm{d}\left(\right(R,r),(u,v\left)\right)=(Ru,Rv+\hat{r}Ru).\end{array}$ Since $\mathfrak{s}\mathfrak{e}\left(3\right)\simeq {\mathbb{R}}^6$, the Ad-action allows us to define the following Lie group action on ${\mathbb{R}}^6$ $\psi :SE\left(3\right)\times {\mathbb{R}}^6\to {\mathbb{R}}^6,\psi \left(\right(R,r),(u,v\left)\right)=(Ru,Rv+\hat{r}Ru).$ We can think of ψ as a Lie group action on $T{S}^2$ since, for any $q\in {\mathbb{R}}^3$, it maps points of $T{S}_{|q|}^2:=\left\{\right(\stackrel{\sim }{q},\stackrel{\sim }{\omega })\in {\mathbb{R}}^3\times {\mathbb{R}}^3:{\stackrel{\sim }{\omega }}^T\stackrel{\sim }{q}=0,|\stackrel{\sim }{q}|=|q|\}\subset {\mathbb{R}}^6$ into other points of $T{S}_{|q|}^2$. Moreover, with standard arguments (see [43]), it is possible to prove that the orbit of a generic point $m=(q,\omega )\in {\mathbb{R}}^6$ with ${\omega }^{T}q=0$ coincides with $\mathrm{O}\mathrm{r}\mathrm{b}\left(m\right)=T{S}_{|q|}^2.$ In particular, when $q\in {\mathbb{R}}^3$ is a unit vector (i.e. $q\in S^2$), ψ allows us to define a transitive Lie group action on $T{S}^2=T{S}_{|q|=1}^2$ which writes $\begin{array}{rl}& \psi :SE\left(3\right)\times T{S}^2\to T{S}^2\\ & \psi \left(\right(A,a),(q,\omega \left)\right):={\psi }_{(A,a)}(q,\omega )=(Aq,A\omega +\hat{a}Aq)=(\overline{q},\overline{\omega }).\end{array}$ To conclude the description of the action, we report here its infinitesimal generator which is fundamental in the Lie group integrators setting ${\left.{\psi }_{\ast }\left(\right(u,v\left)\right)\right|}_{(q,\omega )}=(\hat{u}q,\hat{u}\omega +\hat{v}q).$ We can extend this construction to the case N>1 in a natural way, i.e. through the action of a Lie group obtained from cartesian products of $SE\left(3\right)$ and equipped with the direct product structure. More precisely, we consider the group $G=\left(SE\right(3){)}^N$ and by direct product structure we mean that for any pair of elements ${\delta }^{\left(1\right)}=({\delta }_1^{\left(1\right)},\dots ,{\delta }_N^{\left(1\right)}),{\delta }^{\left(2\right)}=({\delta }_1^{\left(2\right)},\dots ,{\delta }_N^{\left(2\right)})\in G,$ denoted with $\ast$ the semidirect product of $SE\left(3\right)$, we define the product ° on G as ${\delta }^{\left(1\right)}\circ {\delta }^{\left(2\right)}:=({\delta }_1^{\left(1\right)}\ast {\delta }_1^{\left(2\right)},\dots ,{\delta }_N^{\left(1\right)}\ast {\delta }_N^{\left(2\right)})\in G.$ With this group structure defined, we can generalize the action introduced for N = 1 to larger Ns as follows: $\begin{array}{rl}& \psi :\left(SE\right(3){)}^N\times (T{S}^2{)}^N\to (T{S}^2{)}^N,\\ & \psi \left(\right(A_1,a_1,\dots ,A_N,a_n),(q_1,{\omega }_1,\dots ,q_N,{\omega }_N\left)\right)\\ & =(A_1{q}_1,A_1{\omega }_1+{\hat{a}}_1{A}_1{q}_1,\dots ,A_N{q}_N,A_N{\omega }_N+{\hat{a}}_N{A}_N{q}_N),\end{array}$ whose infinitesimal generator writes ${\psi }_{\ast }\left(\xi \right){|}_m=({\hat{u}}_1{q}_1,{\hat{u}}_1{\omega }_1+{\hat{v}}_1{q}_1,\dots ,{\hat{u}}_N{q}_N,{\hat{u}}_N{\omega }_N+{\hat{v}}_N{q}_N),$ where $\xi =[u_1,v_1,\dots ,u_N,v_N]\in \mathfrak{s}\mathfrak{e}(3{)}^N$ and $m=(q_1,{\omega }_1,\dots ,q_N,{\omega }_N)\in (T{S}^2{)}^N$. We have now the only group action we need to deal with the N-fold spherical pendulum. In the following part of this section, we work on the vector field describing the dynamics and adapt it to the Lie group integrators setting.

Figure 4. Threefold pendulum at a fixed time instant, with fixed point placed at the origin.

5.2. Full chain

We consider the vector field $F\in \mathfrak{X}\left(\right(T{S}^2{)}^N)$, describing the dynamics of the N-fold 3D pendulum, and we express it in terms of the infinitesimal generator of the action defined above. More precisely, we find a function $F:(T{S}^2{)}^N\to \mathfrak{s}\mathfrak{e}(3{)}^N$ such that ${\psi }_{\ast }\left(f\right(m\left)\right){|}_m=F{|}_m,\mathrm{\forall }m\in (T{S}^2{)}^N.$ We omit the derivation of F starting from the Lagrangian of the system, which can be found in the section devoted to mechanical systems on $(S^2{)}^N$ of [28]. The configuration manifold of the system is $(S^2{)}^N$, while the Lagrangian, expressed in terms of the variables $(q_1,{\omega }_1,\dots ,q_N,{\omega }_N)\in (T{S}^2{)}^N$, writes $L(q,\omega )=T(q,\omega )-U\left(q\right)=\frac{1}{2}\sum _{i,j=1}^N\left(M_{ij}{\omega }_i^T{\hat{q}}_i^T{\hat{q}}_j{\omega }_j\right)-\sum _{i=1}^N\left(\sum _{j=i}^N{m}_j\right)g{L}_i{e}_3^T{q}_i,$ where $M_{ij}=\left(\sum _{k=\mathrm{m}\mathrm{a}\mathrm{x}\{i,j\}}^N{m}_k\right)L_i{L}_j{I}_3\in {\mathbb{R}}^{3\times 3}$ is the inertia matrix of the system, $I_3$ is the $3\times 3$ identity matrix, and $e_3=[0,0,1{]}^T$. Noticing that when i = j we get ${\omega }_i^T{\hat{q}}_i^T{\hat{q}}_i{\omega }_i={\omega }_i^T(I_3-q_i{q}_i^T){\omega }_i={\omega }_i^T{\omega }_i,$ we simplify the notation writing $T(q,\omega )=\frac{1}{2}\sum _{i,j=1}^N({\omega }_i^{T}R(q{)}_{ij}{\omega }_j)$ where $R\left(q\right)\in {\mathbb{R}}^{3N\times 3N}$ is a symmetric block matrix defined as $\begin{array}{rl}R(q{)}_{ii}& =\left(\sum _{j=i}^N{m}_j\right)L_i^2{I}_3\in {\mathbb{R}}^{3\times 3},\\ R(q{)}_{ij}& =\left(\sum _{k=j}^N{m}_k\right)L_i{L}_j{\hat{q}}_i^T{\hat{q}}_j\in {\mathbb{R}}^{3\times 3}=R(q{)}_{ji}^T,i<j.\end{array}$ The vector field on which we need to work defines the following first-order ODE: $\begin{array}{rl}{\dot{q}}_i& ={\omega }_i\times q_i,i=1,\dots ,N,\\ R\left(q\right)\dot{\omega }& ={\left[\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N{M}_{ij}|{\omega }_j{|}^2{\hat{q}}_i{q}_j-\left(\sum _{j=i}^N{m}_j\right)g{L}_i{\hat{q}}_i{e}_3\right]}_{i=1,\dots ,N}\in {\mathbb{R}}^{3N}\end{array}$ By direct computation, it is possible to see that, for any $q=(q_1,\dots ,q_N)\in (S^2{)}^N$ and $\omega \in T_{q_1}{S}^2\times \cdots \times T_{q_N}{S}^2$, we have $\left(R\right(q)\omega {)}_i\in T_{q_i}{S}^2.$ Therefore, there is a well-defined linear map $A_q:T_{q_1}{S}^2\times \cdots \times T_{q_N}{S}^2\to T_{q_1}{S}^2\times \cdots \times T_{q_N}{S}^2,A_q\left(\omega \right):=R\left(q\right)\omega .$ We can even notice that $R\left(q\right)$ defines a positive-definite bilinear form on this linear space, since ${\omega }^{T}R\left(q\right)\omega =\sum _{i,j=1}^N{\omega }_i^T{\hat{q}}_i^T{M}_{ij}{\hat{q}}_j{\omega }_j=\sum _{i,j=1}^N\left({\hat{q}}_i{\omega }_i{)}^T{M}_{ij}\right({\hat{q}}_j{\omega }_j)=v^{T}Mv>0.$ The last inequality holds because M is the inertia matrix of the system and hence it defines a symmetric positive-definite bilinear form on $T_{q_1}{S}^2\times \cdots \times T_{q_N}{S}^2$, see, e.g. [16].⁶ This implies the map $A_q$ is invertible and hence we are ready to express the vector field in terms of the infinitesimal generator. We can rewrite the ODEs for the angular velocities as follows: $\dot{\omega }=A_q^{-1}([g_1,\dots ,g_N{]}^T)=\left[\begin{array}{c}{h}_1(q,\omega )\\ \cdots \\ h_N(q,\omega )\end{array}\right]=\left[\begin{array}{c}{a}_1(q,\omega )\times q_1\\ \cdots \\ a_N(q,\omega )\times q_N\end{array}\right]$ where $g_i=g_i(q,\omega )=\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^{N}M(q{)}_{ij}|{\omega }_j{|}^2{\hat{q}}_i{q}_j-\left(\sum _{j=i}^N{m}_j\right)g{L}_i{\hat{q}}_i{e}_3,i=1,\dots ,N$ and $a_1,\dots ,a_N:(T{S}^2{)}^N\to {\mathbb{R}}^3$ are N functions whose existence is guaranteed by the analysis done above. Indeed, we can set $a_i(q,\omega ):=q_i\times h_i(q,\omega )$ and conclude that a mapping f from $(T{S}^2{)}^N$ to $\left(\mathfrak{s}\mathfrak{e}\right(3){)}^N$ such that ${\psi }_{\ast }\left(f\right(q,\omega \left)\right){|}_{(q,\omega )}=F{|}_{(q,\omega )}$ is the following one: $f(q,\omega )=\left[\begin{array}{c}{\omega }_1\\ q_1\times h_1\\ \cdots \\ \cdots \\ {\omega }_N\\ q_N\times h_N\end{array}\right]\in \mathfrak{s}\mathfrak{e}(3{)}^N\simeq {\mathbb{R}}^{6N}.$ We will not go into the Hamiltonian formulation of this problem; however, we remark that a similar approach works even in that situation. Indeed, following the derivation presented in [28], we see that for a mechanical system on $(S^2{)}^N$ the conjugate momentum writes $T_{q_1}^{\ast }{S}^2\times \cdots T_{q_N}^{\ast }{S}^2\ni \pi =({\pi }_1,\dots ,{\pi }_N),\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{\pi }_i=-{\hat{q}}_i^2\frac{\mathrm{\partial }L}{\mathrm{\partial }{\omega }_i}$ and its components are still orthogonal to the respective base points $q_i\in S^2$. Moreover, Hamilton's equations take the form $\begin{array}{rl}{\dot{q}}_i& =\frac{\mathrm{\partial }H(q,\pi )}{\mathrm{\partial }{\pi }_i}\times q_i,\\ {\dot{\pi }}_i& =\frac{\mathrm{\partial }H(q,\pi )}{\mathrm{\partial }{q}_i}\times q_i+\frac{\mathrm{\partial }H(q,\pi )}{\mathrm{\partial }{\pi }_i}\times {\pi }_i,\end{array}$ which implies that setting $f(q,\pi )=\left[\begin{array}{ccccc}{\mathrm{\partial }}_{q_1}H(q,\pi ),& {\mathrm{\partial }}_{{\pi }_1}H(q,\pi ),& \dots ,& {\mathrm{\partial }}_{q_N}H(q,\pi ),& {\mathrm{\partial }}_{{\pi }_N}H(q,\pi )\end{array}\right]$ we can represent even the Hamiltonian vector field of the N-fold 3D pendulum in terms of this group action.

5.2.1. Case N = 2

We have seen how it is possible to turn the equations of motion of a N-chain of pendulums into the Lie group integrators setting. Now we focus on the example with N = 2 pendulums. The equations of motion write (35) $\begin{array}{rl}{\dot{q}}_1& ={\hat{\omega }}_1{q}_1,{\dot{q}}_2={\hat{\omega }}_2{q}_2,\\ R\left(q\right)\left[\begin{array}{c}{\dot{\omega }}_1\\ {\dot{\omega }}_2\end{array}\right]& =\left[\begin{array}{c}(-m_2{L}_1{L}_2|{\omega }_2{|}^2{\hat{q}}_2+(m_1+m_2\left)g{L}_1{\hat{e}}_3\right)q_1\\ (-m_2{L}_1{L}_2|{\omega }_1{|}^2{\hat{q}}_1+m_{2}g{L}_2{\hat{e}}_3)q_2\end{array}\right],\end{array}$(35) where $R\left(q\right)=\left[\begin{array}{cc}(m_1+m_2)L_1^2{I}_3& m_2{L}_1{L}_2{\hat{q}}_1^T{\hat{q}}_2\\ m_2{L}_1{L}_2{\hat{q}}_2^T{\hat{q}}_1& m_2{L}_2^2{I}_3\end{array}\right].$ As presented above, the matrix $R\left(q\right)$ defines a linear invertible map of the space $T_{q_1}{S}^2\times T_{q_2}{S}^2$ onto itself: $A_{(q_1,q_2)}:T_{q_1}{S}^2\times T_{q_2}{S}^2\to T_{q_1}{S}^2\times T_{q_2}{S}^2,[{\omega }_1,{\omega }_2{]}^T\to R(q\left)\right[{\omega }_1,{\omega }_2{]}^T.$ We can easily see that it is well defined since $R\left(q\right)\left[\begin{array}{c}{\omega }_1\\ {\omega }_2\end{array}\right]=\left[\begin{array}{cc}(m_1+m_2)L_1^2{I}_3& m_2{L}_1{L}_2{\hat{q}}_1^T{\hat{q}}_2\\ m_2{L}_1{L}_2{\hat{q}}_2^T{\hat{q}}_1& m_2{L}_2^2{I}_3\end{array}\right]\left[\begin{array}{c}{\hat{v}}_1{q}_1\\ {\hat{v}}_2{q}_2\end{array}\right]=\left[\begin{array}{c}{\hat{r}}_1{q}_1\\ {\hat{r}}_2{q}_2\end{array}\right]\in (T{S}^2{)}^2$ with $\begin{array}{rl}{r}_1(q,\omega )& :=(m_1+m_2)L_1^2{v}_1+m_2{L}_1{L}_2{\hat{q}}_2{\hat{v}}_2{q}_2,\\ r_2(q,\omega )& :=m_2{L}_1{L}_2{\hat{q}}_1{\hat{v}}_1{q}_1+m_2{L}_2^2{v}_2.\end{array}$ This map guarantees that if we rewrite the pair of equations for the angular velocities in (35(35) $\begin{array}{rl}{\dot{q}}_1& ={\hat{\omega }}_1{q}_1,{\dot{q}}_2={\hat{\omega }}_2{q}_2,\\ R\left(q\right)\left[\begin{array}{c}{\dot{\omega }}_1\\ {\dot{\omega }}_2\end{array}\right]& =\left[\begin{array}{c}(-m_2{L}_1{L}_2|{\omega }_2{|}^2{\hat{q}}_2+(m_1+m_2\left)g{L}_1{\hat{e}}_3\right)q_1\\ (-m_2{L}_1{L}_2|{\omega }_1{|}^2{\hat{q}}_1+m_{2}g{L}_2{\hat{e}}_3)q_2\end{array}\right],\end{array}$(35) ) as $\begin{array}{rl}\dot{\omega }& =R^{-1}\left(q\right)\left[\begin{array}{c}(-m_2{L}_1{L}_2|{\omega }_2{|}^2{\hat{q}}_2+(m_1+m_2\left)g{L}_1{\hat{e}}_3\right)q_1\\ (-m_2{L}_1{L}_2|{\omega }_1{|}^2{\hat{q}}_1+m_{2}g{L}_2{\hat{e}}_3)q_2\end{array}\right]=R^{-1}\left(q\right)b=\\ & =A_{(q_1,q_2)}^{-1}\left(b\right)=\left[\begin{array}{c}{h}_1\\ h_2\end{array}\right]\in T_{q_1}{S}^2\times T_{q_2}{S}^2,\end{array}$ then we are assured that there exists a pair of functions $a_1,a_2:T{S}^2\times T{S}^2\to {\mathbb{R}}^3$ such that $\dot{\omega }=\left[\begin{array}{c}{a}_1(q,\omega )\times q_1\\ a_2(q,\omega )\times q_2\end{array}\right]=\left[\begin{array}{c}{h}_1\left(q\right)\\ h_2\left(q\right)\end{array}\right].$ Since we want $a_i\times q_i=h_i$, we just impose $a_i=q_i\times h_i$ and hence the whole vector field can be rewritten as $\left[\begin{array}{c}{\dot{q}}_1\\ {\dot{\omega }}_1\\ {\dot{q}}_2\\ {\dot{\omega }}_2\end{array}\right]=\left[\begin{array}{c}{\omega }_1\times q_1\\ (q_1\times h_1)\times q_1\\ {\omega }_2\times q_2\\ (q_2\times h_2)\times q_2\end{array}\right]=F{|}_{(q,\omega )},$ with $h_i=h_i(q,\omega )$ and $\left[\begin{array}{c}{h}_1(q,\omega )\\ h_2(q,\omega )\end{array}\right]=R^{-1}\left(q\right)\left[\begin{array}{c}(-m_2{L}_1{L}_2|{\omega }_2{|}^2{\hat{q}}_2+(m_1+m_2\left)g{L}_1{\hat{e}}_3\right)q_1\\ (-m_2{L}_1{L}_2|{\omega }_1{|}^2{\hat{q}}_1+m_{2}g{L}_2{\hat{e}}_3)q_2\end{array}\right].$ Therefore, we can express the whole vector field in terms of the infinitesimal generator of the action of $SE\left(3\right)\times SE\left(3\right)$ as ${\psi }_{\ast }\left(f\right(q,\omega \left)\right){|}_{(q,\omega )}=F{|}_{(q,\omega )}$ through the function $f:T{S}^2\times T{S}^2\to \mathfrak{s}\mathfrak{e}\left(3\right)\times \mathfrak{s}\mathfrak{e}\left(3\right)\simeq {\mathbb{R}}^{12},(q,\omega )\to ({\omega }_1,q_1\times h_1,{\omega }_2,q_2\times h_2).$

5.3. Numerical experiments

In this section, we present some numerical experiment for the N-chain of pendulums. We start by comparing the various Lie group integrators that we have tested (with the choice N = 2), and conclude by analysing an implementation of variable step size. Lie group integrators allow to keep the evolution of the solution in the correct manifold, which is $T{S}^2\times T{S}^2$ when N = 2. Hence, we briefly report two sets of numerical experiments. In the first one, we show the convergence rate of all the Lie group integrators tested on this model. In the second one, we check how they behave in terms of preserving the two following relations:

• $q_i\left(t{)}^T{q}_i\right(t)=1,\mathrm{i}.\mathrm{e}.q_i(t)\in S^2,i=1,2,$

• $q_i\left(t{)}^T{\omega }_i\right(t)=0,\mathrm{i}.\mathrm{e}.{\omega }_i(t)\in T_{q_i\left(t\right)}{S}^2,i=1,2,$

completing the analysis with a comparison with the classical Runge–Kutta 4 and with ODE45 of MATLAB. The Lie group integrators used to obtain the following experiments are Lie Euler, Lie Euler Heun, three versions of Runge–Kutta–Munthe–Kaas methods of order 4 and one of order 3. The RKMK4 with two commutators mentioned in the plots is the one presented in Section 2, while the other schemes can be found for example in [7].

Figure 5 presents the plots of the errors, in logarithmic scale, obtained considering as a reference solution the one given by the ODE45 method, with strict tolerance. Here, we used an exact expression for the ${\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_{\sigma }^{-1}$ function. However, we could obtain the same results with a truncated version of this function, keeping a sufficiently high number of commutators, or after some clever manipulations of the commutators (as with RKMK4 with 2 commutators, see Section 2.2). The schemes show the right convergence rates, so we can move to the analysis of the time evolution on $T{S}^2\times T{S}^2$.

Figure 5. Convergence rate of the implemented Lie group integrators, based on global error considering as a reference solution the one of ODE45, with strict tolerance.

In Figure 6, we can see the comparison of the time evolution of the 2-norms of $q_1\left(t\right)$ and $q_2\left(t\right)$, for $0\le t\le T=5$. As highlighted above, unlike classical numerical integrators like the one implemented in ODE45 or the Runge–Kutta 4, the Lie group methods preserve the norm of the base components of the solutions, i.e. $|q_1\left(t\right)|=|q_2\left(t\right)|=1$ $\mathrm{\forall }t\in [0,T]$. Therefore, as expected, these integrators preserve the configuration manifold. However, to complete this analysis, we show the plots making a similar comparison but with the tangentiality conditions. Indeed, in Figure 7 we compare the time evolutions of the inner products $q_1\left(t{)}^T{\omega }_1\right(t)$ and $q_2\left(t{)}^T{\omega }_2\right(t)$ for $t\in [0,5]$, i.e. we see if these integrators preserve the geometry of the whole phase space $T{S}^2\times T{S}^2$. As we can see, while for Lie group methods these inner products are of the order of ${10}^{-14}$ and ${10}^{-15}$, the ones obtained with classical integrators show that the tangentiality conditions are not preserved with the same accuracy.

Figure 6. Visualization of the quantity $1-q_i\left(t{)}^T{q}_i\right(t)$, i = 1, 2, for time $t\in [0,5]$. These plots focus on the preservation of the geometry of $S^2$.

We now move to some experiments on variable stepsize. In this last part, we focus on the RKMK pair coming from Dormand–Prince method (DOPRI 5(4) [14]), which we denote with RKMK(5,4). The aim of the plots we show is to compare the same schemes, both with constant and variable stepsize. We start by setting a tolerance and solving the system with the RKMK(5,4) scheme. 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. Since we do not have an available analytical solution to quantify these two schemes' accuracy, we compare them with the solution obtained with a strict tolerance and ODE45. We compute such accuracy, at time T = 3, by means of the Euclidean norm of the ambient space ${\mathbb{R}}^{6N}$.

Figure 7. Visualization of the inner product $q_i\left(t{)}^T{\omega }_i\right(t)$, i = 1, 2, for $t\in [0,5]$. These plots focus on the preservation of the geometry of $T_{q_i\left(t\right)}{S}^2$.

In Figure 8, we compare the performance of the constant and variable stepsize methods, where the structure of the initial condition is always the same, but what changes is the number of connected pendulums. The considered initial condition is $(q_i,{\omega }_i)=(\sqrt{2}/2,0,\sqrt{2}/2,0,1,0),\mathrm{\forall }i=1,\dots ,N$, and all the masses and lengths are set to 1. From these experiments, we can notice situations where the variable step size beats the constant one in terms of accuracy at the final time, like the case N = 2 which we discuss in more detail afterwards.

Figure 8. Comparison of accuracy at final time (on the left) and step adaptation for the case N = 20 (on the right), with all pendulums of length $L_i=1$.

The results presented in Figure 10 (left) do not aim to highlight any particular relation between how the number of pendulums increases or the regularity of the solution. Indeed, as we add more pendulums, we keep incrementing the total length of the chain since $\sum _{i=1}^N{L}_i=N$. Thus here we do not have any appropriate limiting behaviour in the solution as $N\to +\mathrm{\infty }$. The behaviour presented in that figure seems to highlight an improvement in accuracy for the RKMK5 method as N increases. However, this is biased by the fact that when we increase N, to achieve the fixed tolerance of ${10}^{-6}$ with RKMKK(5,4), we need more time steps in the discretization. Thus, this plot does not say that as N increases, the dynamics becomes more regular; it suggests that the number of required timesteps increases faster than the ‘degree of complexity’ of the dynamics.

Figure 9. In these plots, we represent the six components of the solution describing the dynamics of the first mass (on the left) and of the second mass (on the right), for the case N = 2. We compare the behaviour of the solution obtained with constant stepsize RKMK5, the variable stepsize RKMK(5,4) and ODE45. (a) $\left(q_1\right(t),{\omega }_1(t\left)\right)$ (b) $\left(q_2\right(t),{\omega }_2(t\left)\right)$.

Figure 10. Comparison of accuracy at final time (on the left) and step adaptation for the case N = 20 (on the right), with all pendulums of length $L_i=5/N$.

For the case N = 2, we notice a relevant improvement passing to variable stepsize. In Figures 9 and 11 we can see that, for this choice of the parameters, the solution behaves smoothly in most of the time interval, but then there is a peak in the second component of the angular velocities of both the masses, at $t\approx 2.2$. We can observe this behaviour both in the plots of Figure 9, where we project the solution on the 12 components and even in Figure 11(c). In the latter, we plot two of the vector field components, i.e. the second components of the angular accelerations ${\dot{\omega }}_i\left(t\right)$, i = 1, 2. They show an abrupt change in the vector field in correspondence to $t\approx 2.2$, where the step is considerably restricted. Thus, to summarize, the gain we see with variable stepsize when N = 2 is motivated by the unbalance in the length of the time intervals with no abrupt changes in the dynamics and those where they appear. Indeed, we see that apart from a neighbourhood of $t\approx 2.2$, the vector field does not change quickly. On the other hand, for the case N = 20, this is the case. Thus the adaptivity of the stepsize does not bring relevant improvements in the latter situation.

Figure 11. On the left, we compare the adaptation of the stepsize of RKMK(5,4) with the one of ODE45 and with the constant stepsize of RKMK5. In the centre, we plot the second component of the angular velocities ${\omega }_i^{\left(2\right)}$, i = 1, 2, and we zoom in the last time interval $t\in [2.1,3]$ to see that the variable stepsize version of the method better reproduces the reference solution. On the right, we visualize the speed of variation of second component of the angular velocities. (a) Step adaptation, (b) Zoom at final times, (c) Values of ${\dot{\omega }}_i^{\left(2\right)}\left(t\right)$.

Figure 12. Two quadrotors connected to the mass point $m_y$ via massless links of lengths $L_i$.

The motivating application behind our choice of this mechanical system has been some intuitive relation with a beam model, as highlighted in the introduction of this work. However, for this limiting behaviour to make sense, we should fix the length of the entire chain of pendulums to some L (the length of the beam at rest) and then set the size of each pendulum to $L_i=L/N$. In this case, keeping the same tolerance of ${10}^{-6}$ for RKMK(5,4), we get the results presented in the following plot. We do not investigate more in detail this approach, which might be relevant for further work, however, we highlight that here the step adaptivity improves the results as we expected.

6. Dynamics of two quadrotors transporting a mass point

In this section, we consider a multibody system made of two cooperating quadrotor unmanned aerial vehicles (UAV) connected to a point mass (suspended load) via rigid links. This model is described in [28,29].

We consider an inertial frame whose third axis goes in the direction of gravity, but opposite orientation, and we denote with $y\in {\mathbb{R}}^3$ the mass point and with $y_1,y_2\in {\mathbb{R}}^3$ the two quadrotors. We assume that the links between the two quadrotors and the mass point are of a fixed length $L_1,L_2\in {\mathbb{R}}^{+}$. The configuration variables of the system are: the position of the mass point in the inertial frame, $y\in {\mathbb{R}}^3$, the attitude matrices of the two quadrotors, $(R_1,R_2)\in \left(SO\right(3){)}^2$ and the directions of the links which connect the centre of mass of each quadrotor respectively with the mass point, $(q_1,q_2)\in (S^2{)}^2$. The configuration manifold of the system is $Q={\mathbb{R}}^3\times \left(SO\right(3){)}^2\times (S^2{)}^2$.

In order to present the equations of motion of the system, we start by identifying $TSO\left(3\right)\simeq SO\left(3\right)\times \mathfrak{s}\mathfrak{o}\left(3\right)$ via left trivialization. This choice allows us to write the kinematic equations of the system as (36) ${\dot{R}}_i=R_i{\hat{\mathrm{\Omega }}}_i,{\dot{q}}_i={\hat{\omega }}_i{q}_{i}i=1,2,$(36) where ${\mathrm{\Omega }}_1,{\mathrm{\Omega }}_2\in {\mathbb{R}}^3$ represent the angular velocities of each quadrotor, respectively, and ${\omega }_1,{\omega }_2$ express the time derivatives of the orientations $q_1,q_2\in S^2$, respectively, in terms of angular velocities, expressed with respect to the body-fixed frames. From these equations, we define the trivialized Lagrangian $L(y,\dot{y},R_1,{\mathrm{\Omega }}_1,R_2,{\mathrm{\Omega }}_2,q_1,{\omega }_1,q_2,{\omega }_2):{\mathbb{R}}^6\times {\left(SO\left(3\right)\times \mathfrak{s}\mathfrak{o}\left(3\right)\right)}^2\times (T{S}^2{)}^2\to \mathbb{R},$ as the difference of the total kinetic energy of the system and the total potential (gravitational) energy, L = T−U, with: $T=\frac{1}{2}{m}_y\parallel \dot{y}{\parallel }^2+\frac{1}{2}\sum _{i=1}^2(m_i\parallel \dot{y}-L_i{\hat{\omega }}_i{q}_i{\parallel }^2+{\mathrm{\Omega }}_i^T{J}_i{\mathrm{\Omega }}_i),$ and $U=-m_{y}g{e}_3^{T}y-\sum _{i=1}^2{m}_{i}g{e}_3^T(y-L_i{q}_i),$ where $J_1,J_2\in {\mathbb{R}}^{3\times 3}$ are the inertia matrices of the two quadrotors and $m_1,m_2\in {\mathbb{R}}^{+}$ are their respective total masses. In this system, each of the two quadrotors generates a thrust force, which we denote with $u_i=-T_i{R}_i{e}_3\in {\mathbb{R}}^3$, where $T_i$ is the magnitude, while $e_3$ is the direction of this vector in the i th body-fixed frame, i = 1, 2. The presence of these forces make it a nonconservative system. Moreover, the rotors of the two quadrotors generate a moment vector, and we denote with $M_1,M_2\in {\mathbb{R}}^3$ the cumulative moment vector of each of the two quadrotors. To derive the Euler–Lagrange equations, a possible approach is through Lagrange–d'Alambert's principle, as presented in [28]. We write them in matrix form as (37) $A\left(z\right)\dot{z}=h\left(z\right)$(37) where $\begin{array}{rl}z& =[y,v,{\mathrm{\Omega }}_1,{\mathrm{\Omega }}_2,{\omega }_1,{\omega }_2{]}^T\in {\mathbb{R}}^{18},\\ A\left(z\right)& =\left[\begin{array}{cccccc}{I}_3& 0_3& 0_3& 0_3& 0_3& 0_3\\ 0_3& M_q& 0_3& 0_3& 0_3& 0_3\\ 0_3& 0_3& J_1& 0_3& 0_3& 0_3\\ 0_3& 0_3& 0_3& J_2& 0_3& 0_3\\ 0_3& -\frac{1}{L_1}{\hat{q}}_1& 0_3& 0_3& I_3& 0_3\\ 0_3& -\frac{1}{L_2}{\hat{q}}_2& 0_3& 0_3& 0_3& I_3\end{array}\right],\\ h\left(z\right)& =\left[\begin{array}{c}{h}_1\left(z\right)\\ h_2\left(z\right)\\ h_3\left(z\right)\\ h_4\left(z\right)\\ h_5\left(z\right)\\ h_6\left(z\right)\end{array}\right]=\left[\begin{array}{c}v\\ -\sum _{i=1}^2{m}_i{L}_i\parallel {\omega }_i{\parallel }^2{q}_i+M_{q}g{e}_3+\sum _{i=1}^2{u}_i^{\parallel }\\ -{\mathrm{\Omega }}_1\times J_1{\mathrm{\Omega }}_1+M_1\\ -{\mathrm{\Omega }}_2\times J_2{\mathrm{\Omega }}_2+M_2\\ -\frac{1}{L_1}g{\hat{q}}_1{e}_3-\frac{1}{m_1{L}_1}{q}_1\times u_1^{\perp }\\ -\frac{1}{L_2}g{\hat{q}}_2{e}_3-\frac{1}{m_2{L}_2}{q}_2\times u_2^{\perp }\end{array}\right],\end{array}$ where $M_q=m_y{I}_3+\sum _{i=1}^2{m}_i{q}_i{q}_i^T,$ and $u_i^{\parallel },u_i^{\perp }$ are respectively the orthogonal projection of $u_i$ along $q_i$ and to the plane $T_{q_i}{S}^2$, i = 1, 2, i.e. $u_i^{\parallel }=q_i{q}_i^T{u}_i$, $u_i^{\perp }=(I-q_i{q}_i^T)u_i$. These equations, coupled with the kinematic equations in (36(36) ${\dot{R}}_i=R_i{\hat{\mathrm{\Omega }}}_i,{\dot{q}}_i={\hat{\omega }}_i{q}_{i}i=1,2,$(36) ), describe the dynamics of a point $P=\left[y,v,R_1,{\mathrm{\Omega }}_1,R_2,{\mathrm{\Omega }}_2,q_1,{\omega }_1,q_2,{\omega }_2\right]\in M=TQ.$ Since the matrix $A\left(z\right)$ is invertible, we pass to the following set of equations: (38) $\dot{z}=A^{-1}\left(z\right)h\left(z\right):=\stackrel{\sim }{h}\left(z\right):=\overline{h}\left(P\right)=\left[{\overline{h}}_1\right(P),\dots ,{\overline{h}}_7(P){]}^T.$(38)

6.1. Analysis via transitive group actions

We identify the phase space M with $M\simeq T{\mathbb{R}}^3\times \left(TSO\right(3){)}^2\times (T{S}^2{)}^2$. The group we consider is $\overline{G}={\mathbb{R}}^6\times \left(TSO\right(3){)}^2\times (SE\left(3\right){)}^2,$ where the groups are combined with a direct-product structure and ${\mathbb{R}}^6$ is the additive group. For a group element $g=\left(\right(a_1,a_2),((B_1,b_1),(B_2,b_2)),((C_1,c_1),(C_2,c_2)\left)\right)\in \overline{G}$ and a point $P\in M$ in the manifold, we consider the following left action: $\begin{array}{rl}{\psi }_g\left(P\right)& =[y+a_1,v+a_2,B_1{R}_1,{\mathrm{\Omega }}_1+b_1,B_2{R}_2,{\mathrm{\Omega }}_2+b_2,\\ & C_1{q}_1,C_1{\omega }_1+c_1\times C_1{q}_1,C_2{q}_2,C_2{\omega }_2+c_2\times C_2{q}_2].\end{array}$ The well-definiteness and transitivity of this action come from standard arguments, see, e.g. [43]. The infinitesimal generator associated to $\xi =\left[{\xi }_1,{\xi }_2,{\eta }_1,{\eta }_2,{\eta }_3,{\eta }_4,{\mu }_1,{\mu }_2,{\mu }_3,{\mu }_4\right]\in \overline{\mathfrak{g}},$ where $\overline{\mathfrak{g}}=T_e\overline{G}$, writes $\begin{array}{rl}{\psi }_{\ast }\left(\xi \right){|}_P& =[{\xi }_1,{\xi }_2,{\hat{\eta }}_1{R}_1,{\eta }_2,{\hat{\eta }}_3{R}_2,{\eta }_4,\\ & {\hat{\mu }}_1{q}_1,{\hat{\mu }}_1{\omega }_1+{\hat{\mu }}_2{q}_1,{\hat{\mu }}_3{q}_2,{\hat{\mu }}_3{\omega }_2+{\hat{\mu }}_4{q}_2].\end{array}$ We can now focus on the construction of the function $f:M\to \overline{\mathfrak{g}}$ such that ${\psi }_{\ast }\left(f\right(P\left)\right){|}_P=F{|}_P$, where $\begin{array}{rl}F{|}_P& =\left[{\overline{h}}_1\right(P),{\overline{h}}_2(P),R_1{\hat{\mathrm{\Omega }}}_1,{\overline{h}}_3(P),R_2{\hat{\mathrm{\Omega }}}_2,\\ & {\overline{h}}_4\left(P\right),{\hat{\omega }}_1{q}_1,{\overline{h}}_5\left(P\right),{\hat{\omega }}_2{q}_2,{\overline{h}}_6\left(P\right)]\in T_{P}M\end{array}$ is the vector field obtained combining Equations (36(36) ${\dot{R}}_i=R_i{\hat{\mathrm{\Omega }}}_i,{\dot{q}}_i={\hat{\omega }}_i{q}_{i}i=1,2,$(36) ) and (38(38) $\dot{z}=A^{-1}\left(z\right)h\left(z\right):=\stackrel{\sim }{h}\left(z\right):=\overline{h}\left(P\right)=\left[{\overline{h}}_1\right(P),\dots ,{\overline{h}}_7(P){]}^T.$(38) ). We have $\begin{array}{rl}f\left(P\right)& =\left[{\overline{h}}_1\right(P),{\overline{h}}_2(P),R_1{\mathrm{\Omega }}_1,{\overline{h}}_3(P),R_2{\mathrm{\Omega }}_2,{\overline{h}}_4(P),\\ & {\omega }_1,q_1\times {\overline{h}}_5\left(P\right),{\omega }_2,q_2\times {\overline{h}}_6\left(P\right)]\in \overline{\mathfrak{g}}.\end{array}$ We have obtained the local representation of the vector field $F\in \mathfrak{X}\left(M\right)$ in terms of the infinitesimal generator of the transitive group action ψ, hence we can solve for one time step $\mathrm{\Delta }t$ the IVP $\left\{\begin{array}{l}\dot{\sigma }\left(t\right)={\mathrm{d}\mathrm{e}\mathrm{x}\mathrm{p}}_{\sigma \left(t\right)}^{-1}(f(\psi (\exp (\sigma \left(t\right)),P(t\left)\right)))\\ \sigma \left(0\right)=0\in \overline{\mathfrak{g}}\end{array}\right.$ and then update the solution $P(t+\mathrm{\Delta }t)=\psi (\exp (\sigma \left(\mathrm{\Delta }t\right)),P(t\left)\right)$.

The above construction is completely independent of the control functions $\{u_i^{\parallel },u_i^{\perp },M_i{\}}_{i=1,2}$ and hence it is compatible with any choice of these parameters.

6.2. Numerical experiments

We tested Lie group numerical integrators for a load transportation problem presented in [29]. The control inputs $\{u_i^{\parallel },u_i^{\perp },M_i{\}}_{i=1,2}$ are constructed such that the point mass asymptotically follows a given desired trajectory $y_d\in {\mathbb{R}}^3$, given by a smooth function of time, and the quadrotors maintain a prescribed formation relative to the point mass. In particular, the parallel components $u_i^{\parallel }$ are designed such that the payload follows the desired trajectory $y_d$ (load transportation problem), while the normal components $u_i^{\perp }$ are designed such that $q_i$ converge to desired directions $q_{id}$ (tracking problem in $S_2$). 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 neglect the construction of the controllers $M_i$ for the attitude dynamics of the quadrotors. However, the full dynamics model can also be easily integrated, once the expressions for the attitude controllers are available.

In Figure 13, we show the convergence rate of four different RKMK methods compared with the reference solution obtained with ODE45 in MATLAB.

Figure 13. Convergence rate of the numerical schemes compared with ODE45.

In Figures 14–18, we show results in the tracking of a parabolic trajectory, obtained by integrating the system (37(37) $A\left(z\right)\dot{z}=h\left(z\right)$(37) ) with a RKMK method of order 4.

Figure 14. Snapshots at $0\le t\le 5$.

Figure 15. Components of the load position (in blue) and the desired trajectory (in red) as a function time.

Figure 16. Deviation of the load position from the target trajectory.

Figure 17. Direction error of the links.

Figure 18. Preservation of the norms of $q_1,q_2\in S^2$.

7. Summary and outlook

In this paper, we have considered Lie group integrators with a particular focus on problems from mechanics. In mathematical terms, this means that the Lie groups and manifolds of particular interest are $SO\left(n\right),n=2,3$, $SE\left(n\right),n=2,3$ as well as the manifolds $S^2$ and $T{S}^2$. The abstract formulations by, e.g. Crouch and Grossman [11], Munthe-Kaas [41] and Celledoni et al. [6] have often been demonstrated on small toy problems in the literature, such as the free rigid body or the heavy top systems. But in papers like [4], hybrid versions of Lie group integrators have been applied to more complex beam and multi-body problems. The present paper is attempting to move in the direction of more relevant examples without causing the numerical solution to depend on how the manifold is embedded in an ambient space, or the choice of local coordinates.

It will be the subject of future work to explore more examples and to aim for a more systematic approach to applying Lie group integrators to mechanical problems. In particular, it is of interest to the authors to consider models of beams that could be seen as a generalization of the N-fold pendulum discussed here.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by Marie Sklodowska-Curie [860124].

Notes

1 If the Lie group action is smooth, a map f of the same regularity as F can be found [53].

2 ${\omega }_{(g,\mu )}$ is obtained from the natural symplectic form on $T^{\ast }G$ (which is a differential two-form), defined as ${\mathrm{\Omega }}_{(g,p_g)}\left(\right(\delta v_1,\delta {\pi }_1),(\delta v_2,\delta {\pi }_2\left)\right)=〈\delta {\pi }_2,\delta v_1〉-〈\delta {\pi }_1,\delta v_2〉,$ by right trivialization.

3 The Euler–Poincaré equations are Euler–Lagrange equations with respect to a Lagrange–d'Alembert principle obtained taking constraint variations.

4 In this paper, we will assume $\tilde{p}<p$ in which case the local error estimate is relevant for the approximation ${\tilde{y}}_{n+1}$.

5 There are many options for how to do this in practice, and the choice may also depend on the application. E.g. a Riemannian metric is a natural and robust alternative here.

6 It follows from the definition of the inertia tensor, i.e. $0\le \stackrel{\sim }{T}(q,\dot{q})=\frac{1}{2}\sum _{i=1}^N\left(\sum _{j\ge i}{m}_j\right)L_i{L}_j{\dot{q}}_i^T{\dot{q}}_j:=\frac{1}{2}{\dot{q}}^{T}M\dot{q}.$ Moreover, in this situation it is even possible to explicitly find the Cholesky factorization of the matrix M with an iterative algorithm.