# Lie Group integrators for mechanical systems

Elena Celledoni  Ergys Çokaj  Andrea Leone  Davide Murari  Brynjulf Owren

(Revised version July 2021, last compiled: )

###### 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(3)$
to $TS^{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 modeled on
$\mathbb{R}^{6}\times\left(SO(3)\times\mathfrak{so}(3)\right)^{2}\times(TS^{2})^{2}$
and put in a format where Lie group integrators can be applied.

##  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 [11, 40, 6]. But many numerical methods have been suggested that
add structure or utilise 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 [31, 13], 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 [49] and the symplectic and energy-momentum
methods by Lewis and Simo [31]. 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 [51, 37]
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 [27, 20, 29]. The majority of Lie
group methods found in the literature are one-step type generalisations for
classical methods, such as Runge–Kutta type formulas. In mechanical
engineering, the classical BDF methods have played an important role, and were
recently generalised [54] to Lie groups. Similarly, the celebrated
$\alpha$-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 generalised to a Lie group integrator.

The literature on Lie group integrators is rich and diverse, the interested
reader may consult the surveys [26, 10, 7, 44] 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 [39, 40] 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 the 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
succesfully applied for the simulation of mechanical systems, and in problems
of control, bio-mechanics and other engineering applications, see for example
[46], [27] [9], [25]. The present work is motivated by applications in
modeling and simulation of slender structures like Cosserat rods and beams
[49], 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(t)\in\mathcal{M}$ whose tangent at any point coincides with a vector
field $F\in\mathcal{X}(\mathcal{M})$ and passing through a designated initial
value $y_{0}$ at $t=t_{0}$

$\dot{y}(t)=F|_{y(t)},\qquad y(t_{0})=y_{0}.$ (1)

Before addressing numerical methods for solving (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}\rightarrow\mathcal{M}$ by some Lie
group $G$ of dimension $d\geq\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}\rightarrow G$ the exponential
map. We define $\psi_{*}:\mathfrak{g}\rightarrow\mathcal{X}(\mathcal{M})$ to
be the infinitesimal generator of the action, i.e.

$\left.F_{\xi}\right|_{m}=\left.\psi_{*}(\xi)\right|_{m}=\left.\frac{d}{dt}\right|_{t=0}\psi(\exp(t\xi),m)$
(2)

The transitivity of the action now ensures that
$\left.\psi_{*}(\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}=\left.\psi_{*}(\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}(\mathcal{M})$ there exists a map
$f:\mathcal{M}\rightarrow\mathfrak{g}$111If the Lie group action is smooth, a
map $f$ of the same regularity as $F$ can be found [53] such that

$F|_{m}=\left.\psi_{*}(f(m))\right|_{m},\quad\text{for all}\;m\in\mathcal{M}$
(3)

This is the original tool [40] 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},\ldots,E_{d}$ in $\mathcal{X}(\mathcal{M})$ was
introduced assuming that for every $m\in\mathcal{M}$,

$\text{span}\\{\left.E_{1}\right|_{m},\ldots,\left.E_{d}\right|_{m}\\}=T_{m}\mathcal{M}.$

Then, for any vector field $F\in\mathcal{X}(\mathcal{M})$ there are, in
general non-unique, functions $f_{i}:\mathcal{M}\rightarrow\mathbb{R}$, which
can be chosen with the same regularity as $F$, such that

$F|_{m}=\sum_{i=1}^{d}f_{i}(m)\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)

$\left.F_{\xi}\right|_{m}=\sum_{i=1}^{d}\xi_{i}E_{i}|_{m}$ (4)

If $\xi_{i}=f_{i}(p)$ 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 [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},\ldots,e_{d}$ is a basis
of the Lie algebra $\mathfrak{g}$, then we can simply define frame vector
fields as $E_{i}=\psi_{*}(e_{i})$ and the vector field we aim to describe is,

$F|_{m}=\left.\psi_{*}(f(m))\right|_{m}=\left.\psi_{*}(\sum_{i}f_{i}(m)e_{i})\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}(\mathcal{M})$ satisfying

  1. 1.

$F_{p}|_{p}=F|_{p}$

  2. 2.

$\nabla_{X}F_{p}=0$ for any $X\in\mathcal{X}(M)$.

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 $\nabla$.
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. [32].

######  Example 1.

For mechanical systems on Lie groups, two important constructions are the
adjoint and coadjoint representatons. For every $g\in G$ there is an
automorphism $\mathrm{Ad}_{g}:\mathfrak{g}\rightarrow\mathfrak{g}$ defined as

$\mathrm{Ad}_{g}(\xi)=TL_{g}\circ TR_{g^{-1}}(\xi)$

where $L_{g}$ and $R_{g}$ are the left and right multiplications respectively,
$L_{g}(h)=gh$ and $R_{g}(h)=hg$. Since $\mathrm{Ad}$ is a representation, i.e.
$\mathrm{Ad}_{gh}=\mathrm{Ad}_{g}\circ\mathrm{Ad}_{h}$ it also defines a left
Lie group action by $G$ on $\mathfrak{g}$. From this definition and a duality
pairing $\langle\cdot,\cdot\rangle$ between $\mathfrak{g}$ and
$\mathfrak{g}^{*}$, we can also derive a representation on $\mathfrak{g}^{*}$
denoted $\mathrm{Ad}_{g}^{*}$, simply by

$\langle\mathrm{Ad}_{g}^{*}(\mu),\xi\rangle=\langle\mu,\mathrm{Ad}_{g}(\xi)\rangle,\quad\xi\in\mathfrak{g},\
\mu\in\mathfrak{g}^{*}.$

The action $g\cdot\mu=\mathrm{Ad}_{g^{-1}}^{*}(\mu)$ has infinitesimal
generator given as

$\left.\psi_{*}(\xi)\right|_{\mu}=-\mathrm{ad}_{\xi}^{*}\mu$

Following [34], for a Hamiltonian $H:T^{*}G\rightarrow\mathbb{R}$, define
$H^{-}$ to be its restriction to $\mathfrak{g}^{*}$. Then the Lie-Poisson
reduction of the dynamical system is defined on $\mathfrak{g}^{*}$ as

$\dot{\mu}=-\mathrm{ad}^{*}_{\frac{\partial H^{-}}{\partial\mu}}\mu$

and this vector field is precisely of the form (3) with $f(\mu)=\frac{\partial
H^{-}}{\partial\mu}(\mu)$. 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(3)$. The
Lie algebra $\mathfrak{so}(3)$ can be modeled 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

$\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{so}(3)^{*}$ as vectors in
$\mathbb{R}^{3}$ with duality pairing $\langle\mu,\xi\rangle=\mu^{T}\xi$. With
these representations, we find that the coadjoint action can be expressed as

$g\cdot\mu=\psi(g,\mu)=\mathrm{Ad}_{g^{-1}}^{*}\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)^{*}$ with a
left invariant Hamiltonian which reduces to the function

$H^{-}(\mu)=\frac{1}{2}\langle\mu,\mathbb{I}^{-1}\mu\rangle$

on $\mathfrak{so}(3)^{*}$ where
$\mathbb{I}:\mathfrak{so}(3)\rightarrow\mathfrak{so}(3)^{*}$ is the inertia
tensor. From this, we can now set $f(\mu)=\partial
H^{-}/\partial\mu=\mathbb{I}^{-1}\mu$. We then recover the Euler free rigid
body equation as

$\dot{\mu}=\left.\psi_{*}(f(\mu)\right|_{\mu}=-\mathrm{ad}_{\mathbb{I}^{-1}\mu}^{*}\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(y)$, $y(0)=y_{0}$ 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 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(t)=\exp(\sigma(t))\cdot y_{0},\quad\sigma(0)=0.$

The transformed differential equation for $\sigma(t)$ makes use of the
derivative of the exponential mapping, the reader should consult [40] for
details about the derivation, we give the final result

$\dot{\sigma}(t)=\mathrm{dexp}_{\sigma(t)}^{-1}(f(\exp(\sigma(t))\cdot
y_{0}))$ (6)

The map $v\mapsto\mathrm{dexp}_{u}(v)$ is linear and invertible when $u$
belongs to some sufficiently small neighborhood of $0\in\mathfrak{g}$. It has
an expansion in nested Lie brackets [21]. Using the operator
$\mathrm{ad}_{u}(v)=[u,v]$ and its powers $\mathrm{ad}_{u}^{2}v=[u,[u,v]]$
etc, one can write

$\mathrm{dexp}_{u}(v)=\left.\frac{e^{z}-1}{z}\right|_{z=\mathrm{ad}_{u}}(v)=v+\frac{1}{2}[u,v]+\frac{1}{6}[u,[u,v]]+\cdots$
(7)

and the inverse is

$\mathrm{dexp}_{u}^{-1}(v)=\left.\frac{z}{e^{z}-1}\right|_{z=\mathrm{ad}_{u}}(v)=v-\frac{1}{2}[u,v]+\frac{1}{12}[u,[u,v]]+\cdots$
(8)

The RKMK methods are now obtained simply by applying some standard Runge–Kutta
method to the transformed equation (6) with a time step $h$, using initial
value $\sigma(0)=0$. This leads to an output $\sigma_{1}\in\mathfrak{g}$ and
one simply sets $y_{1}=\exp(\sigma_{1})\cdot y_{0}$. Then one repeats the
procedure replacing $y_{0}$ by $y_{1}$ in the next step etc. While solving (6)
one needs to evaluate $\mathrm{dexp}_{u}^{-1}(v)$ as a part of the process.
This can be done by truncating the series (8) since $\sigma(0)=0$ implies that
we always evaluate $\mathrm{dexp}_{u}^{-1}$ with $u=\mathcal{O}(h)$, and thus,
the $k$th iterated commutator $\mathrm{ad}_{u}^{k}=\mathcal{O}(h^{k})$. For a
given Runge–Kutta method, there are some clever tricks that can be done to
minimise the total number of commutators to be included from the expansion of
$\mathrm{dexp}_{u}^{-1}v$, see [5, 41]. We give here one concrete example of
an RKMK method proposed in [5]

$\displaystyle f_{n,1}$ $\displaystyle=hf(y_{n}),$ $\displaystyle f_{n,2}$
$\displaystyle=hf(\exp(\tfrac{1}{2}f_{n,1})\cdot y_{n}),$ $\displaystyle
f_{n,3}$
$\displaystyle=hf(\exp(\tfrac{1}{2}f_{n,2}-\tfrac{1}{8}[f_{n,1},f_{n,2}])\cdot
y_{n}),$ $\displaystyle f_{n,4}$ $\displaystyle=hf(\exp(f_{n,3})\cdot y_{n}),$
$\displaystyle y_{n+1}$
$\displaystyle=\exp(\tfrac{1}{6}(f_{n,1}+2f_{n,2}+2f_{n,3}+f_{n,4}-\tfrac{1}{2}[f_{n,1},f_{n,4}]))\cdot
y_{n}.$

The other option is to compute the exact expression for
$\mathrm{dexp}_{u}^{-1}(v)$ for the particular Lie algebra we use. For
instance, it was shown in [8] that for the Lie algebra $\mathfrak{so}(3)$ one
has

$\mathrm{dexp}_{u}^{-1}(v)=v-\frac{1}{2}u\times
v+\alpha^{-2}(1-\tfrac{\alpha}{2}\cot\tfrac{\alpha}{2})\;u\times(u\times v)$

We will present the corresponding formula for $\mathfrak{se}(3)$ 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

$\displaystyle Y_{n,r}$
$\displaystyle=\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}$ $\displaystyle f_{n,r}$ $\displaystyle=f(Y_{n,r})$ $\displaystyle
y_{n+1}$
$\displaystyle=\exp\left(h\sum_{k}\beta_{J}^{k}f_{n,k}\right)\cdots\exp\left(h\sum_{k}\beta_{1}^{k}f_{n,k}\right)\cdot
y_{n}$

Here the Runge–Kutta coefficients $\alpha_{r,j}^{k}$, $\beta_{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}\beta_{j}^{r}$. The $\alpha_{r,j}^{k}$, $\beta_{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
generalisation of the classical 4th order Runge–Kutta method.

$\displaystyle\begin{split}Y_{n,1}&=y_{n}\\\
Y_{n,2}&=\exp(\tfrac{1}{2}hf_{n,1})\cdot y_{n}\\\
Y_{n,3}&=\exp(\tfrac{1}{2}hf_{n,2})\cdot y_{n}\\\
Y_{n,4}&=\exp(hf_{n,3}-\tfrac{1}{2}hf_{n,1})\cdot Y_{n,2}\\\
y_{n+\frac{1}{2}}&=\exp(\tfrac{1}{12}h(3f_{n,1}+2f_{n,2}+2f_{n,3}-f_{n,4}))\cdot
y_{n}\\\
y_{n+1}&=\exp(\tfrac{1}{12}h(-f_{n,1}+2f_{n,2}+2f_{n,3}+3f_{n,4}))\cdot
y_{n+\frac{1}{2}}\end{split}$ (9)

where $f_{n,i}=f(Y_{n,i})$. 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{dexp}_{u}^{-1}(v)$ in
$\mathfrak{se}(3)$

As an alternative to using a truncated version of the infinite series for
$\mathrm{dexp}_{u}^{-1}$ (8), one can consider exact expressions obtained for
certain Lie algebras. Since $\mathfrak{se}(3)$ is particularly important in
applications to mechanics, we give here its exact expression. For this, we
represent elements of $\mathfrak{se}(3)$ 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) and $a$
is the translational part. Now, let $\varphi(z)$ be a real analytic function
at $z=0$. We define

$\varphi_{+}(z)=\frac{\varphi(iz)+\varphi(-iz)}{2},\qquad\varphi_{-}(z)=\frac{\varphi(iz)-\varphi(-iz)}{2i}$

We next define the four functions

$g_{1}(z)=\frac{\varphi_{-}(z)}{z},\
\tilde{g}_{1}(z)=\frac{g_{1}^{\prime}(z)}{z},\quad
g_{2}(z)=\frac{\varphi(0)-\varphi_{+}(z)}{z^{2}},\
\tilde{g}_{2}(z)=\frac{g_{2}^{\prime}(z)}{z}$

and the two scalars $\rho=A^{T}a$, $\alpha=\|A\|_{2}$. One can show that for
any $(A,a)$ and $(B,b)$ in $\mathfrak{se}(3)$, it holds that

$\varphi(\mathrm{ad}_{(A,a)})(B,b)=(C,c)$

where

$\displaystyle C$ $\displaystyle=\varphi(0)B+g_{1}(\alpha)A\times
B+g_{2}(\alpha)\,A\times(A\times B)$ $\displaystyle c$
$\displaystyle=\varphi(0)b+g_{1}(\alpha)\,(a\times B+A\times
b)+\rho\tilde{g}_{1}(\alpha)\,A\times
B+\rho\tilde{g}_{2}(\alpha)\,A\times(A\times B)$
$\displaystyle+g_{2}(\alpha)\,(a\times(A\times B)+A\times(a\times
B)+A\times(A\times b))$

Considering for instance (8), we may now use $\varphi(z)=\frac{z}{e^{z}-1}$ to
calculate

$g_{1}(z)=-\frac{1}{2},\ \tilde{g}_{1}(z)=0,\
g_{2}(z)=\frac{1-\tfrac{z}{2}\cot\tfrac{z}{2}}{z^{2}},\
\tilde{g}_{2}(z)=\frac{1}{z}\frac{d}{dz}g_{2}(z),\ \varphi(0)=1.$

and thereby obtain an expression for $\mathrm{dexp}_{(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{ad}$-operator for several other Lie groups
of small and medium dimension. For instance in [38] a variety of coordinate
mappings for rigid body motions are discussed. For Lie algebras of larger
dimension, both the exponential mapping and $\mathrm{dexp}_{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
$\phi:\mathfrak{g}\rightarrow 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 the various cases with different formulations of the heavy top
system.

###  3.1 Semi-direct products

The coadjoint action by $G$ on $\mathfrak{g}^{*}$ is denoted
$\mathrm{Ad}_{g}^{*}$ defined for any $g\in G$ as

$\langle\mathrm{Ad}_{g}^{*}\mu,\xi\rangle=\langle\mu,\mathrm{Ad}_{g}\xi\rangle,\quad\forall\xi\in\mathfrak{g},$
(10)

where $\mathrm{Ad}:\mathfrak{g}\rightarrow\mathfrak{g}$ is the adjoint
representation and for a duality pairing $\langle{\cdot},{\cdot}\rangle$
between $\mathfrak{g}^{*}$ and $\mathfrak{g}$.

We consider the cotangent bundle of a Lie group $G$, $T^{*}G$ and identify it
with $G\times\mathfrak{g}^{*}$ using the right multiplication
$R_{g}:G\rightarrow G$ and its tangent mapping $R_{g*}:=TR_{g}$. The cartesian
product $G\times\mathfrak{g}^{*}$ can be given a semi-direct product structure
that turns it into a Lie group $\mathbf{G}:=G\ltimes\mathfrak{g}^{*}$ where
the group multiplication is

$(g_{1},\mu_{1})\cdot(g_{2},\mu_{2})=(g_{1}\cdot
g_{2},\mu_{1}+\mathrm{Ad}^{*}_{g_{1}^{-1}}\mu_{2}).$ (11)

Acting by left multiplication any vector field $F\in\mathcal{X}(\mathbf{G})$
is expressed by means of a map $f\colon\mathbf{G}\rightarrow T_{e}\mathbf{G}$,

$F(g,\mu)=T_{e}R_{(g,\mu)}f(g,\mu)=(R_{g*}f_{1},f_{2}-\mathrm{ad}_{f_{1}}^{*}\mu),$
(12)

where $f_{1}=f_{1}(g,\mu)\in\mathfrak{g}$,
$f_{2}=f_{2}(g,\mu)\in\mathfrak{g}^{*}$ are the two components of $f$.

###  3.2 Symplectic form and Hamiltonian vector fields

The right trivialised222$\omega_{(g,\mu)}$ is obtained from the natural
symplectic form on $T^{*}G$ (which is a differential two-form), defined as
$\Omega_{(g,p_{g})}((\delta v_{1},\delta\pi_{1}),(\delta
v_{2},\delta\pi_{2}))=\langle\delta\pi_{2},\delta
v_{1}\rangle-\langle\delta\pi_{1},\delta v_{2}\rangle,$ by right
trivialization.  symplectic form pulled back to $\mathbf{G}$ reads

$\begin{split}\omega_{(g,\mu)}((R_{g*}\xi_{1},\delta\nu_{1}),(R_{g*}\xi_{2},\delta\nu_{2}))&=\langle\delta\nu_{2},\xi_{1}\rangle+\\\
&-\langle\delta\nu_{1},\xi_{2}\rangle-\langle\mu,[\xi_{1},\xi_{2}]\rangle,\qquad\xi_{1},\xi_{2}\in\mathfrak{g}.\end{split}$
(13)

See [31] 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^{*}G\rightarrow\mathbb{R}$, where
$\mathrm{i}_{F}$ is defined as $\mathrm{i}_{F}(X):=\omega(F,X)$ for any vector
field $X$. This implies that the map $f$ for such a Hamiltonian vector field
gets the form

$f(g,\mu)=\left(\frac{\partial H}{\partial\mu}(g,\mu),-R_{g}^{*}\frac{\partial
H}{\partial g}(g,\mu)\right).$ (14)

The following is a one-parameter family of symplectic Lie group integrators on
$T^{*}G$:

$\displaystyle M_{\theta}$
$\displaystyle=\mathrm{dexp}_{-\xi}^{*}(\mu_{0}+\mathrm{Ad}^{*}_{\exp(\theta\xi)}(\bar{n}))-\theta\mathrm{dexp}_{-\theta\xi}^{*}\mathrm{Ad}^{*}_{\exp(\theta\xi)}(\bar{n}),$
(15) $\displaystyle(\xi,\bar{n})$
$\displaystyle=hf\bigl{(}\exp(\theta\xi)\cdot g_{0},M_{\theta}\bigr{)},$ (16)
$\displaystyle(g_{1},\mu_{1})$
$\displaystyle=(\exp(\xi),\mathrm{Ad}_{\exp((\theta-1)\xi)}^{*}\bar{n})\cdot(g_{0},\mu_{0}).$
(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^{*}G$ with a left or
right invariant Hamiltonian can be reduced to a system on $\mathfrak{g}^{*}$
[33]. In fact for a Hamiltonian $H$ right invariant under the left action of
$G$, $\frac{\partial H}{\partial g}=0$, and from (12) and (14) we get for the
second equation

$\dot{\mu}=\mp\mathrm{ad}^{*}_{\frac{\partial H}{\partial\mu}}\mu,$ (18)

where the positive sign is used in case of left invariance (see e.g. section
13.4 in [35]). The solution to this system preserves coadjoint orbits, thus
using the Lie group action

$g\cdot\mu=\mathrm{Ad}_{g^{-1}}^{*}\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é equations333The Euler–Poincaré equations are Euler–Lagrange
equations with respect to a Lagrange–d’Alembert principle obtained taking
constraint variations., [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^{*}SO(3)$ where the
equations are canonical Hamiltonian, a second point of view is that the system
is a Lie–Poisson system on $\mathfrak{se}(3)^{*}$, 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.

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

####  3.4.1 Heavy top equations on $T^{*}SO(3)$.

The heavy top is a rigid body with a fixed point in a gravitational field. The
phase space of this mechanical system is $T^{*}SO(3)$ where the equations of
the heavy top are in canonical Hamiltonian form. Assuming $(Q,p)$ are
coordinates for $T^{*}SO(3)$, $\Pi=(T_{e}L_{Q})^{*}(p)$ is the left
trivialized or body momentum. The Hamiltonian of the heavy top is given in
terms of $(Q,\Pi)$ as

$H\colon SO(3)\ltimes\mathfrak{so}(3)^{*}\to\mathbb{R},\quad
H(Q,\Pi)=\frac{1}{2}\langle\Pi,\mathbb{I}^{-1}\Pi\rangle+Mg\ell\,\mathbf{\Gamma}\cdot\mathcal{X},\quad\mathbf{\Gamma}=Q^{-1}\mathbf{\Gamma}_{0},$

where $\mathbb{I}\colon\mathfrak{so}(3)\rightarrow\mathfrak{so}(3)^{*}$ 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 center of mass of the body, $\ell$ is the length of this segment.
The equations of motion on $SO(3)\ltimes\mathfrak{so}(3)^{*}$ are

$\displaystyle\dot{\Pi}$
$\displaystyle=\Pi\times\mathbb{I}^{-1}\Pi+Mg\ell\,\mathbf{\Gamma}\times\mathcal{X},$
(19) $\displaystyle\dot{Q}$ $\displaystyle=Q\,\widehat{\mathbb{I}^{-1}\Pi}.$
(20)

The identification of $T^{*}SO(3)$ with $SO(3)\ltimes\mathfrak{so}(3)^{*}$ via
right trivialization leads to the spatial momentum variable
$\pi=(T_{e}R_{Q})^{*}(p)=Q\Pi$. The equations written in the space variables
$(Q,\pi)$ get the form

$\displaystyle\dot{\pi}$ $\displaystyle=Mg\ell\,\mathbf{\Gamma}_{0}\times
Q\mathcal{X},$ (21) $\displaystyle\dot{Q}$
$\displaystyle=\hat{\omega}Q\qquad\omega=Q\mathbb{I}^{-1}Q^{T}\pi.$ (22)

where, the first equation states that the component of $\pi$ parallel to
$\mathbf{\Gamma}_{0}$ is constant in time. These equations can be obtained
from (12) and (14) on the right trivialized $T^{*}SO(3)$,
$SO(3)\ltimes\mathfrak{so}(3)^{*}$, with the heavy top Hamiltonian and the
symplectic Lie group integrators (16)-(17) can be applied in this case.
Similar methods were proposed in [31] and [48].

####  3.4.2 Heavy top equations on $\mathfrak{se}^{*}(3)$

The Hamiltonian of the heavy top is not invariant under the action of $SO(3)$,
so the equations (19)-(20) given in section (3.4.1) cannot be reduced to
$\mathfrak{so}^{*}(3)$, nevertheless the heavy top equations are Lie–Poisson
on $\mathfrak{se}^{*}(3)$, [52, 17, 47].

Observe that the equations of the heavy top on $T^{*}SO(3)$ (19)-(20) can be
easily modified eliminating the variable $Q\in SO(3)$ and replacing it with
$\mathbf{\Gamma}\in\mathbb{R}^{3}$ $\mathbf{\Gamma}=Q^{-1}\mathbf{\Gamma}_{0}$
to obtain

$\displaystyle\dot{\Pi}$
$\displaystyle=\Pi\times\mathbb{I}^{-1}\Pi+Mg\ell\,\mathbf{\Gamma}\times\mathcal{X},$
(23) $\displaystyle\dot{\mathbf{\Gamma}}$
$\displaystyle=\mathbf{\Gamma}\times(\mathbb{I}^{-1}\Pi).$ (24)

We will see that the solutions of these equations evolve on
$\mathfrak{se}^{*}(3)$. In what follows, we consider elements of
$\mathfrak{se}^{*}(3)$ to be pairs of vectors in $\mathbb{R}^{3}$, e.g.
$(\Pi,\mathbf{\Gamma})$. Correspondingly the elements of $SE(3)$ are
represented as pairs $(g,\mathbf{u})$ with $g\in SO(3)$ and
$\mathbf{u}\in\mathbb{R}^{3}$. The group multiplication in $SE(3)$ 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(3)$ 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{se}^{*}(3)$ take the form

$\mathrm{Ad}^{*}_{(g,\mathbf{u})}(\Pi,\mathbf{\Gamma})=(g^{-1}(\Pi-\mathbf{u}\times\mathbf{\Gamma}),g^{-1}\mathbf{\Gamma}),\;\;\mathrm{ad}^{*}_{(\xi,\mathbf{u})}(\Pi,\mathbf{\Gamma})=(-\xi\times\Pi-\mathbf{u}\times\mathbf{\Gamma},-\xi\times\mathbf{\Gamma}).$

Using this expression for $\mathrm{ad}^{*}_{(\xi,\mathbf{u})}$ with
$(\xi=\frac{\partial H}{\partial\Pi},\mathbf{u}=\frac{\partial
H}{\partial\mathbf{\Gamma}})$, it can be easily seen that the equations (18)
in this setting reproduce the heavy top equations (23)-(24). Therefore the
equations are Lie–Poisson equations on $\mathfrak{se}^{*}(3)$. 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^{*}SE(3)$.
Moreover the Hamiltonian on $\mathfrak{se}(3)^{*}$ 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 $\Pi$, $Q$, $\mathbf{p}$ and $\mathbf{q}$:

$\displaystyle H\colon
T^{*}SO(3)\times{\mathbb{R}^{3}}^{*}\times{\mathbb{R}^{3}}\to\mathbb{R},$
$\displaystyle
H((\Pi,Q),(\mathbf{p},\mathbf{q}))=\frac{1}{2}\langle\Pi,\mathbb{I}^{-1}\Pi\rangle+\frac{1}{2}\|\mathbf{p}-Q^{-1}\mathbf{\Gamma}_{0}\|^{2}-\frac{1}{2}\|Q^{-1}\mathbf{\Gamma}_{0}\|^{2},$

see [23, section 8.5]. This Hamiltonian is invariant under the left action of
$SO(3)$. The corresponding equations are canonical on $T^{*}S\equiv
S\ltimes\mathfrak{s}^{*}$ where $S=SO(3)\times\mathbb{R}^{3}$ with Lie algebra
$\mathfrak{s}:=\mathfrak{so}(3)\times{\mathbb{R}^{3}}$ and $T^{*}S$ can be
identified with $T^{*}SO(3)\times{\mathbb{R}^{3}}^{*}\times\mathbb{R}^{3}$.
The equations are

$\displaystyle\dot{\Pi}$
$\displaystyle=\Pi\times\mathbb{I}^{-1}\Pi-(Q^{-1}\mathbf{\Gamma}_{0})\times\mathbf{p},$
(25) $\displaystyle\dot{Q}$ $\displaystyle=Q\,\widehat{\mathbb{I}^{-1}\Pi},$
(26) $\displaystyle\dot{\mathbf{p}}$ $\displaystyle=\mathbf{0},$ (27)
$\displaystyle\dot{\mathbf{q}}$
$\displaystyle=\mathbf{p}-Q^{-1}\mathbf{\Gamma}_{0}.$ (28)

and in the spatial momentum variables

$\displaystyle\dot{\pi}$ $\displaystyle=-\mathbf{\Gamma}_{0}\times
Q\mathbf{p},$ (29) $\displaystyle\dot{Q}$
$\displaystyle=\hat{\omega}Q,\quad\omega=Q\mathbb{I}^{-1}Q^{T}\pi,$ (30)
$\displaystyle\dot{\mathbf{p}}$ $\displaystyle=\mathbf{0},$ (31)
$\displaystyle\dot{\mathbf{q}}$
$\displaystyle=\mathbf{p}-Q^{-1}\mathbf{\Gamma}_{0}.$ (32)

Similar formulations were considered in [30] 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(0)=I$, $\ell=2$,
$M=15$
$\mathbb{I}=\mathrm{diag}(0.234375,0.46875,0.234375)$,
$\pi(0)=\mathbb{I}(0,150,-4.61538)$, $\mathcal{X}=(0,1,0)$
$\Gamma_{0}=(0,0,-9.81)$.

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

In Figure 2 we report the performance of the symplectic Lie group integrators
(15)-(17) applied both on the equations (21)-(22) with $\theta=0$,
$\theta=\frac{1}{2}$ and $\theta=1$ (SLGI), and to the equations (29)-(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 3: Symplectic Lie group integrators, long time integration, $h=0.01$,
$6000$ steps.. Top: energy error, bottom 3D plot of $M\ell Q^{-1}\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 summarise 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 [36] and a detailed treatment can be found in [19, p.
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{tol}$ in
every step. Suppose the main method is of order $p$ and the auxiliary method
is of order $\tilde{p}\neq p.$ 444In this paper we will assume $\tilde{p}<p$
in which case the local error estimate is relevant for the approximation
$\tilde{y}_{n+1}$ 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 measure555There 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. 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}=Ch_{n+1}^{\tilde{p}+1}+\mathcal{O}(h^{\tilde{p}+2})$. Aiming at
$e_{n+1}\approx\mathrm{tol}$ in every step, one may use a formula of the type

$h_{n+1}=\theta\left(\frac{\mathrm{tol}}{e_{n+1}}\right)^{\tfrac{1}{\tilde{p}+1}}\,h_{n}$
(33)

where $\theta$ is a ‘safety factor’, typically chosen between $0.8$ and $0.9$.
In case the step is rejected because $e_{n}>\mathrm{tol}$ we can redo the step
with a step size obtained by the same formula. We summarise the approach in
the following algorithm

---
Given $y_{n}$, $h_{n}$, $\mathrm{tol}$
Let $h:=h_{n}$
repeat
| Compute $y_{n+1}$, $\tilde{y}_{n+1}$, $e_{n+1}$ from $y_{n}$, $h$
| Update stepsize
$h:=\theta\left(\frac{\mathrm{tol}}{e_{n+1}}\right)^{\alpha}h$
| accepted := $e_{n+1}<\mathrm{tol}$
|  if accepted
|  | update step index: $n:=n+1$
|  |  $h_{n}:=h$
until accepted

Here we have used again the safety factor $\theta$, and the parameter $\alpha$
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) 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) to
obtain local Lie algebra approximations $\sigma_{1}$, $\tilde{\sigma}_{1}$ and
one uses e.g. $e_{n+1}=\|\sigma_{1}-\tilde{\sigma}_{1}\|$ (note that the
equation itself depends on $y_{n}$). For methods which use a truncated version
of the series for $\mathrm{dexp}_{u}^{-1}$ one may also try to optimise
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 [43] 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

$\displaystyle Y_{n,1}$ $\displaystyle=y_{n}$ $\displaystyle Y_{n,2}$
$\displaystyle=\exp(\tfrac{1}{3}hf_{n,1})\cdot y_{n}$ $\displaystyle Y_{n,3}$
$\displaystyle=\exp(\tfrac{2}{3}hf_{n,2})\cdot y_{n}$ $\displaystyle y_{n+1}$
$\displaystyle=\exp(h(-\tfrac{1}{12}f_{n,1}+\tfrac{3}{4}f_{n,3}))\cdot
Y_{n,2}$ $\displaystyle\tilde{y}_{n+1}$
$\displaystyle=\exp(\tfrac{1}{2}h(f_{n,2}+f_{n,3}))\cdot y_{n}$

One could also have reused the third stage $Y_{n,3}$ in the update, rather
than $Y_{n,2}$.

$\displaystyle Y_{n,1}$ $\displaystyle=y_{n}$ $\displaystyle Y_{n,2}$
$\displaystyle=\exp(\tfrac{2}{3}hf_{n,1})\cdot y_{n}$ $\displaystyle Y_{n,3}$
$\displaystyle=\exp(h(\tfrac{5}{12}f_{n,1}+\tfrac{1}{4}f_{n,2})\cdot y_{n}$
$\displaystyle y_{n+1}$
$\displaystyle=\exp(h(-\tfrac{1}{6}f_{n,1}-\tfrac{1}{2}f_{n,2}+f_{n,3}))\cdot
Y_{n,3}$ $\displaystyle\tilde{y}_{n+1}$
$\displaystyle=\exp(\tfrac{1}{4}h(f_{n,1}+3f_{n,3}))\cdot y_{n}$

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 minimise 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. 1.

Compute $Y_{n,i},\;i=1\ldots,4$ and $y_{n+1}$ from (9)

  2. 2.

Compute an additional stage $\bar{Y}_{n,3}$ and then $\tilde{y}_{n+1}$ as

$\displaystyle\begin{split}\bar{Y}_{n,3}&=\exp(\tfrac{3}{4}hf_{n,2})\cdot
y_{n}\\\
\tilde{y}_{n+1}&=\exp(\tfrac{h}{9}(-f_{n,1}+3f_{n,2}+4\bar{f}_{n,3}))\cdot\exp(\tfrac{h}{3}f_{n,1})\cdot
y_{n}\end{split}$ (34)

##  5 The $N$-fold 3D pendulum

In this section, we present a model for a system of $N$ connected
3-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.

Figure 4: $3-$fold pendulum at a fixed time instant, with fixed point placed
at the origin.

The Lagrangian we consider is a function from $(TS^{2})^{N}$ to $\mathbb{R}$.
Instead of the coordinates $(q_{1},...,q_{N},\dot{q}_{1},...,\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\\}=\langle
q_{i}\rangle^{\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},...,q_{N},\omega_{N})\in(TS^{2})^{N}$ to
describe the dynamics. In this section we denote with $m_{1},...,m_{N}$ the
masses of the pendulums and with $L_{1},...,L_{N}$ their lengths. Figure 4
shows the case $N=3$. We organize the section into three parts:

  1. 1.

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

  2. 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. 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 $(TS^{2})^{N}$

We characterize a transitive action for $(TS^{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{se}(3)$, the Lie algebra of $SE(3)$, and
$\mathbb{R}^{6}$. We start from the Ad-action of $SE(3)$ on $\mathfrak{se}(3)$
(see [23]), which writes

$\mathrm{Ad}:SE(3)\times\mathfrak{se}(3)\rightarrow\mathfrak{se}(3),$
$\mathrm{Ad}((R,r),(u,v))=(Ru,Rv+\hat{r}Ru).$

Since $\mathfrak{se}(3)\simeq\mathbb{R}^{6}$, the Ad-action allows us to
define the following Lie group action on $\mathbb{R}^{6}$

$\psi:SE(3)\times\mathbb{R}^{6}\rightarrow\mathbb{R}^{6},\;\;\psi((R,r),(u,v))=(Ru,Rv+\hat{r}Ru).$

We can think of $\psi$ as a Lie group action on $TS^{2}$ since, for any
$q\in\mathbb{R}^{3}$, it maps points of

$TS_{|q|}^{2}:=\\{(\tilde{q},\tilde{\omega})\in\mathbb{R}^{3}\times\mathbb{R}^{3}:\;\tilde{\omega}^{T}\tilde{q}=0,\;|\tilde{q}|=|q|\\}\subset\mathbb{R}^{6}$

into other points of $TS_{|q|}^{2}$. Moreover, with standard arguments (see
[42]), 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

$\text{Orb}(m)=TS_{|q|}^{2}.$

In particular, when $q\in\mathbb{R}^{3}$ is a unit vector (i.e. $q\in S^{2}$),
$\psi$ allows us to define a transitive Lie group action on
$TS^{2}=TS_{|q|=1}^{2}$ which writes

$\psi:SE(3)\times TS^{2}\rightarrow TS^{2}$
$\psi((A,a),(q,\omega)):=\psi_{(A,a)}(q,\omega)=(Aq,A\omega+\hat{a}Aq)=(\bar{q},\bar{\omega}).$

To conclude the description of the action, we report here its infinitesimal
generator which is fundamental in the Lie group integrators setting

$\left.\psi_{*}((u,v))\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(3)$
and equipped with the direct product structure. More precisely, we consider
the group $G=(SE(3))^{N}$ and by direct product structure we mean that for any
pair of elements

$\delta^{(1)}=(\delta^{(1)}_{1},...,\delta^{(1)}_{N}),\quad\delta^{(2)}=(\delta^{(2)}_{1},...,\delta^{(2)}_{N})\in
G,$

denoted with $*$ the semidirect product of $SE(3)$, we define the product
$\circ$ on $G$ as

$\delta^{(1)}\circ\delta^{(2)}:=(\delta^{(1)}_{1}*\delta^{(2)}_{1},...,\delta^{(1)}_{N}*\delta^{(2)}_{N})\in
G.$

With this group structure defined, we can generalize the action introduced for
$N=1$ to larger $N$s as follows

$\psi:(SE(3))^{N}\times(TS^{2})^{N}\rightarrow(TS^{2})^{N},$
$\begin{split}\psi&((A_{1},a_{1},...,A_{N},a_{n}),(q_{1},\omega_{1},...,q_{N},\omega_{N}))=\\\
&=(A_{1}q_{1},A_{1}\omega_{1}+\hat{a}_{1}A_{1}q_{1},...,A_{N}q_{N},A_{N}\omega_{N}+\hat{a}_{N}A_{N}q_{N}),\end{split}$

whose infinitesimal generator writes

$\psi_{*}(\xi)|_{m}=(\hat{u}_{1}q_{1},\hat{u}_{1}\omega_{1}+\hat{v}_{1}q_{1},...,\hat{u}_{N}q_{N},\hat{u}_{N}\omega_{N}+\hat{v}_{N}q_{N}),$

where $\xi=[u_{1},v_{1},...,u_{N},v_{N}]\in\mathfrak{se}(3)^{N}$ and
$m=(q_{1},\omega_{1},...,q_{N},\omega_{N})\in(TS^{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.

###  5.2 Full chain

We consider the vector field $F\in\mathfrak{X}((TS^{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:(TS^{2})^{N}\rightarrow\mathfrak{se}(3)^{N}$ such that

$\psi_{*}(f(m))|_{m}=F|_{m},\;\;\forall m\in(TS^{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},...,q_{N},\omega_{N})\in(TS^{2})^{N}$, writes

$L(q,\omega)=T(q,\omega)-U(q)=\frac{1}{2}\sum_{i,j=1}^{N}\Big{(}M_{ij}\omega_{i}^{T}\hat{q}_{i}^{T}\hat{q}_{j}\omega_{j}\Big{)}-\sum_{i=1}^{N}\Big{(}\sum_{j=i}^{N}m_{j}\Big{)}gL_{i}e_{3}^{T}q_{i},$

where

$M_{ij}=\Big{(}\sum_{k=\text{max}\\{i,j\\}}^{N}m_{k}\Big{)}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}\Big{(}\omega_{i}^{T}R(q)_{ij}\omega_{j}\Big{)}$

where $R(q)\in\mathbb{R}^{3N\times 3N}$ is a symmetric block matrix defined as

$R(q)_{ii}=\Big{(}\sum_{j=i}^{N}m_{j}\Big{)}L_{i}^{2}I_{3}\in\mathbb{R}^{3\times
3},$
$R(q)_{ij}=\Big{(}\sum_{k=j}^{N}m_{k}\Big{)}L_{i}L_{j}\hat{q}_{i}^{T}\hat{q}_{j}\in\mathbb{R}^{3\times
3}=R(q)_{ji}^{T},\;i<j.$

The vector field on which we need to work defines the following first-order
ODE

$\begin{split}\dot{q}_{i}&=\omega_{i}\times q_{i},\;i=1,...,N,\\\
R(q)\dot{\omega}&=\left[\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{N}M_{ij}|\omega_{j}|^{2}\hat{q}_{i}q_{j}-\Big{(}\sum_{j=i}^{N}m_{j}\Big{)}gL_{i}\hat{q}_{i}e_{3}\right]_{i=1,...,N}\in\mathbb{R}^{3N}\end{split}$

By direct computation it is possible to see that, for any
$q=(q_{1},...,q_{N})\in(S^{2})^{N}$ and $\omega\in
T_{q_{1}}S^{2}\times...\times T_{q_{N}}S^{2}$, we have

$(R(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...\times T_{q_{N}}S^{2}\rightarrow
T_{q_{1}}S^{2}\times...\times T_{q_{N}}S^{2},A_{q}(\omega):=R(q)\omega.$

We can even notice that $R(q)$ defines a positive-definite bilinear form on
this linear space, since

$\omega^{T}R(q)\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}(\hat{q}_{i}\omega_{i})^{T}M_{ij}(\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...\times T_{q_{N}}S^{2}$, see e.g. [16] 666It follows
from the definition of the inertia tensor, i.e.
$0\leq\tilde{T}(q,\dot{q})=\frac{1}{2}\sum_{i=1}^{N}\Big{(}\sum_{j\geq
i}m_{j}\Big{)}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.. 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}\Big{(}[g_{1},...,g_{N}]^{T}\Big{)}=\begin{bmatrix}h_{1}(q,\omega)\\\
...\\\ h_{N}(q,\omega)\end{bmatrix}=\begin{bmatrix}a_{1}(q,\omega)\times
q_{1}\\\ ...\\\ a_{N}(q,\omega)\times q_{N}\end{bmatrix}$

where

$g_{i}=g_{i}(q,\omega)=\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{N}M(q)_{ij}|\omega_{j}|^{2}\hat{q}_{i}q_{j}-\Big{(}\sum_{j=i}^{N}m_{j}\Big{)}gL_{i}\hat{q}_{i}e_{3},\;i=1,...,N$

and $a_{1},...,a_{N}:(TS^{2})^{N}\rightarrow\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 $(TS^{2})^{N}$ to $(\mathfrak{se}(3))^{N}$ such that

$\psi_{*}(f(q,\omega))|_{(q,\omega)}=F|_{(q,\omega)}$

is the following one

$f(q,\omega)=\begin{bmatrix}\omega_{1}\\\ q_{1}\times h_{1}\\\ ...\\\ ...\\\
\omega_{N}\\\ q_{N}\times
h_{N}\end{bmatrix}\in\mathfrak{se}(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}}^{*}S^{2}\times...T_{q_{N}}^{*}S^{2}\ni\pi=(\pi_{1},...,\pi_{N}),\text{
where }\pi_{i}=-\hat{q}_{i}^{2}\frac{\partial L}{\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{split}\dot{q}_{i}&=\frac{\partial H(q,\pi)}{\partial\pi_{i}}\times
q_{i},\\\ \dot{\pi}_{i}&=\frac{\partial H(q,\pi)}{\partial q_{i}}\times
q_{i}+\frac{\partial H(q,\pi)}{\partial\pi_{i}}\times\pi_{i},\end{split}$

which implies that setting

$f(q,\pi)=\begin{bmatrix}\partial_{q_{1}}H(q,\pi),&\partial_{\pi_{1}}H(q,\pi),&\dots,&\partial_{q_{N}}H(q,\pi),&\partial_{\pi_{N}}H(q,\pi)\end{bmatrix}$

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

$\dot{q}_{1}=\hat{\omega}_{1}q_{1},\quad\dot{q}_{2}=\hat{\omega}_{2}q_{2},$
$R(q)\begin{bmatrix}\dot{\omega}_{1}\\\
\dot{\omega}_{2}\end{bmatrix}=\begin{bmatrix}(-m_{2}L_{1}L_{2}|\omega_{2}|^{2}\hat{q}_{2}+(m_{1}+m_{2})gL_{1}\hat{e}_{3})q_{1}\\\
(-m_{2}L_{1}L_{2}|\omega_{1}|^{2}\hat{q}_{1}+m_{2}gL_{2}\hat{e}_{3})q_{2}\end{bmatrix},$
(35)

where

$R(q)=\begin{bmatrix}(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{bmatrix}.$

As presented above, the matrix $R(q)$ 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}\rightarrow
T_{q_{1}}S^{2}\times T_{q_{2}}S^{2},\;[\omega_{1},\omega_{2}]^{T}\rightarrow
R(q)[\omega_{1},\omega_{2}]^{T}.$

We can easily see that it is well defined since

$R(q)\begin{bmatrix}\omega_{1}\\\
\omega_{2}\end{bmatrix}=\begin{bmatrix}(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{bmatrix}\begin{bmatrix}\hat{v}_{1}q_{1}\\\
\hat{v}_{2}q_{2}\end{bmatrix}=\begin{bmatrix}\hat{r}_{1}q_{1}\\\
\hat{r}_{2}q_{2}\end{bmatrix}\in(TS^{2})^{2}$

with

$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}.$

This map guarantees that if we rewrite the pair of equations for the angular
velocities in (35) as

$\begin{split}\dot{\omega}&=R^{-1}(q)\begin{bmatrix}(-m_{2}L_{1}L_{2}|\omega_{2}|^{2}\hat{q}_{2}+(m_{1}+m_{2})gL_{1}\hat{e}_{3})q_{1}\\\
(-m_{2}L_{1}L_{2}|\omega_{1}|^{2}\hat{q}_{1}+m_{2}gL_{2}\hat{e}_{3})q_{2}\end{bmatrix}=R^{-1}(q)b=\\\
&=A_{(q_{1},q_{2})}^{-1}(b)=\begin{bmatrix}h_{1}\\\ h_{2}\end{bmatrix}\in
T_{q_{1}}S^{2}\times T_{q_{2}}S^{2},\end{split}$

then we are assured that there exists a pair of functions
$a_{1},a_{2}:TS^{2}\times TS^{2}\rightarrow\mathbb{R}^{3}$ such that

$\dot{\omega}=\begin{bmatrix}a_{1}(q,\omega)\times q_{1}\\\
a_{2}(q,\omega)\times q_{2}\end{bmatrix}=\begin{bmatrix}h_{1}(q)\\\
h_{2}(q)\end{bmatrix}.$

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

$\begin{bmatrix}\dot{q}_{1}\\\ \dot{\omega}_{1}\\\ \dot{q}_{2}\\\
\dot{\omega}_{2}\end{bmatrix}=\begin{bmatrix}\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{bmatrix}=F|_{(q,\omega)},$

with $h_{i}=h_{i}(q,\omega)$ and

$\begin{bmatrix}h_{1}(q,\omega)\\\
h_{2}(q,\omega)\end{bmatrix}=R^{-1}(q)\begin{bmatrix}(-m_{2}L_{1}L_{2}|\omega_{2}|^{2}\hat{q}_{2}+(m_{1}+m_{2})gL_{1}\hat{e}_{3})q_{1}\\\
(-m_{2}L_{1}L_{2}|\omega_{1}|^{2}\hat{q}_{1}+m_{2}gL_{2}\hat{e}_{3})q_{2}\end{bmatrix}.$

Therefore, we can express the whole vector field in terms of the infinitesimal
generator of the action of $SE(3)\times SE(3)$ as

$\psi_{*}(f(q,\omega))|_{(q,\omega)}=F|_{(q,\omega)}$

through the function

$f:TS^{2}\times
TS^{2}\rightarrow\mathfrak{se}(3)\times\mathfrak{se}(3)\simeq\mathbb{R}^{12},\;\;(q,\omega)\rightarrow(\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 analyzing an
implementation of variable step size. Lie group integrators allow to keep the
evolution of the solution in the correct manifold, which is $TS^{2}\times
TS^{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}(t)^{T}q_{i}(t)=1,\text{ i.e. }q_{i}(t)\in S^{2},\;i=1,2,$

  * •

$q_{i}(t)^{T}\omega_{i}(t)=0,\text{ i.e. }\omega_{i}(t)\in
T_{q_{i}(t)}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 four and one of order three. 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{dexp}_{\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 $TS^{2}\times TS^{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}(t)$ and $q_{2}(t)$, for $0\leq t\leq 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}(t)|=|q_{2}(t)|=1$ $\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.

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

Indeed, in Figure 7 we compare the time evolutions of the inner products
$q_{1}(t)^{T}\omega_{1}(t)$ and $q_{2}(t)^{T}\omega_{2}(t)$ for $t\in[0,5]$,
i.e. we see if these integrators preserve the geometry of the whole phase
space $TS^{2}\times TS^{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 7: Visualization of the inner product $q_{i}(t)^{T}\omega_{i}(t)$,
$i=1,2$, for $t\in[0,5]$. These plots focus on the preservation of the
geometry of $T_{q_{i}(t)}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 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$.

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})=\left(\sqrt{2}/2,0,\sqrt{2}/2,0,1,0\right),\quad\forall
i=1,...,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.

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\rightarrow+\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.

(a) $(q_{1}(t),\omega_{1}(t))$

(b) $(q_{2}(t),\omega_{2}(t))$

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.

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 twelve 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}(t)$,
$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.

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 details this approach, which might be relevant for further
work, however we highlight that here the step adaptivity improves the results
as we expected.

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$.

(a) Step adaptation

(b) Zoom at final times

(c) Values of $\dot{\omega}_{i}^{(2)}(t)$

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 center
we plot the second component of the angular velocities $\omega_{i}^{(2)}$,
$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.

##  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, 50].

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(SO(3))^{2}$ and the directions of the links which connect
the center 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(SO(3))^{2}\times(S^{2})^{2}$.

Figure 12: Two quadrotors connected to the mass point $m_{y}$ via massless
links of lengths $L_{i}$.

In order to present the equations of motion of the system we start by
identifying $TSO(3)\simeq SO(3)\times\mathfrak{so}(3)$ via left-
trivialization. This choice allows us to write the kinematic equations of the
system as

$\dot{R}_{i}=R_{i}\hat{\Omega}_{i},\quad\dot{q}_{i}=\hat{\omega}_{i}q_{i}\quad\quad
i=1,2,$ (36)

where $\Omega_{1},\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},\Omega_{1},R_{2},\Omega_{2},q_{1},\omega_{1},q_{2},\omega_{2}):\mathbb{R}^{6}\times\left(SO(3)\times\mathfrak{so}(3)\right)^{2}\times(TS^{2})^{2}\rightarrow\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}\|\dot{y}\|^{2}+\frac{1}{2}\sum_{i=1}^{2}(m_{i}\|\dot{y}-L_{i}\hat{\omega}_{i}q_{i}\|^{2}+\Omega_{i}^{T}J_{i}\Omega_{i}),$

and

$U=-m_{y}ge_{3}^{T}y-\sum_{i=1}^{2}m_{i}ge_{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 non
conservative 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

$A(z)\dot{z}=h(z)$ (37)

where

$z=[y,v,\Omega_{1},\Omega_{2},\omega_{1},\omega_{2}]^{T}\in\mathbb{R}^{18},$
$A(z)=\begin{bmatrix}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{bmatrix},$
$h(z)=\begin{bmatrix}h_{1}(z)\\\ h_{2}(z)\\\ h_{3}(z)\\\ h_{4}(z)\\\
h_{5}(z)\\\ h_{6}(z)\end{bmatrix}=\begin{bmatrix}v\\\
-\sum_{i=1}^{2}m_{i}L_{i}\|\omega_{i}\|^{2}q_{i}+M_{q}ge_{3}+\sum_{i=1}^{2}u_{i}^{\parallel}\\\
-\Omega_{1}\times J_{1}\Omega_{1}+M_{1}\\\ -\Omega_{2}\times
J_{2}\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{bmatrix},$

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), describe the dynamics of a point

$P=\left[y,\;\;v,\;\;R_{1},\;\;\Omega_{1},\;\;R_{2},\;\;\Omega_{2},\;\;q_{1},\;\;\omega_{1},\;\;q_{2},\;\;\omega_{2}\right]\in
M=TQ.$

Since the matrix $A(z)$ is invertible, we pass to the following set of
equations

$\dot{z}=A^{-1}(z)h(z):=\tilde{h}(z):=\bar{h}(P)=[\bar{h}_{1}(P),...,\bar{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(TSO(3))^{2}\times(TS^{2})^{2}$. The group we consider is

$\bar{G}=\mathbb{R}^{6}\times(TSO(3))^{2}\times(SE(3))^{2},$

where the groups are combined with a direct-product structure and
$\mathbb{R}^{6}$ is the additive group. For a group element

$g=((a_{1},a_{2}),((B_{1},b_{1}),(B_{2},b_{2})),((C_{1},c_{1}),(C_{2},c_{2})))\in\bar{G}$

and a point $P\in M$ in the manifold, we consider the following left action

$\begin{split}\psi_{g}(P)=[y+a_{1},\;\;v+a_{2},\;\;&B_{1}R_{1},\;\;\Omega_{1}+b_{1},\;\;B_{2}R_{2},\;\;\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{split}$

The well-definiteness and transitivity of this action come from standard
arguments, see for example [42]. 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\mathfrak{\bar{g}},$

where $\mathfrak{\bar{g}}=T_{e}\bar{G}$, writes

$\begin{split}\psi_{*}(\xi)|_{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{split}$

We can now focus on the construction of the function
$f:M\rightarrow\bar{\mathfrak{g}}$ such that $\psi_{*}(f(P))|_{P}=F|_{P}$,
where

$\begin{split}F|_{P}=[\bar{h}_{1}(P),\;\;\bar{h}_{2}(P),\;\;R_{1}&\hat{\Omega}_{1},\;\;\bar{h}_{3}(P),\;\;R_{2}\hat{\Omega}_{2},\;\;\\\
&\bar{h}_{4}(P),\;\;\hat{\omega}_{1}q_{1},\;\;\bar{h}_{5}(P),\;\;\hat{\omega}_{2}q_{2},\;\;\bar{h}_{6}(P)]\in
T_{P}M\end{split}$

is the vector field obtained combining the equations (36) and (38). We have

$\begin{split}f(P)=[\bar{h}_{1}(P),\;\;\bar{h}_{2}(P),\;\;R_{1}\Omega_{1},\;\;&\bar{h}_{3}(P),\;\;R_{2}\Omega_{2},\;\;\bar{h}_{4}(P),\\\
\;\;&\omega_{1},\;\;q_{1}\times\bar{h}_{5}(P),\;\;\omega_{2},\;\;q_{2}\times\bar{h}_{6}(P)]\in\bar{\mathfrak{g}}.\end{split}$

We have obtained the local representation of the vector field
$F\in\mathfrak{X}(M)$ in terms of the infinitesimal generator of the
transitive group action $\psi$, hence we can solve for one time step $\Delta
t$ the IVP

$\begin{cases}\dot{\sigma}(t)=\mathrm{dexp}_{\sigma(t)}^{-1}\Big{(}f\big{(}\psi(\exp(\sigma(t)),P(t))\big{)}\Big{)}\\\
\sigma(0)=0\in\bar{\mathfrak{g}}\end{cases}$

and then update the solution $P(t+\Delta t)=\psi(\exp(\sigma(\Delta
t)),P(t))$.

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 [50]. 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 15-18 we show results in the tracking of a parabolic trajectory,
obtained by integrating the system (37) with a RKMK method of order 4.

Figure 14: Snapshots at $0\leq t\leq 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(n),\ n=2,3$, $SE(n),\
n=2,3$ as well as the manifolds $S^{2}$ and $TS^{2}$. The abstract
formulations by e.g. Crouch and Grossman [11], Munthe-Kaas [40] 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 generalisation of the $N$-fold pendulum
discussed here.

## References

  * [1] M. Arnold and O. Brüls.  Convergence of the generalized-$\alpha$ scheme for constrained mechanical systems.  Multibody Syst. Dyn., 18(2):185–202, 2007.
  * [2] M. Arnold, O. Brüls, and A. Cardona.  Error analysis of generalized-$\alpha$ 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 Lie group methods.  Foundations of Computational Mathematics, pages 1–38, 2015.
  * [4] O. Bruls and A. Cardona.  On the Use of Lie Group Time Integrators in Multibody dynamics.  J. Computational Nonlinear Dynamics, 5(3):031002, 2010.
  * [5] F. Casas and B. Owren.  Cost Efficient Lie Group Integrators in the RKMK Class.  BIT Numerical Mathematics, 43(4):723–742, 2003.
  * [6] E. Celledoni, A. Marthinsen, and B. Owren.  Commutator-free Lie group methods.  Future Generation Computer Systems, 19:341–352, 2003.
  * [7] E. Celledoni, H. Marthinsen, and B. Owren.  An introduction to Lie group integrators—basics, new developments and applications.  J. Comput. Phys., 257(part B):1040–1061, 2014.
  * [8] E. Celledoni and B. Owren.  Lie group methods for rigid body dynamics and time integration on manifolds.  Comput. Methods Appl. Mech. Engrg., 192(3-4):421–438, 2003.
  * [9] J. Ćesić, V. Jukov, I. Petrovic, and D. Kulić.  Full body human motion estimation on lie groups using 3D marker position measurements.  In Proceedings of the IEEE-RAS International Conference on HumanoidRobotics, Cancun, Mexico, 15–17 November 2016, 2016.
  * [10] S.H. Christiansen, H.Z. Munthe-Kaas, and B. Owren.  Topics in structure-preserving discretization.  Acta Numerica, 20:1–119, 2011.
  * [11] P. E. Crouch and R. Grossman.  Numerical integration of ordinary differential equations on manifolds.  J. Nonlinear Sci., 3:1–33, 1993.
  * [12] C. Curry and B. Owren.  Variable step size commutator free Lie group integrators.  Numer. Algorithms, 82(4):1359–1376, 2019.
  * [13] F. Diele, L. Lopez, and R. Peluso.  The Cayley transform in the numerical solution of unitary differential systems.  Adv. Comput. Math., 8(4):317–334, 1998.
  * [14] 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\.
  * [15] K. Engø and S. Faltinsen.  Numerical integration of Lie–Poisson systems while preserving coadjoint orbits and energy.  SIAM J. Numer. Anal., 39(1):128–145, 2001.
  * [16] H. Goldstein, C.P. Poole, and J. Safko.  Classical Mechanics.  Pearson, 2013.
  * [17] V. Guillemin and S. Sternberg.  The moment map and collective motion.  Ann. Physics, 127:220–253, 1980.
  * [18] 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 differential equations, Reprint of the second (2006) edition.
  * [19] E. Hairer, S. P. Nørsett, and G. Wanner.  Solving Ordinary Differential Equations I, Nonstiff Problems.  Springer-Verlag, Second revised edition, 1993.
  * [20] J. Hall and M. Leok.  Lie group spectral variational integrators.  Found. Comput. Math., 17(1):199–257, 2017.
  * [21] F. Hausdorff.  Die symbolische Exponentialformel in der Gruppentheorie.  Leipziger Ber., 58:19–48, 1906.
  * [22] H. M. Hilber, T. J. R. Hughes, and R. L. Taylor.  Improved numerical dissipation for time integration algorithms in structural dynamics.  Earthquake Engineering & Structural Dynamics, 5(3):283–292, 1977\.  cited By 1683.
  * [23] D. Holm.  Geometric Mechanics: Part II: Rotating, Translating and Rolling.  World Scientific Publishing Company, 2008.
  * [24] D. Holm, J. Marsden, and T. Ratiu.  The Euler–poincaré equations and semidirect products with applications to continuum theories.  Adv. in Math., 137:1–81, 1998.
  * [25] S. Holzinger and J. Gerstmayr.  Time integration of rigid bodies modelled with three rotation parameters.  Multibody Sys Dyn, 1–34, 2021.
  * [26] A. Iserles, H. Z. Munthe-Kaas, S. P. Nørsett, and A. Zanna.  Lie-group methods.  Acta Numerica, 9:215–365, 2000.
  * [27] T. Lee, M. Leok, and N. H. McClamroch.  Lie group variational integrators for the full body problem.  Comput. Methods Appl. Mech. Engrg., 196(29-30):2907–2924, 2007\.
  * [28] T. Lee, M. Leok, and N. H. McClamroch.  Global formulations of Lagrangian and Hamiltonian dynamics on manifolds.  Interaction of Mechanics and Mathematics. Springer, Cham, 2018.  A geometric approach to modeling and analysis.
  * [29] T. Leitz and S. Leyendecker.  Galerkin Lie-group variational integrators based on unit quaternion interpolation.  Comput. Methods Appl. Mech. Engrg., 338:333–361, 2018.
  * [30] N. E. Leonard and J. E. Marsden.  Stability and drift of underwater vehicle dynamics: Mechanical systems with rigid motion symmetry.  Physica D, 105(1–3):130–162, 1997.
  * [31] D. Lewis and J. C. Simo.  Conserving algorithms for the dynamics of Hamiltonian systems of Lie groups.  J. Nonlinear Sci., 4:253–299, 1994.
  * [32] A. Lundervold and H. Z. Munthe-Kaas.  On algebraic structures of numerical integration on vector spaces and manifolds.  In Faà di Bruno Hopf algebras, Dyson-Schwinger equations, and Lie-Butcher series, volume 21 of IRMA Lect. Math. Theor. Phys., pages 219–263. Eur. Math. Soc., Zürich, 2015.
  * [33] J. E. Marsden, T. Ratiu, and A. Weinstein.  Semi-direct products and reduction in mechanics.  Transactions of the American Mathematical Society, 281:147–77, 1884\.
  * [34] J. E. Marsden and T. S. Ratiu.  Introduction to Mechanics and Symmetry.  Springer-Verlag, 1994.
  * [35] J. E. Marsden and T. S. Ratiu.  Introduction to Mechanics and Symmetry.  Number 17 in Texts in Applied Mathematics. Springer-Verlag, second edition, 1999.
  * [36] R. H. Merson.  An operational method for the study of integration processes.  In Proc. Symp. Data Processing, 1957.
  * [37] J. Moser and A. P. Veselov.  Discrete versions of some classical integrable systems and factorization of matrix polynomials.  Comm. Math. Phys., 139(2):217–243, 1991.
  * [38] A. Müller.  Coordinate mappings for rigid body motions.  ASME Journal of Computational and Nonlinear Dynamics, 12, 2017.
  * [39] H. Munthe-Kaas.  Runge–Kutta methods on Lie groups.  BIT, 38(1):92–111, 1998.
  * [40] H. Munthe-Kaas.  High order Runge–Kutta methods on manifolds.  Appl. Numer. Math., 29:115–127, 1999.
  * [41] H. Munthe-Kaas and B. Owren.  Computations in a free Lie algebra.  Phil. Trans. Royal Soc. A, 357:957–981, 1999.
  * [42] P. J. Olver.  Applications of Lie groups to differential equations, volume 107\.  Springer Science & Business Media, 2000.
  * [43] B. Owren.  Order conditions for commutator-free Lie group methods.  J. Phys. A, 39(19):5585–5599, 2006.
  * [44] B. Owren.  Lie group integrators.  In Discrete mechanics, geometric integration and Lie-Butcher series, volume 267 of Springer Proc. Math. Stat., pages 29–69. Springer, Cham, 2018.
  * [45] B. Owren and A. Marthinsen.  Integration methods based on canonical coordinates of the second kind.  Numer. Math., 87(4):763–790, 2001.
  * [46] J. Park and W. Chung.  Geometric integration on Euclidean group with application to articulated multibody systems.  J. CAM, 21, 2005.
  * [47] T. Ratiu.  Euler-Poisson equations on Lie algebras and the n-dimensional heavy rigid body.  Proc. Nat. Acad. Sci. USA, pages 1327–1328, 1981.
  * [48] A. Saccon.  Midpoint rule for variational integrators on Lie groups.  Internat. J. Numer. Methods Engrg., 78(11):1345–1364, 2009.
  * [49] J. C. Simo and L. Vu-Quoc.  On the dynamics of finite-strain rods undergoing large motions — a geometrically exact approach.  Comput. Methods Appl. Mech. Engrg., 66:125–161, 1988.
  * [50] Lee T., K. Sreenath, and V. Kumar.  Geometric control of cooperating multiple quadrotor uavs with a suspended payload.  52nd IEEE Conference on Decision and Control, pages 5510–5515, 2013\.
  * [51] A. P. Veselov.  Integrable systems with discrete time, and difference operators.  Funktsional. Anal. i Prilozhen., 22(2):1–13, 96, 1988.
  * [52] A. M. Vinogradov and B. Kupershmidt.  The structure of Hamiltonian mechanics.  Funktsional. Anal. i Prilozhen., 32:177–243, 1977.
  * [53] F. W. Warner.  Foundations of Differentiable Manifolds and Lie Groups.  GTM 94. Springer-Verlag, 1983.
  * [54] V. Wieloch and M. Arnold.  BDF integrators for constrained mechanical systems on Lie groups.  J. CAM, 387:112517, 2019, 2019.