The double spherical pendulum
In the paper (Celledoni, Çokaj, Leone, Murari and Owren, (2021) International Journal of Computer Mathematics) we discuss in detail the N-fold spherical pendulum, which is an example of the initial value problem discussed in the first section.
We consider here the special case of the double spherical pendulum, \(N = 2\).
We then test our Lie group integrators on this particular case, showing the rate of convergence, the preservation of the geometry of the manifold \(S^2\) and the phase space \(T_{q_{i}(t)}S^2\).
The Lagrangian we consider is a function from \((TS^2)^2\) to \(\mathbb{R}\). Instead of the coordinates \((q_1, q_2,\dot{q}_1, \dot{q}_2)\), where \(\dot{q}_i\in T_{q_i}S^2\), we choose to work with the angular velocities.
Precisely,
(1)\[\begin{align}
T_{q_i}S^2 = \{v\in\mathbb{R}^3:\;v^Tq_i=0\} = \langle q_i\rangle ^{\perp} \subset \mathbb{R}^3, \quad i = 1,2
\end{align}\]
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^Tq_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)\in (TS^2)^2\) to describe the dynamics.
We denote with \(m_1, m_2\) the masses of the pendulums and with \(L_1, L_2\) their lengths.
Transitive group action on \((TS^2)^2\)
In (Celledoni, Çokaj, Leone, Murari and Owren, (2021) International Journal of Computer Mathematics) 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)\), which is
(2)\[\begin{align}
\textrm{Ad} : SE(3)\times \mathfrak{se}(3) \rightarrow \mathfrak{se}(3),
\end{align}\]
(3)\[\begin{align}
\textrm{Ad}((R,r),(u,v)) = (Ru,Rv+\hat{r}Ru).
\end{align}\]
Since \(\mathfrak{se}(3)\simeq \mathbb{R}^6\), the Ad-action allows us to define the following Lie group action on \(\mathbb{R}^6\)
(4)\[\begin{align}
\Psi: SE(3)\times\mathbb{R}^6\rightarrow \mathbb{R}^6,\;\;\Psi((R,r),(u,v)) = (Ru,Rv+\hat{r}Ru).
\end{align}\]
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
(5)\[\begin{align}
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
\end{align}\]
into other points of \(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 is
(6)\[\begin{align}
\Psi : SE(3)\times TS^2 \rightarrow TS^2
\end{align}\]
(7)\[\begin{align}
\Psi((A,a),(q,\omega)) := \Psi_{(A,a)}(q,\omega) = (Aq,A\omega + \hat{a}Aq)=(\bar{q},\bar{\omega}).
\end{align}\]
To conclude the description of the action, we report here its infinitesimal generator which is fundamental in the Lie group integrators setting
(8)\[\begin{align}
\Psi_*((u,v))\|_{(q,\omega)} =(\hat{u}q,\hat{u}\omega + \hat{v}q).
\end{align}\]
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.
Here we limit ourselves to the case \(N=2\) for which we also show numerical experiments.
More precisely, we consider the group \(G=(SE(3))^2\) and by direct product structure we mean that for any pair of elements
(9)\[\begin{align}
\delta^{(1)}=(\delta^{(1)}_1, \delta^{(1)}_2), \delta^{(2)}=(\delta^{(2)}_1, \delta^{(2)}_2)\in G,
\end{align}\]
denoted with \(*\) the semidirect product of \(SE(3)\), we define the product \(\circ\) on \(G\) as
(10)\[\begin{align}
\delta^{(1)}\circ \delta^{(2)} := (\delta^{(1)}_1 * \delta^{(2)}_1, \delta^{(1)}_2 * \delta^{(2)}_2)\in G.
\end{align}\]
With this group structure defined, we can write the action follows
(11)\[\begin{align}
\Psi : (SE(3))^2\times (TS^2)^2 \rightarrow (TS^2)^2,
\end{align}\]
(12)\[\begin{split}\begin{align}
\begin{split}
\Psi&((A_1,a_1, A_2,a_2),(q_1,\omega_1, q_2,\omega_2)) =\\ &=(A_1q_1,A_1\omega_1+\hat{a}_1A_1q_1, A_2q_2,A_2\omega_2+\hat{a}_2A_2q_2),
\end{split}
\end{align}\end{split}\]
whose infinitesimal generator is
(13)\[\begin{align}
\Psi_*(\xi)\vert_m =(\hat{u}_1q_1,\hat{u}_1\omega_1+\hat{v}_1q_1, \hat{u}_2q_2,\hat{u}_2\omega_2+\hat{v}_2q_2),
\end{align}\]
where \(\xi=[u_1,v_1, u_2,v_2]\in\mathfrak{se}(3)^2\) and \(m=(q_1,\omega_1, q_2,\omega_2)\in (TS^2)^2\).
We have now the only group action we need to deal with the double 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.
The equations of motion and the vector field
We consider the vector field \(F\in\mathfrak{X}((TS^2)^2)\), describing the dynamics of the double spherical 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)^2\rightarrow \mathfrak{se}(3)^2\) such that
(14)\[\begin{align}
\Psi_*(f(m))\vert_m = F\vert_m,\;\;\forall m\in (TS^2)^2.
\end{align}\]
The derivation of \(F\) starting from the Lagrangian of the system can be found in the section devoted to mechanical systems on \((S^2)^2\) of (Lee, Leok and McClamroch, (2018)).
The configuration manifold of the system is \((S^2)^2\), while the Lagrangian, expressed in terms of the variables \((q_1,\omega_1, q_2,\omega_2)\in (TS^2)^2\), is
(15)\[\begin{align}
L(q,\omega) = T(q,\omega)-U(q) =\frac{1}{2}\sum_{i,j=1}^2\Big(M_{ij}\omega_i^T\hat{q}_i^T\hat{q}_j\omega_j\Big) - \sum_{i=1}^2\Big(\sum_{j=i}^2 m_j\Big)gL_ie_3^Tq_i,
\end{align}\]
where
(16)\[\begin{align}
M_{ij} =\Big(\sum_{k=\textrm{max}\{i,j\}}^2 m_k\Big)L_iL_j I_3\in\mathbb{R}^{3\times 3}
\end{align}\]
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
(17)\[\begin{align}
\omega_i^T\hat{q}_i^T\hat{q}_i\omega_i = \omega_i^T(I_3-q_iq_i^T)\omega_i = \omega_i^T\omega_i,
\end{align}\]
we simplify the notation writing
(18)\[\begin{align}
T(q,\omega) = \frac{1}{2}\sum_{i,j=1}^2\Big(\omega_i^TR(q)_{ij}\omega_j\Big)
\end{align}\]
where \(R(q)\in\mathbb{R}^{6\times 6}\) is a symmetric block matrix defined as
(19)\[\begin{align}
R(q)_{ii} = \Big(\sum_{j=i}^2m_j\Big)L_i^2I_3\in\mathbb{R}^{3\times 3},
\end{align}\]
(20)\[\begin{align}
R(q)_{ij} = \Big(\sum_{k=j}^2 m_k\Big)L_iL_j\hat{q}_i^T\hat{q}_j\in\mathbb{R}^{3\times 3} = R(q)_{ji}^T,\; i<j.
\end{align}\]
Precisely, the equations of motion write:
(21)\[\begin{align}
\dot{q}_1 = \hat{\omega}_1q_1,\quad \dot{q}_2 = \hat{\omega}_2q_2,
\end{align}\]
(22)\[\begin{split}\begin{align}
R(q)\begin{bmatrix}
\dot{\omega}_1 \\ \dot{\omega}_2
\end{bmatrix}=
\begin{bmatrix}
(-m_2L_1L_2|\omega_2|^2\hat{q}_2 + (m_1+m_2)gL_1\hat{e}_3)q_1 \\
(-m_2L_1L_2|\omega_1|^2\hat{q}_1 + m_2gL_2\hat{e}_3)q_2
\end{bmatrix},
\end{align}\end{split}\]
where
(23)\[\begin{split}\begin{align}
R(q) = \begin{bmatrix}
(m_1+m_2)L_1^2I_3 & m_2L_1L_2\hat{q}_1^T\hat{q}_2 \\
m_2L_1L_2\hat{q}_2^T\hat{q}_1 & m_2L_2^2I_3
\end{bmatrix}.
\end{align}\end{split}\]
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:
(24)\[\begin{align}
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.
\end{align}\]
We can easily see that it is well defined since
(25)\[\begin{split}\begin{align}
R(q)\begin{bmatrix}
\omega_1 \\ \omega_2
\end{bmatrix} = \begin{bmatrix}
(m_1+m_2)L_1^2I_3 & m_2L_1L_2\hat{q}_1^T\hat{q}_2 \\
m_2L_1L_2\hat{q}_2^T\hat{q}_1 & m_2L_2^2I_3
\end{bmatrix}\begin{bmatrix}
\hat{v}_1q_1 \\ \hat{v}_2q_2
\end{bmatrix} = \begin{bmatrix}
\hat{r}_1q_1\\ \hat{r}_2q_2
\end{bmatrix}\in (TS^2)^2
\end{align}\end{split}\]
with
(26)\[\begin{align}
r_1(q,\omega):=(m_1+m_2)L_1^2v_1+m_2L_1L_2\hat{q}_2\hat{v}_2q_2,
\end{align}\]
(27)\[\begin{align}
r_2(q,\omega):=m_2L_1L_2\hat{q}_1\hat{v}_1q_1+m_2L_2^2v_2.
\end{align}\]
This map guarantees that if we rewrite the pair of equations for the angular velocities in (22) as
(28)\[\begin{split}\begin{align}
\begin{split}
\dot{\omega}&= R^{-1}(q)\begin{bmatrix}
(-m_2L_1L_2|\omega_2|^2\hat{q}_2 + (m_1+m_2)gL_1\hat{e}_3)q_1 \\
(-m_2L_1L_2|\omega_1|^2\hat{q}_1 + m_2gL_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}
\end{align}\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
(29)\[\begin{split}\begin{align}
\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}.
\end{align}\end{split}\]
Since we want \(a_i\times q_i = h_i\), \(i = 1,2\), we just impose \(a_i=q_i\times h_i\), \(i = 1,2\), and hence the whole vector field can be rewritten as
(30)\[\begin{split}\begin{align}
\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\vert_{(q,\omega)},
\end{align}\end{split}\]
with \(h_i=h_i(q,\omega)\), \(i = 1,2\), and
(31)\[\begin{split}\begin{align}
\begin{bmatrix}
h_1(q,\omega) \\ h_2(q,\omega)
\end{bmatrix} = R^{-1}(q)\begin{bmatrix}
(-m_2L_1L_2|\omega_2|^2\hat{q}_2 + (m_1+m_2)gL_1\hat{e}_3)q_1 \\
(-m_2L_1L_2|\omega_1|^2\hat{q}_1 + m_2gL_2\hat{e}_3)q_2
\end{bmatrix}.
\end{align}\end{split}\]
Therefore, we can express the whole vector field in terms of the infinitesimal generator of the action of \(SE(3)\times SE(3)\) as
(32)\[\begin{align}
\Psi_*(f(q,\omega))\vert_{(q,\omega)}=F\vert_{(q,\omega)}
\end{align}\]
through the function
(33)\[\begin{align}
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).
\end{align}\]