Category: Attitude matrix

Circles of eternal return

Posted on by arkajad

Let me recall the particular trajectory that we are discussing. In fact the particular part of this particular trajectory where the flip occurs.


Click on the image to open gif animation. It will take time to load as it is 2.5 MB.

The flip is much like one of these flips that Russian cosmonaut Dzhanibekov observed in zero gravity, and was very much perplexed by the phenomenon.

We are investigating the jungle of math that we have to apply in order to simulate such a behavior of an ideal Platonic rigid body. At present we are analyzing the special conditions when only one flip happens. Apart of a rather short (as compared to eternity) period of time when the flip occurs, the body rotates uniformly about its middle axis. Its trajectory in the rotation group, represented by quaternions of unit norm, follows a circle. One circle before flip, another circle after the flip.

In this note I am taking a close look at these two circles. Can we find their exact coordinate representation in 3D after stereographic projection?
Yes, we can, and we will do it now.
We use the formulas from Meeting with remarkable circles.
For the body of our choice, with moments of inertia I_1=1, I_2=2,I_3=3, and on the trajectory where d=1/I_2, the constants that enter the solution have values

(1)   \begin{eqnarray*} A_1&=&\sqrt{\frac{I_1(I_3-I_2)}{I_2(I_3-I_1)}}=1/2,\\ A_2&=&1,\\ A_3&=&\sqrt{\frac{I_3(I_2-I_1)}{I_2(I_3-I_1)}}=\frac{\sqrt{3}}{2},\\ B&=&\frac{1}{I_2}\,\sqrt{\frac{(I_2-I_1)(I_3-I_2)}{I_1I_3}}=\frac{1}{2\sqrt{3}},\\ \delta&=&\frac{\sqrt{I_2(I_3-I_1)}-\sqrt{I_1(I_3-I_2)}}{\sqrt{I_3(I_2-I_1)}}=\frac{1}{\sqrt{3}}. \end{eqnarray*}

In the solution of the Euler’s equations we have functions \tanh (Bt) and \mathrm{sech} (Bt)

(2)   \begin{eqnarray*} l_1(t)&=& A_1\,\mathrm{sech }(B t),\\ l_2(t)&=&A_2\,\tanh( B t),\\ l_3(t)&=&A_3\,\mathrm{sech }(B t), \end{eqnarray*}

with \mathrm{sech }(x)=1/\cosh(x).
Here are the plots of these functions:


For |t|>100 we have l_1(t) smaller than 10^{-12}, and for |t|>200 smaller than 10^{-25}. We can consider it being zero for all practical purposes.
As for l_2(t), it becomes practically constant and equal 1 for t>100, and -1 for t<-100.
We then have the function \psi(t)

(3)   \begin{equation*} \psi(t)=\frac{t}{I_2}+2\arctan\left(\delta \tanh(B t/2)\right). \end{equation*}

Again for all practical purposes for t>100 we have

(4)   \begin{equation*}\psi(t)\approx t/2+\pi/3,\end{equation*}

and for t<-100,

(5)   \begin{equation*}\psi(t)\approx t/2-\pi/3.\end{equation*}

Let us substitute these formulas into the solution for q(t)

(6)   \begin{eqnarray*} q_0(t)&=&\frac{\sqrt{1+l_1(t)}\cos\frac{\psi(t)}{2}}{\sqrt{2}},\\ q_1(t)&=&\frac{\sqrt{1+l_1(t)}\sin\frac{\psi(t)}{2}}{\sqrt{2}},\\ q_2(t)&=&\frac{l_3(t)\cos\frac{\psi(t)}{2}+l_2(t)\sin\frac{\psi(t)}{2}}{\sqrt{2(1+l_1(t))}},\\ q_3(t)&=&\frac{-l_2(t)\cos\frac{\psi(t)}{2}+l_3(t)\sin\frac{\psi(t)}{2}}{\sqrt{2(1+l_1(t))}}. \end{eqnarray*}

We obtain

(7)   \begin{eqnarray*} q_0(t)&=&\frac{\cos\frac{\psi(t)}{2}}{\sqrt{2}},\\ q_1(t)&=&\frac{\sin\frac{\psi(t)}{2}}{\sqrt{2}},\\ q_2(t)&=&\pm\frac{\sin\frac{\psi(t)}{2}}{\sqrt{2}},\\ q_3(t)&=&\frac{\mp\cos\frac{\psi(t)}{2}}{\sqrt{2}}. \end{eqnarray*}

The upper signs concern t>100 the lower signs concern t<-100.
In stereographic projection we are plotting (q_1,q_2,q_3) divided by 1-q_0. We see that the “future circle” is in the plane x=y, while the past circle in the orthogonal plane x=-y.

I did my geometric exercises and I have found that “the future circle” has origin at (0,0,-1), while “the past circle” has the origin at (0,0,1). Both circles have radius \sqrt{2}.

So, we have complete information about these “circles of eternal return”.
One more comment: in Another geodesic line I was showing what I thought was a subtle structure of the trajectory approaching the circle:

But now I know that it was a numerical artefact. As noticed in today’s post, for |t|>100 computer will see only one line, no fine structure. To be sure I asked my computer to take more points on the plot (100000 instead of 10000) and all the fine structure disappeared.

Quaternion evolution

Posted on by arkajad

In this note I will repeat some stuff that I was already explaining in the past, but we will also find some important new stuff.
The time evolution differential equation for the attitude matrix Q(t) is

(1)   \begin{equation*}\dot{Q}=QW(\boldsymbol{\Omega}),\end{equation*}

where \boldsymbol{\Omega}=\boldsymbol{\Omega}(t) is the angular velocity vector in the body frame, and W maps any vector \mathbf{v} into skew-symmetric matrix as follows:

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

It has the property that for any vector \mathbf{r} we have

(3)   \begin{equation*} W(\mathbf{v})\mathbf{r}=\mathbf{v}\times\mathbf{r}.\end{equation*}

Equation (1) follows immediately from the definition of the angular velocity. The angular velocity \boldsymbol{\omega} in the laboratory frame is defined by

(4)   \begin{equation*}W(\boldsymbol{\omega})=\dot{Q}Q^*.\end{equation*}

But Q maps coordinates of vectors in the body frame into their coordinates in the laboratory frame. Thus \boldsymbol{\omega}=Q\boldsymbol{\Omega}, and therefore

(5)   \begin{equation*}W(Q\boldsymbol{\Omega})=\dot{Q}Q^*.\end{equation*}

But from Eq. (3) we have W(Q\boldsymbol{\Omega})=QW(\boldsymbol{\Omega}) Q^*, therefore

    \[QW(\boldsymbol{\Omega}) Q^*=\dot{Q}Q^*,\]

and so Eq. (1) follows.
While description of rotations in terms of orthogonal 3\times3 matrices in principle suffices in classical mechanics, sometimes it is convenient to use the group of quaternions of unit norm, the group isomorphic to the group \mathrm{SU}(2) used in quantum mechanical description of half-integer spin particles. Quaternions have some advantages in numerical procedures (as for instance in 3D computer games), but they are also convenient for graphical representation of the intrinsic geometry of the rotation group. And this is what interests us, when we plot trajectories representing history of a spinning rigid body.


Quaternions of unit norm, representing rotations, form a 3-dimensional sphere in 4-dimensional Euclidean space, and we are projecting stereographically this sphere onto our familiar three-dimensional space, where we orient ourselves in a usual way known from everyday experience.
The question therefore arises: how the equation describing the time evolution looks like when represented in the quaternion setting?
To derive it we have to return to the fundamental relation between quaternions and rotations of vectors in space.
For every vector \mathbf{v}, with components (v_1,v_2,v_3) denote by \hat{\mathbf{v}} the pure imaginary quaternion defined as:

(6)   \begin{equation*}\hat{\mathbf{v}}=v_1\mathbf{i}+v_2\mathbf{j}+v_3\mathbf{k}.\end{equation*}

Then to unit quaternion q, that is such that qq^*=1, there corresponds rotation matrix R=R(q) such that for all \mathbf{v} the following identity holds

(7)   \begin{equation*}q\hat{\mathbf{v}}q^*=\widehat{R\mathbf{v}}.\end{equation*}

The second interesting property of the hat map is that for all \mathbf{v},\mathbf{w} we have

(8)   \begin{equation*}[\hat{\mathbf{v}},\hat{\mathbf{w}}]=2\widehat{\mathbf{v}\times\mathbf{w}},\end{equation*}

where [\cdot,\cdot] is the commutator. The property follows directly form the definitions and from the quaternion multiplication rules. Every pure imaginary quaternion is of the form \hat{\mathbf{v}} for some \mathbf{v}.

Derivation:
As shown on plaque by the Royal Canal at Broome Bridge in Dublin,

the essence of the quaternions is contained in the four equations

(9)   \begin{equation*}\mathbf{i}^2=\mathbf{j}^2=\mathbf{k}^2=-1,\quad \mathbf{i}\mathbf{j}\mathbf{k}=-1.\end{equation*}

From these, multiplying from the left and/or from the right by \mathbf{i} or \mathbf{j} or \mathbf{k} we derive

(10)   \begin{eqnarray*} \mathbf{i}\mathbf{j}&=&\mathbf{k},\quad \mathbf{j}\mathbf{i}=-\mathbf{k},\notag\\ \mathbf{j}\mathbf{k}&=&\mathbf{i},\quad \mathbf{k}\mathbf{j}=-\mathbf{i},\\ \mathbf{k}\mathbf{i}&=&\mathbf{j},\quad \mathbf{i}\mathbf{k}=-\mathbf{j}\notag. \end{eqnarray*}

From the definition of the hat map we have

(11)   \begin{eqnarray*}\hat{\mathbf{v}}&=&v_1\mathbf{i}+v_2\mathbf{j}+v_3\mathbf{k},\\ \hat{\mathbf{w}}=w_1\mathbf{i}+w_2\mathbf{j}+w_3\mathbf{k}. \end{eqnarray*}

Therefore

    \[ [\hat{\mathbf{v}},\hat{\mathbf{w}}]=\hat{\mathbf{v}}\hat{\mathbf{w}}-\hat{\mathbf{w}},\hat{\mathbf{v}}\]

    \[=(v_1\mathbf{i}+v_2\mathbf{j}+v_3\mathbf{k})(w_1\mathbf{i}+w_2\mathbf{j}+w_3\mathbf{k})-(w_1\mathbf{i}+w_2\mathbf{j}+w_3\mathbf{k})(v_1\mathbf{i}+v_2\mathbf{j}+v_3\mathbf{k}).\]

When we expand the parenthesis and use the multiplication rules in Eqs. (9),(10), we find that the terms with squares cancel out, while the terms with products like \mathbf{i}\mathbf{j} etc. can be organized as follows

    \[ [\hat{\mathbf{v}},\hat{\mathbf{w}}]=2\left[(v_2w_3-w_2v_3)\mathbf{i}+(v_3w_1-w_3v_1)\mathbf{j}+(v_1w_2-w_2v_1)\mathbf{k}\right].\]

The coefficients in the parentheses on the right are now exactly the components of the cross product \mathbf{v}\times\mathbf{w}.
Thus

    \[[\hat{\mathbf{v}},\hat{\mathbf{w}}]=2\widehat{\mathbf{v}\times\mathbf{w}}.\]

Suppose now that q(t) is a trajectory in the space of unit quaternions, while Q(t) is the corresponding trajectory in the rotation group. Thus we have

(12)   \begin{equation*}q(t)\hat{\mathbf{v}}q(t)^*=\widehat{Q(t)\mathbf{v}},\end{equation*}

for all t and all \mathbf{v}. We now differentiate both sides with respect to t. On the left we use the standard product rule for differentiation, but we pay attention so as to preserve the order, because multiplication of quaternions is non commutative. On the right we enter with the differentiation under the “hat”, because it is a linear operation. We obtain:

(13)   \begin{equation*}\dot{q}\hat{\mathbf{v}}q^*+q\hat{\mathbf{v}}\dot{q}^*=\widehat{\dot{Q}\mathbf{v}}.\end{equation*}

Assume now that Q(t) is a solution of Eq. (1), so that \dot{Q}=QW(\boldsymbol{\Omega}). We have

(14)   \begin{equation*}\dot{q}\hat{\mathbf{v}}q^*+q\hat{\mathbf{v}}\dot{q}^*=\widehat{(QW(\boldsymbol{\Omega})\mathbf{v})}.\end{equation*}

We now use Eq. (12) on the right to obtain

(15)   \begin{equation*}\dot{q}\hat{\mathbf{v}}q^*+q\hat{\mathbf{v}}\dot{q}^*=q\widehat{(W(\boldsymbol{\Omega})\mathbf{v})}q^*.\end{equation*}

In order to get rid of \dot{q}^* we differentiate the defining equation of unit quaternions qq^*=1 to obtain

(16)   \begin{equation*}\dot{q}q^*+q\dot{q}^*=0,\end{equation*}

or

(17)   \begin{equation*}\dot{q}^*=-q^*\dot{q}q^*.\end{equation*}

Putting this into Eq. (15) and multiplying both sides on the right by q we obtain

(18)   \begin{equation*}\dot{q}\hat{\mathbf{v}}-q\hat{\mathbf{v}}q^*\dot{q}=q\widehat{(W(\boldsymbol{\Omega})\mathbf{v})}.\end{equation*}

Multiplying from the left by q^*:

(19)   \begin{equation*}q^*\dot{q}\hat{\mathbf{v}}-\hat{\mathbf{v}}q^*\dot{q}=\widehat{(W(\boldsymbol{\Omega})\mathbf{v})}.\end{equation*}

From Eq. (16) it follows that the quaternion q^*\dot{q} is pure imaginary: (q^*\dot{q})^*=-q^*\dot{q}. Therefore there exists vector \mathbf{w} such that

(20)   \begin{equation*}q^*\dot{q}=\hat{\mathbf{w}}.\end{equation*}

Eq. (19) can now be written as

(21)   \begin{equation*}[\hat{\mathbf{w}},\hat{\mathbf{v}}]=\widehat{\boldsymbol{\Omega}\times\mathbf{v}},\end{equation*}

where on the right we have used Eq. (3). Applying Eq. (8) to the left we end with

(22)   \begin{equation*}2\widehat{\mathbf{w}\times\mathbf{v}}=\widehat{\boldsymbol{\Omega}\times\mathbf{v}}.\end{equation*}

Since the above holds for any \mathbf{v}, we deduce that

(23)   \begin{equation*}\hat{\mathbf{w}}=\frac{1}{2}\hat{\boldsymbol{\Omega}}.\end{equation*}

Using (20) we finally obtain:

Evolution equation for quaternions

(24)   \begin{equation*} \dot{q}=\frac{1}{2}q\hat{\boldsymbol{\Omega}}.\end{equation*}

Note: It may be that the final formula can be derived in a shorter way, but I do not know how. It is this last formula that I was using when verifying that the algorithm provided in Meeting with remarkable circles gives indeed a solution of the evolution equation.

New wine in new bottles

Posted on by arkajad

In my recent post I was presenting new formulas for the attitude matrix, different from those of a month ago. And I did not explain in detail how I got them, and why they look like being in disagreement with the previous ones. This caused a legitimate concern from Bjab, when he noticed the apparent disagreement. And indeed putting new wine into old bottles would lead to a trouble.

Mar 2:21 No man also seweth a piece of new cloth on an old garment: else the new piece that filled it up taketh away from the old, and the rent is made worse.
Mar 2:22 And no man putteth new wine into old bottles: else the new wine doth burst the bottles, and the wine is spilled, and the bottles will be marred: but new wine must be put into new bottles.


What we did is this: we poured new wine into new bottles. And here is how it is done.

We consider free rigid body with normalized (unit length) angular momentum vector, and the doubled kinetic energy equal to d=1/I_2. Euler’s equations can be then solved explicitly and we choose the following solution:

(1)   \begin{eqnarray*} A_1&=&\sqrt{\frac{I_1(I_3-I_2)}{I_2(I_3-I_1)}}=1/2,\\ A_2&=&1,\\ A_3&=&\sqrt{\frac{I_3(I_2-I_1)}{I_2(I_3-I_1)}}=\frac{\sqrt{3}}{2},\\ B&=&\frac{1}{I_2}\,\sqrt{\frac{(I_2-I_1)(I_3-I_2)}{I_1I_3}}=\frac{1}{2\sqrt{3}},\\ \end{eqnarray*}

(2)   \begin{eqnarray*} L_1(t)&=& A_1\,\mathrm{sech }(B t),\\ L_2(t)&=&A_2\,\tanh( B t),\\ L_3(t)&=&A_3\,\mathrm{sech }(B t), \end{eqnarray*}

where \mathrm{sech }(x)=1/\cosh(x).
We will need the angular velocity vector with components

(3)   \begin{equation*}\Omega_i(t)=\frac{L_i(t)}{I_i}.\end{equation*}

We want to solve now the differential equation for the attitude matrix Q(t):

(4)   \begin{equation*}\frac{dQ(t)}{dt}=Q(t)W(t),\end{equation*}

where

(5)   \begin{equation*}W(t)=\begin{bmatrix}0&-\Omega_3(t)&\Omega_2(t)\\\Omega_3(t)&0&-\Omega_1(t)\\-\Omega_2(t)&\Omega_1(t)&0\end{bmatrix}.\end{equation*}

The strategy is this: We seek the solution of the form Q(t)=Q_2(t)Q_1(t), where Q_1(t) is a certain matrix cleverly constructed out of L_i(t), and Q_2(t) is a simple rotation matrix by some angle \psi(t).
There are two methods: old one, and new one.

The old method.

In the old method, as used for instance in Taming the T-handle continued (we use now Q_1,Q_2 instead of Q_0,Q_1 there), we set:

(6)   \begin{equation*} L_p(t)=\sqrt{L_1(t)^2+L_2(t)^2},\end{equation*}

(7)   \begin{equation*}Q_1(t)=\begin{bmatrix}\frac{L_1(t)L_3(t)}{L_p(t)}&\frac{L_2(t)L_3(t)}{L_p(t)}&-L_p(t)\\-\frac{L_2(t)}{L_p(t)}&\frac{L_1(t)}{L_p(t)}&0\\L_1(t)&L_2(t)&L_3(t)\end{bmatrix},\end{equation*}

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

Then we solve for \psi(t):

(9)   \begin{equation*}\psi(t)=\frac{t}{I_2}+\arctan\left( \sqrt{\frac{(I_2-I_1)I_3}{I_1(I_3-I_2)}}\tanh\sqrt{\frac{(I_2-I_1)(I_3-I_2)}{I_1I_2^2I_3}}t \right).\end{equation*}

The new method.

(10)   \begin{equation*} Q_1(t) =\begin{bmatrix} L_1(t)& L_2(t)& L_3(t)\\-L_2(t)& 1 - \frac{L_2(t)^2}{1 + L_1(t)}& -\frac{L_2(t) L_3(t)}{1 + L_1(t)}\\ -L_3(t)& -\frac{L_2(t) L_3(t)}{1 + L_1(t)}& 1 - \frac{L_3(t)^2}{1 + L_1(t)},\end{bmatrix} \end{equation*}

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

(12)   \begin{equation*} \psi(t)=\frac{t}{I_2}+2\arctan\left(\delta \tanh(B t/2)\right), \end{equation*}

where

(13)   \begin{equation*}\delta&=&\frac{\sqrt{I_2(I_3-I_1)}-\sqrt{I_1(I_3-I_2)}}{\sqrt{I_3(I_2-I_1)}}.\end{equation*}

Comparing old and new wines

In both cases we get a solution of the attitude differential equation (4). But these are different solutions. If Q(t) is a solution, and if R is a fixed rotation matrix of determinant one, then RQ(t) is another solution. This is the same body rotating, but observed from two different laboratory frames.
Let Q_{old}(t) and Q_{new}(t) be the two solutions, one obtained with the old, the other one with the new method.
With some little algebra it can be verified that Q_{new}(0)=RQ_{old}(0), where

(14)   \begin{equation*}R=\begin{bmatrix}0&0&1\\0&1&0\\-1&0&0\end{bmatrix}.\end{equation*}

.
Therefore for all t

(15)   \begin{equation*}Q_{new}(t)=R Q_{old}(t).\end{equation*}

I verified it numerically. The above comparison was for rotation matrices. For quaternions it is similar, though the quaternion case deserves a comment. I will give it in another note.

Update: Mathematica notebook verifying Eq. (15) numerically is available here.