Giants that stumble

Posted on by arkajad

User friendly means “easy to use”. We can all agree on that, I think. And it’s what we mean when we say it. So far, so good.

Except, not. When people say “easy to use”, what they REALLY mean is “easy to learn”, or “easy to use without learning”.
From What User friendly means


I am trying to make my blog user friendly. Sometimes I think I succeed. But not always. We are following giants. It is not easy. Sometimes we need to overcome our fears.

I am studying the paper and the code:

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

I am trying to understand it and to reproduce the results. If I am able to reproduce the results, then there is a good chance that I really understand what I am studying, and then I can share my understanding here, on the blog. Not all I am able to understand, but what to do? Sometimes I am completely lost.

Here is one such case. I am quoting one paragraph from Standing on the shoulders of giants – Reboot

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

In the same post we introduced the matrix Q_0 that transforms vector with components (m_1,m_2,m_3) into (0,0,1). In Attitudes and Behavior we wrote it as (and I am replacing L_1(t) with m_1 etc.):

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

We have the formula unit mapping quaternions into rotation matrices (see Chiromancy in the rotation group).
If

(2)   \begin{equation*}q=W+\mathbf{i}X+\mathbf{j}Y+\mathbf{k}Z\end{equation*}

Then

(3)   \begin{equation*} R(q)=\begin{bmatrix} W^2+X^2-Y^2-Z^2&2(XY-WZ)&2(WY+XZ)\\ 2(WZ+XY)&W^2-X^2+Y^2-Z^2&2(YZ-WX)\\ 2(XZ-WY)&2(WX+YZ)&W^2-X^2-Y^2+Z^2 \end{bmatrix}, \end{equation*}

If we use the formula (21) from the paper (see image above), then

(4)   \begin{eqnarray*} W&=&\sqrt{1+m_3}/\sqrt{2},\\ X&=&\frac{m_2}{\sqrt{2(1+m_3)}},\\ Y&=&-\frac{m_1}{\sqrt{2(1+m_3)}},\\ Z&=&0. \end{eqnarray*}

Now I am using Mathematica to calculate R(q)

We obtain Q_0. Two details will catch the attention of the Reader:
In the Mathematica code I use, for instance, 1-2(Y^2+Z^2) while in the formula for R(q) I have W^2+X^2-Y^2-Z^2. But for the unit quaternion, with W^2+X^2+Y^2+Z^2, these are the same! Similarly the Mathematica code gives as the component (3,3) of the matrix

    \[\frac{1-m1^2-m2^2+m3}{1+m_3},\]

while in Q_0 we have simply m_3. But again, with m_1^2+m_2^2+m_3^2=1, the Reader will easily check that these are the same!
So it seems that we are progressing with our understanding of the language and the science of giants! Success. Therefore we must not stop! We must continue!
After the case m>1 Celledoni and Zanna move to the case of m<1:

So, if we understand it correctly, for m<1 we should replace Q_0 matrix with the new one, constructed the same way as before, but with the third axis being replaced by the first. The quaternion that implements such a rotation, according to the formula (23) from the paper, should be given by

(5)   \begin{eqnarray*} W&=&\sqrt{1+m_1}/\sqrt{2},\\ X&=&0,\\ Y&=&-\frac{m_3}{\sqrt{2(1+m_1)}},\\ Z&=&\frac{m_2}{\sqrt{2(1+m_1)}}.\\ \end{eqnarray*}

I calculate R(q) for this quaternion, and I get …. mess. Certainly not the matrix that acting on (m_1,m_2,m_3) would give (1,0,0). So, something is wrong! I am trying to compare the formulas in the paper with those in the Fortran code. But then I am getting confused even more. For the case m>1 the code has these lines:

APPP = 1.0D0 + EM30
APP = SQRT(APPP)
ADUE=SQRT(2.0D0)
Q2(1) = APP/ADUE
Q2(4) = -0.0D0
Q2(2) = -Q2(1)*EM20/APPP
Q2(3) = Q2(1)*EM10/APPP

We recognize that Q2(1) is our W=\sqrt{1+m_3}/\sqrt{2} in Eqs. (\label{eq{wxyz3}), Q2(4) is our Z. EM10,EM20,EM30 is probably our m_1,m_2,m_3. If so, then Q2(2) is minus our X and Q2(3) is minus our Y. I have no idea why there are minuses, but there are there!

What about the case m<1? The code has these lines:

APPP = 1.0D0 + EM10
APP = SQRT(APPP)
ADUE=SQRT(2.0D0)
Q2(1) = APP/ADUE
Q2(2) = 0.0D0
Q2(3) = -Q2(1)*EM30/APPP
Q2(4) = Q2(1)*EM20/APPP

Again we recognize Q2(1) as our W=\sqrt{1+m_1}/\sqrt{2} we recognize Q2(2) as X=0. But now the signs in Q2(2) and Q2(3) agree with those in the paper. Something weird is going on. So, as an experiment I change the signs in formula (23):

(6)   \begin{eqnarray*} W&=&\sqrt{1+m_1}/\sqrt{2},\\ X&=&0,\\ Y&=&\frac{m_3}{\sqrt{2(1+m_1)}},\\ Z&=&-\frac{m_2}{\sqrt{2(1+m_1)}}.\\ \end{eqnarray*}

So that, consistently, the two signs in the paper are opposite to those in the code. Now R(q) gives what it is supposed to do: rotation matrix that transforms (m_1,m_2,m_3) into (1,0,0).

So, my guess is that the formula (23) in the paper has wrong signs. What to do?
Should I write to the authors? My thinking goes into the direction: it would be of no use.
The paper is 7 years old. In fact the preprint tells me that the work has been essentially finished in 2008. Now the authors are working on a different project, they have other things to do than worrying about signs in a formula. Perhaps the formula is wrong, so what? Who cares? No one is reading paper anyway. Scientists write papers in order to document that they were busy with something, to show that they are experts, to get new grants, to be quoted. Anyway intelligent reader should be able to fix typos. It should even be advised to make some errors. Fixing someone’s errors takes our time, true. But, on the other hand, what a wonderful opportunity to think and to verify that we really understand the subject! Finding errors in pear reviewed papers gives us reassurance that referees in these papers, as a rule, do not pay attention to the question as to whether a given paper is correct or incorrect. It should simply give a strong impression that it is correct.. That is all.
True, some authors take mistakes in their published works seriously. I like to quote this particular very good example: ERRATA TO ANALYSIS, MANIFOLDS AND PHYSICS. PART I, the monograph by Yvonne Choquet-Bruhat and Cecile DeWitt-Morette, is 11 pages long. With all kinds of errors corrected, including correction to wrong proofs. That is life.

So, let us continue. We take what is good. We correct what can be corrected. And we move forward.

I was going to make this post about the case m>1 – to give the complete algorithm. But somehow I got deflected, I stumbled. I will recover in the next post.

6 thoughts on “Giants that stumble”

  1. How is the quaternion p(t) produced.
    Given a vector m=(m1(t),m2(t),m3(t)), from Standing on the shoulders of giants -Reboot
    cos(theta)=m3(t) and sin(theta)=mp(t) assuming you are rotating the vector up to the z axis. The Rodrigues formula is simple to apply.
    But quaternions use 1/2 angles so would need \cos \left( 1/2\,\arccos \left( {\it m3} \left( t \right) \right) \right) etc. Don’t see how to convert it at the moment to produce the p(t) in the paper.

  4. Correct me is my intrepertation is incorrect.
    The quaternion is the encoding of the axis and angle of rotation. Changing the sign 0f say (W,X,-Y,0) to (W,-X,Y,0) changes the dirction vector by 180 degrees so the rotation rotates the vector in the opposite direction. If found(-W,-X,Y,0) rotates on opposite direction and works corectly. Have not tested othe x axis yet. Have just derived the formula on paper so far and was unclear which way the direction vector should be pointing. Some of these formulas are easier to simplify be hand than get Maple too.

Leave a Reply