Standing on the shoulders of giants – Reboot

Posted on by arkajad

Standing on the shoulders of giants may be good, even very-very good. But at the same time it can be very dangerous. First of all it is a dangerous balancing act. It is easy to fall.

But what if the giant gets a hiccup? Or stumble? You may die. Or someone can die because of aftershocks.

In Standing on the shoulders of giants I quoted two papers that I used in working out my computer simulations of Dzhanibekov’s effect. These were

Ramses van Zon, Jeremy Schofield, “Numerical implementation of the exact dynamics of free rigid bodies“, J. Comput. Phys. 225, 145-164 (2007)
and
Celledoni, Elena; Zanna, Antonella.E Celledoni, F Fassò, N Säfström, A Zanna, The exact computation of the free rigid body motion and its use in splitting methods, SIAM Journal on Scientific Computing 30 (4), 2084-2112 (2008)

The algorithms were working smoothly. But recently I decided to use quaternions instead of rotation matrices. And in recent post Introducing geodesics, I have presented a little animation:

It looks impressive, but if we pay attention to details, we see something strange:

We have a nice curly trajectory, but there are also strange straight line spikes. On animation they are, perhaps, harmless. But if something like this happens with the software controlling the flight of airplanes, space rockets or satellites – people may die. There is a BUG in the algorithm. Yes, I did stay on shoulders of giants, but perhaps they were not the right giants for my aims? And indeed this happens. The method I used is not appropriate for working with quaternions.

That is why I have to do the reboot, and look for another giant. Here is my giant:
Antonella Zanna, professor of Mathematics, University of Bergen, specialty Numerical Analysis/Geometric Integration.

I am going to use the following paper:

Elena Celledoni, Antonella Zanna, et al. “Algorithm 903: FRB–Fortran routines for the exact computation of free rigid body motions.” ACM Transactions on Mathematical Software, Volume 37 Issue 2, April 2010, Article No. 23, doi:10.1145/1731022.1731033.

The paper is only seven years old. And it has all what we need. As in Standing on the shoulders of giants we will split the attitude matrix Q into a product: Q=Q_1Q_0, but this time we will do it somewhat smarter. Celledoni and Zanna denote the angular momentum vector in the body frame using letter \mathbf{m}=(m_1,m_2,m_3), and they assume that \mathbf{m} is normalized: m=\sqrt{\mathbf{m}^2}=1. They assume, as we do, I_1<I_2<I_3, and use T to denote the kinetic energy. Let us first look at the following part of their paper (which is available here):

Though it is not very important, I must say that it is very nice of the two authors that they acknowledge Kosenko. I was trying to study Kosenko’s paper, but I gave up not being able to understand it. Then wen notice that the case considered is that of d>1/I_2 using the notation in my posts. That means “high energy regime”.
The angular momentum \mathbf{m} is constant (fixed) in the laboratory frame owing to the conservation of angular momentum. Suppose we orient the laboratory frame so that its z axis is oriented along the angular momentum vector. Then, in the laboratory frame angular momentum has coordinates (0,0,1). Suppose we construct a matrix Q_0 that transforms \mathbf{m}=(m_1,m_2,m_3) into (0,0,1). The attitude matrix Q transforms vector components in the body frame into their components in the laboratory frame. In particular it transforms \mathbf{m}=(m_1,m_2,m_3) into (0,0,1). Therefore,if we write Q=Q_1Q_0, the matrix Q_0 must transform (0,0,1) into itself. Therefore Q_0 must be a rotation about the third axis, that is it is of the form:

(1)   \begin{equation*}Q_1= \begin{bmatrix}\cos\psi&-\sin\psi&0\\\sin\psi&\cos\psi&0\\0&0&1\end{bmatrix}.\end{equation*}

All of that we have discussed already in Standing on the shoulders of giants. But now we are going to choose Q_0 differently.

So suppose we have vector \mathbf{m} of unit length. What would be the most natural way of constructing a rotation matrix that rotates \mathbf{m} into (0,0,1)? Simple. Project \mathbf{m} vertically onto (x,y) plane. Draw a line \mathbf{n} on (x,y) plane that is perpendicular to the projection. Rotate about this line by the angle \theta between \mathbf{m} and \mathbf{k}=(0,0,1).
Let us now do the calculations. Vector \mathbf{m} has components (m_1,m_2,m_3). Its projection on (x,y) plane has components (m_1,m_2,0). Orthogonal vector in (x,y) plane has components (m_2,-m_1,0) (to check orthogonality calculate scalar product), So the normalized vector \mathbf{n} has components

    \[\mathbf{n}=(\frac{m_2}{m_p}, -\frac{m_1}{m_p},0),\]

where

    \[m_p=\sqrt{m_1^2+m_2^2}.\]

The cosinus of the angle \theta is m_3. Its sinus is \sqrt{1-m_3^2}=m_p.
In Spin – we know that we do not know we have derived a general formula for a rotation about an angle \theta about axis defined by a unit vector \mathbf{n}:

    \[R=I+\sin\theta\, W(\vec{n})+(1-\cos\theta)\,W(\vec{n})^2,\]

where for any vector \mathbf{v}

(2)   \begin{equation*}W=W({\vec{v})=\begin{pmatrix}0& -v_3&v_2\\v_3&0&-v_1\\-v_2&v_1&0.\end{pmatrix}\end{equation*}

Applying the formula to our case, simple algebra leads to:

(3)   \begin{equation*}Q_0=\begin{bmatrix}(m_2^2 + m_1^2 m_3)/m_p^2& (m_1 m_2 (-1 + m_3))/m_p^2& -m_1\\ (m_1 m_2 (-1 + m_3))/m_p^2& ( m_1^2 + m_2^2 m_3)/m_p^2& -m_2\\m_1& m_2& m_3\end{bmatrix}.\end{equation*}

We can easily check that indeed Q_0 is an orthogonal matrix with determinant 1, Q0 acting on (m_1,m_2,m_3) vector gives (0,0,1). One can also check that the vector (-m_2,m_1,0) is invariant under Q_0 – as it should be as it is on the rotation axis. REDUCE program that checks it all is here.

At this point we have to make a break. It is not yet clear what is the relation of my Q_0 to the paper quoted, what is the difference between this and the old version, and why is this version better?

We will continue climbing on the shoulders of giants in the following notes. But, please, remember, Standing on shoulders of giants (like “Nobel Prize Winner” Obama), is risky:

Photo of Kellyanne Conway kneeling on Oval Office couch sparks…

7 thoughts on “Standing on the shoulders of giants – Reboot”

  1. Ark,
    today the World delivered me the thought:
    “Did you finish with your sequence?”

    Under your post Never twice the same
    I asked:
    “\pi/4,\;\pi/3,\;\pi/6
    Where is the limit of this sequence?”

    Now the sequence looks like:
    \pi/4,\;\pi/3,\;\pi/6;\pi/3

    But is it finished?
    From symmetry there should be fifth term:
    \pi/3 or \frac{4}{3}\pi ?

    Do you agree?

    1. Have to think about it. Today I have first to understand what the giants say. Could not sleep.

      Why don’t YOU finish the sequence while I am busy with the giants?

      1. Take a break from the giants.

        I gave the proposition of the fifth term.

        A bug in your post “Never twice the same” is deeply hidden and I couldn’t find it so far.

  2. “Applying the formula to our case, simple algebra leads to:”

    I applied formulas to our case and I got something wich is not as your Q_0 but Q_0^T

    So maybe we should take not Rodrigues formula for R but for R^T

    1. So, we have two different results. And only one of them is correct. I can’t find mistake in my calculations. But will continue searching until the issue get resolved.

Leave a Reply