Canonical rotations that do not produce computational singularities

Intro

On the topic of dynamical systems associated with 3-dimensional rotation of rigid bodies, you will always encounter singularities in the equations of motion that will produce computational errors at certain rotations. At least, these singularities will always occur if the following approach is attempted (see next section).

Problem

Suppose you want to describe the coordinates $(x,y,z)$ of a rigid body by spherical parameterization using two canonical coordinates $q_1$ and $q_2$ as spherical coordinates

$$vec r = left( {begin{array}{*{20}{c}}
{rsin q_1 cos q_2 }
{rsin q_1 sin q_2 }
{rcos q_1 }
end{array}} right).$$

You can easily obtain the velocity vector $vec {dot r}$ and then associate it with the kinetic energy
$$T = frac{1}{2}m{left| {vec {dot r}} right|^2}.$$
To further simplify the problem, assume the potential energy to be 0. The Lagrangian $L = T – V$ and Hamiltonian $H = T + V$ can then be obtained.

From the Lagrangian we can obtain the equations of motion
$$begin{array}{l}
{{ddot q}_1} = frac{1}{2}dot q_2^2sin left( {2{q_1}} right)
{{ddot q}_2} = 2{{dot q}_1}{{dot q}_2}sin left( {2{q_1}} right)frac{1}{{cos left( {2{q_1}} right) – 1}}.
end{array}$$
It is clear that the equations of motion are ill-defined (in a computational sense) for $q_1 = npi$, which is at the poles of the sphere.

Let us examine the Hamiltonian equations to verify that the issue persists. Using canonical momenta $p_1$ and $p_2$, we find:

$$begin{array}{l}
{p_1} = {{dot q}_1}m{r^2}
{p_2} = – frac{1}{2}{{dot q}_2}m{r^2}left( {cos left( {2{q_1}} right) – 1} right)
H = – left( {p_1^2 – p_1^2cos left( {2{q_1}} right) + 2p_2^2} right)frac{1}{{2m{r^2}( cos left( {2{q_1} } right) – 1)}}
frac{{partial H}}{{d{q_1}}} = – 2p_2^2sin left( {2{q_1}} right)frac{1}{{m{r^2}{{left( {cos left( {2{q_1}} right)} -1 right)}^2}}}
frac{{partial H}}{{d{q_2}}} = 0
frac{{partial H}}{{d{p_1}}} = frac{{{p_1}}}{{m{r^2}}}
frac{{partial H}}{{d{p_2}}} = – frac{{2{p_2}}}{{m{r^2}left( {cos left( {2{q_1}} right) – 1} right)}}.
end{array}$$

The problem persists: The equations of motions are always ill-defined (from a computational point of view) at the poles. Note that the limit is never diverging and the solution makes perfect sense in a physical point of view. From a computational point of view, any numerical simulation will produce large numerical errors if the rotation $q_1$ approaches $npi$. This would be a problem if the system was a spherical pendulum for instance.

The same issue occurs using cylindrical coordinates. Using other alternative representations of $vec r$ does not appear to solve the issue either.

Question

Is it possible to describe the equations of motion for rigid body rotation without having to deal with the threat of singularities?

I am aware you can use conditional coordinate transformations during simulation to address the issue. This poses some other issues for multistep integrators, and I therefore do not desire this option.

How would you approach the problem of describing the equations for rotating rigid bodies without computational singularities, specifically preventing that the Hamiltonian gradient and the canonical coordinates approaches infinity at any point?

• Perhaps using a different parameterization of the canonical coordinates?
• Perhaps using more than two canonical coordinates to describe the system?
• Perhaps making the system time-dependent?
• Perhaps taking advantage of the cyclic variable – in this case $q_2$ with $frac{{partial H}}{{d{q_2}}} = 0$ to somehow reduce the system
• Perhaps some transformations of coordinates
• Perhaps using quaternions (suggested by G. Smith, JEB and suggested transformation by Eli). How would you describe the equations of motions with quaternions without introducing singularities?

Open for any ideas

To further specify what I seek: A set of generalized coordinates from which the equations of motion associated with the rotation can be expressed without intrusive divergence or singularities.

I pursue suggestions that can be utilized also for more complicated rigid body systems with non-zero potential energy. For instance the double spherical pendulum and other systems with cyclic properties.

Progress

1) Quaternions (proposed in comments)

Quaternions coordinates were tested. The Hamiltonian and Lagrangian equations of motion can be expressed in terms of three canonical/generalized quaternion coordinates. Unfortunately divergence occurs $a^2 + b^2 + c^2 = 1$ due to division of the fourth coordinate $d = 0$. Divergence occurs more often than systems of two canonical coordinates, and this approach is therefore not pursued any further.

2) epsilon to remove the singularity (proposed by Eli)

A small value $epsilon$ can be added to the divisor in either the Hamiltonian or Lagrangian equations of motion, preventing division by 0. This can be done using both spherical and cylindrical coordinates.

Unfortunately this method changes the system noticeable when $epsilon$ is large. If $epsilon$ is small, an error similar to the effect of the singularity will be produced. This solution is very practical, but not ideal for accurate depictions of the system.

3) Variable axis of singularity

We can construct the coordinate system of the spherical pendulum (or other cyclic systems) using two pairs of coordinate $vec r_x$ and $vec r_y$. When the x-coordinate approaches its maximum and minimum value, a singularity is appraoched for coordinate system $vec r_x$. Likewise for the y-coordinate. The idea is then to transform the coordinates from the $vec r_x$ to $vec r_y$ whenever the x-singularity is approached, and transform form $vec r_y$ to $vec r_x$ whenever the y-singularity is approached.

$${vec r_x} = left( {begin{array}{*{20}{c}}
{rcos left( {{q_1}} right)}
{rsin left( {{q_1}} right)sin left( {{q_2}} right)}
{rsin left( {{q_1}} right)cos left( {{q_2}} right)}
end{array}} right),{vec r_y} = left( {begin{array}{*{20}{c}}
{rsin left( {{q_1}} right)sin left( {{q_2}} right)}
{rcos left( {{q_1}} right)}
{rsin left( {{q_1}} right)cos left( {{q_2}} right)}
end{array}} right)$$

You can transform coordinates from one system to the other using
$$left( {begin{array}{*{20}{c}}
{{q_1}}
{{q_2}}
end{array}} right) to left( {begin{array}{*{20}{c}}
{{{cos }^{ – 1}}left( {sin left( {{q_1}} right)sin left( {{q_2}} right)} right)}
{{rm{atan2}}left( {cosleft( {{q_1}} right),sin left( {{q_1}} right)cos left( {{q_2}} right)} right)}
end{array}} right)$$

https://youtu.be/5bh_dMn-Plc

Above video above illustrates this method, where the simulation switches between $vec r_x$ (red sphere) and $vec r_y$ (blue sphere).
The simulation never reaches singularities, but the transformation between the two systems is not symplectic and produces a linear energy dissipation (energy error depicted top-right in the video).
The transformation needs to be done in a different computational way. Investigating this approach further in hope of a proper solution.

4) Other ideas?

If you have another idea – either a new approach or one based on the above, please share it. This post will be updated when further progress is made.