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


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}} {r\sin q_1 \cos q_2 }\\ {r\sin q_1 \sin q_2 }\\ {r\cos 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^2\sin \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 = n\pi$, 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^2\cos \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^2\sin \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 $n\pi$. 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.


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.


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}} {r\cos \left( {{q_1}} \right)}\\ {r\sin \left( {{q_1}} \right)\sin \left( {{q_2}} \right)}\\ {r\sin \left( {{q_1}} \right)\cos \left( {{q_2}} \right)} \end{array}} \right),{\vec r_y} = \left( {\begin{array}{*{20}{c}} {r\sin \left( {{q_1}} \right)\sin \left( {{q_2}} \right)}\\ {r\cos \left( {{q_1}} \right)}\\ {r\sin \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( {cos\left( {{q_1}} \right),\sin \left( {{q_1}} \right)\cos \left( {{q_2}} \right)} \right)} \end{array}} \right)$$


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.

  • 1
    $\begingroup$ Aren’t quaternions good for avoiding such problems? $\endgroup$
    – G. Smith
    Jul 27, 2020 at 5:01
  • 2
    $\begingroup$ Would Computational Science be a better home for this question? $\endgroup$ Jul 27, 2020 at 7:07
  • $\begingroup$ If I am to follow G. Smith's and JEBs suggestion, the challenge is to express the Hamiltonian or Lagrangian with quaternion canonical coordinates. So it feels more like a physics problem to me. But this problem apparently evaded all physics textbooks and internet resources on rigid body systems... $\endgroup$
    – user270876
    Jul 27, 2020 at 7:17
  • 1
    $\begingroup$ Read page 7 (above) for references to a discussion on rotational representations with and without singularities. $\endgroup$ Aug 1, 2020 at 21:00
  • 1
    $\begingroup$ malcolmdshuster.com/Pub_1993c_J_RotVec_IEEE.pdf is a paper that I think solves this problem, or at least provides useful insight into past efforts. $\endgroup$ Aug 1, 2020 at 21:25

2 Answers 2


It's not just pendulums. It's a major issue for spacecraft navigation, especially if, for instance, your trajectory involves the deployment of a supersonic parachute on Mars.

Hence: quaternions are the standard.

They are also computationally faster, which can be a factor with space-qualified CPUs, which are not fast.

  • 1
    $\begingroup$ Answer to both JEB and G. Smith. You are correct that quaternions are used to describe the rotations without nasty singularities and gimbal lock. However, I am yet to see Hamiltonian systems with non-zero potential being expressed in terms of quaternion coordinates. Do you have an idea how to describe the kinetic energy and potential energy (such as the spherical pendulum with $V = mgz$) in terms of quaternion canonical coordinates? $\endgroup$
    – user270876
    Jul 27, 2020 at 5:25

I don't think you can avoid this singularity by using Quaternions because they are good for rotation .

but for numerical simulation you can use this concept:

your geodetic equations are:

$$\ddot{\theta}=-2\,{\frac {\dot{\theta}\dot{\phi}\cos \left( \phi \right) }{\sin \left( \phi \right) }} \tag 1$$


$$\ddot{\phi}={\frac {\sin \left( \phi \right) \left( r{\dot \theta }^{2}\cos \left( \phi \right) +g \right) }{r}}\tag 2$$

you have singularity if you want to simulate equation (1) ans (2) with initial condition $\phi(0)=0$ , to avoid this singularity substitute in equation (1) $\sin(\phi)\mapsto \sin(\phi+\epsilon)$ where $\epsilon$ is a small number

$$\ddot{\theta}=-2\,{\frac {\dot{\theta}\dot{\phi}\cos \left( \phi \right) }{\sin \left( \phi+\epsilon \right) }} \tag 3$$

if the initial condition $\dot{\theta}(0)$ equal zero or $\dot{\phi}(0)$ equal zero you don't get problem with the initial condition $\phi(0)=0$

This is the simulation result ,geodetic line ,with initial conditions all zero and gravitation g, I choose $\epsilon=0.001$

enter image description here


Sphere Position Vector:

$$\vec{R}=\left[ \begin {array}{c} r\cos \left( \theta \right) \sin \left( \phi \right) \\ r\sin \left( \theta \right) \sin \left( \phi \right) \\ r\cos \left( \phi \right) \end {array} \right] \tag 1$$

where $\theta$ the azimuth coordinate $0\le\theta\le 2\pi$ , $\phi$ the polar coordinate $0\le \phi\le \pi$ and r is the radius of the sphere .

How to avoid the singularity at $\phi=0$ and $\phi=\pi$

substitute in equation (1) $\sin(\phi)\mapsto \sin(\phi+\epsilon)~,$ ( with r=1 )

the kinetic energy is now

$$T=1/2\, \left( \left( \cos \left( \phi+\epsilon \right) \right) ^{2}+1 - \left( \cos \left( \phi \right) \right) ^{2} \right) m{\dot\phi }^{2}+ 1/2\,m \left( 1- \left( \cos \left( \phi+\epsilon \right) \right) ^{2 } \right) {\dot\theta }^{2} $$

to inspect the singularity you obtain the Mass Matrix from the kinetic energy $$M(\theta \,,\phi )=\left[ \begin {array}{cc} {\frac {\partial ^{2}}{\partial {\theta p}^ {2}}}T \left( \theta p,\phi p \right) &{\frac {\partial ^{2}}{ \partial \theta p\partial \phi p}}T \left( \theta p,\phi p \right) \\ {\frac {\partial ^{2}}{\partial \theta p\partial \phi p}}T \left( \theta p,\phi p \right) &{\frac {\partial ^{2}}{ \partial {\phi p}^{2}}}T \left( \theta p,\phi p \right) \end {array} \right] $$

where $\theta p=\dot{\theta}$ and $\phi p=\dot{\phi}$

$$M=\left[ \begin {array}{cc} -m \left( -1+ \left( \cos \left( \phi+ \epsilon \right) \right) ^{2} \right) &0\\ 0& \left( \left( \cos \left( \phi+\epsilon \right) \right) ^{2}+1- \left( \cos \left( \phi \right) \right) ^{2} \right) m\end {array} \right] $$ thus the determinate is: $$\det(M)=-{m}^{2} \left( -1+ \left( \cos \left( \phi+\epsilon \right) \right) ^{2} \right) \left( \left( \cos \left( \phi+\epsilon \right) \right) ^{2}+1- \left( \cos \left( \phi \right) \right) ^{2} \right) $$ take the determinate equal zero and solve for $\phi$ you obtain two real solutions

$$\phi_0=-\epsilon ~,\phi_0=-\epsilon+\pi$$

thus for $\epsilon\ll$ the simulation will work perfect

  • $\begingroup$ Your first suggestion on modifying the equations of motion with eps will not do the trick, as my concern is not the invalid initial condition, but rather the spurious Hamiltonian dissipation that arise when the singularity is approached. $\endgroup$
    – user270876
    Jul 27, 2020 at 19:09
  • 1
    $\begingroup$ @Egeris The Transformation Matrix is correct but the equation of motion is not, because d is also a function of the time, so it is a little bit complicate, I can give you the answer for that but I need your quotation "How to implement the Quanternios " $\endgroup$
    – Eli
    Jul 28, 2020 at 20:06
  • 1
    $\begingroup$ Typo: It is Quaternions not Quaternios. $\endgroup$ Jul 29, 2020 at 12:03
  • 1
    $\begingroup$ I obtained the correct equations from the quaternion rotation matrix though. But it isn't useful unless I can make it work with transformations during the simulation to skip the singularities. Currently having some luck doing simulations transforming between two coordinate pairs of different spherical coordinates. $\endgroup$
    – user270876
    Jul 30, 2020 at 8:32
  • 1
    $\begingroup$ @Eli It seems like a practical approach for approximating the time evolution for simple spherical systems. The challenge for me, is the to obtain $exactly$ the phase space from the time evolution, but I am yet to find a method that achieves this. I'll do a comparative analysis of this approach against my other techniques involving moving axes, once I manage to achieve symplecticity. $\endgroup$
    – user270876
    Aug 3, 2020 at 13:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.