Taking square root of space rotation

While taking squares is easy, taking square roots may be dangerous. Taking square root of -1 leads to imaginary numbers. Paul Dirac took square root of the Schrödinger equation, and he has created antimatter. Very dangerous. From Wikipedia:

The reaction of 1 kg of antimatter with 1 kg of matter would produce 1.8×1017 J (180 petajoules) of energy (by the mass–energy equivalence formula, E = mc2), or the rough equivalent of 43 megatons of TNT – slightly less than the yield of the 27,000 kg Tsar Bomb, the largest thermonuclear weapon ever detonated.

Now we want to take square roots of rotations. We will see how dangerous it is. We know that quaternions (equivalently elements of the special unitary group \mathrm{SU}(2) produce rotations. We even know the exact formula. Given and quaternion

    \[q= W+X\mathbf{i}+Y\mathbf{j}+Z\mathbf{k},\]

or, equivalently, U\in\mathrm{SU}(2)

(1)   \begin{equation*}U=W I+X\mathbf{i}+Y\mathbf{j}+Z\mathbf{k}=\begin{bmatrix}W+iZ&iX-Y\\iX+Y&W-iZ}\end{bmatrix},\end{equation*}

we have the rotation matrix R(U)=R(q) given as (see Putting a spin on mistakes):

(2)   \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*}

The matrix R(q) is determined by products of the quaternion components. Therefore R(q)=R(-q). Taking squares is easy. But how to go from R to q? This amounts to taking square root of the rotation. How do we do it? The problem can be solved, for instance, with a computer algebra system. Here is such a solution:

Given matrix R\in \mathrm{SO}(3) define

(3)   \begin{equation*} s=2\sqrt{R_{11}+R_{22}+R_{33}+1}.\end{equation*}


(4)   \begin{align*} W&=s/4,\\ X&=(R_{32}-R_{23})/s.\\ Y&=(R_{13}-R_{31})/s.\\ Z&=(R_{21}-R_{12})/s. \end{align*}

Then, with q= W+X\mathbf{i}+Y\mathbf{j}+Z\mathbf{k}, we have R(q)=R.
The following REDUCE code verifies this statement:

FOR I:=1:3 DO
FOR J:=1:3 DO
LET X**2=1-Y^2-Z^2-W^2;

At first sight it looks innocent. We have verified that we have a solution. What else do we need?
Let’s do a test. We take rotations about z-axis (see Mathematics and sex).


    \[R_0(t)=\exp tW(\vec{k})=\begin{bmatrix}\cos t&-\sin t&0\\\sin t&\cos t&0\\0&0&1\end{bmatrix}\]

and tilt by \pi/6 angle about y axis \vec{n}=(0,1,0), to obtain

(5)   \begin{equation*}R(t)=\begin{bmatrix}  \frac{1}{2} \sqrt{3} \cos (t) & -\frac{1}{2} \sqrt{3} \sin (t) & \frac{1}{2} \\  \sin (t) & \cos (t) & 0 \\  -\frac{\cos (t)}{2} & \frac{\sin (t)}{2} & \frac{\sqrt{3}}{2} \\ \end{bmatrix}.\end{equation*}

We now use the above map from rotation matrices to quaternions, and then stereographic projection. Mathematica draws the following path:

But we know from More recollections that it should look as follows:

Something is evidently going wrong at the bottom pf the plot.
What is going wrong is this: in the formula (3) we divide by s=2\sqrt{R_{11}+R_{22}+R_{33}+1}. When s is very close to 0 the calculations become unstable. That is what happens at the bottom of the graph. Mathematica skips the segment for which the precision becomes questionable. Instead of the nice arc, we see a straight line segment there.

In the next post we will be drawing the trajectory of the asymmetric top that is rotating and nutating close to its z axis. In order to avoid instabilities as the one above, we will take special precautions.

More recollections

The recollections in the last post were incomplete. So, here is the continuation.
We have derived the property

(1)   \begin{equation*}R\,W(\vec{v})\,R^t = W(R\vec{v}).\end{equation*}

On the other hand we stressed the fact that

given a unit vector \boldsymbol{\omega}:

(2)   \begin{equation*} Q(\boldsymbol{\omega},\theta)=\exp(\theta W(\boldsymbol{\omega}))\end{equation*}

is the rotation matrix that describes the rotation about the axis in the direction \boldsymbol{\omega} by an angle \theta.

From these two properties we can deduce the third one:

If R\in\mathrm{SO}(3) is a rotation, then rotation about the vector R\boldsymbol{\omega} by the angle \theta is given by the matrix

(3)   \begin{equation*} Q(R\boldsymbol{\omega},\theta)=R\, Q(\boldsymbol{\omega},\theta)\,R^t.\end{equation*}

To deduce this last property we use the fact that for any X and any invertible R we have:

    \[ e^{RXR^{-1}}=Re^XR^{-1}.\]

The proof of the above follows by expanding e^X as a power series, and noticing that (RXR^{-1})^n=RX^nR^{-1}.
I will use spin equivalent of Eq. (3) when explaining how exactly I got the two images from the last post that I am showing again below:

The advantage in using Eq. (3) is that it may be easy to calculate the exponential \exp(\theta W(\boldsymbol{\omega})), but not so easy the exponential \exp(\theta W(R\boldsymbol{\omega})).

Our aim is to draw trajectories in the rotation group \mathrm{SO}(3) made by the asymmetric top spinning about its largest moment of inertia axis, which in our setting is the third axis. As the rotation group \mathrm{SO}(3) is not very graphics friendly, we are using the double covering group \mathrm{SU}(2) or, equivalently, the group of unit quaternions. We then use stereographic projection to project trajectories from the 3-dimensional sphere S^3 to 3-dimensional Euclidean space \mathbf{R}^3.

But first things first. Before doing anything more we first have to translate the above properties into the language of \mathrm{SU}(2). Therefore another recollection is necessary at this point. In Pauli, rotations and quaternions we have defined three spin matrices:

(4)   \begin{equation*}s_1=\begin{bmatrix}0&1\\1&0\end{bmatrix},\quad s_2=\begin{bmatrix}0&i\\-i&0\end{bmatrix},\quad s_3=\begin{bmatrix}1&0\\0&-1\end{bmatrix}.\end{equation*}

In Putting a spin on mistakes we have explained that we have the map U\mapsto R(U) that associates rotation matrix R(U)\in\mathrm{SO}(3) to every matrix U\in\mathrm{SU}(2) in such a way that for every vector \vec{v} we have

(5)   \begin{equation*}U\,\vec{v}\cdot\vec{s}\,U^*=(R\vec{v})\cdot\vec{s}.\end{equation*}

It follows immediately from the above formula that U\mapsto R(U) is a group homomorphism, that is we have

(6)   \begin{equation*}R(UU')=R(U)R(U'),\quad R(U^*)=R(U)^t.\end{equation*}

At the end of Putting a spin on mistakes we have arrived at the formula:

(7)   \begin{equation*}\exp(\theta W(\vec{k}))=R\left(\exp(i\frac{\theta}{2}\vec{k}\cdot \vec{s})\right)\end{equation*}

which tells us that the matrix U(\vec{k},\theta) defined as

(8)   \begin{equation*}U(\vec{k},\theta)\stackrel{df}{=}\exp(i\frac{\theta}{2}\vec{k}\cdot \vec{s})\end{equation*}

describes the rotation about the direction of the unit vector \vec{k} by the angle \theta at the \mathrm{SU}(2) level.
The map U\mapsto R(U) is 2:1. We have R(U)=R(-U). Matrices U and -U implement the same rotation. When \theta changes from 0 to 2 \pi the matrix R in Eq. (7) describes the full 2\pi turn, but the matrix U becomes -I rather than I.. We need to make 4 \pi rotation at the \mathbf{R}^3 vector level to have 2\pi rotation at the spin level. That is described by \frac{\theta}{2} in the exponential in Eq. (7).
At the spin level we then obtain the following analogue of Eq. (\ref:rqr}):

(9)   \begin{equation*}V U(\vec{k},\theta)\,V^*=U(R(V)\vec{k},\theta),\quad V\in\mathrm{SU}(2).\end{equation*}

Now we can finally look back at the images of trajectories. In Seeing spin like an artist we noticed that rotation about the z-axis is described by the matrix

(10)   \begin{equation*}U(t)=\begin{bmatrix}\cos t/2+i\sin t/2&0\\0&\cos t/2 -i\sin t/2\end{bmatrix}=\begin{bmatrix}e^{i t/2}&0\\0&e^{-i t/2}\end{bmatrix},\end{equation*}

while rotation by the angle \phi about y-axis is described by the matrix

(11)   \begin{equation*}V(\phi)=\begin{bmatrix} \cos \left(\frac{\phi}{2}\right) & -\sin \left(\frac{\phi}{2}\right) \\ \sin \left(\frac{\phi}{2}\right) & \cos \left(\frac{\phi}{2}\right)\end{bmatrix}.\end{equation*}

U(t) describes uniform rotation of the top about its z-axis, when the z-axis of the top and of the laboratory frame are aligned. As observed in Towards the road less traveled with spin that is not an interested path for drawing. To make the path more interesting we tilt the laboratory frame. To this end we use, for instance, the matrix V(\phi) with \phi=\pi/6

(12)   \begin{equation*}V(\pi/6)=\frac{1}{2\sqrt{2}}\begin{bmatrix}1+\sqrt{3}&1-\sqrt{3}\\ \sqrt{3}-1&\sqrt{3}+1\end{bmatrix}.\end{equation*}

Now the axis of rotation of the top is tilted by 30 degrees with respect to the z-axis of the laboratory frame. We get new path of the evolution in the group \mathrm{SU}(2)

(13)   \begin{equation*} U'(t)=V(\phi/6)U(t)=\frac{1}{2\sqrt{2}}\begin{bmatrix}(\sqrt{3}+1)e^{ît/2}&(1-\sqrt{3})e^{-it/2}\\(\sqrt{3}-1)e^{it/2}&(\sqrt{3}+1)e^{-it/2}.\end{bmatrix}\end{equation*}

From Eq. (1) in Chiromancy in the rotation group we can now calculate real parameters W,X,Y,Z. For a general matrix U\in\mathrm{SU}(2) the formulas are:

(14)   \begin{align*} W&=\frac{1}{2}(U_{11}+U_{22}),\\  X&=\frac{1}{2i}(U_{12}+U_{21}),\\  Y&=\frac{1}{2}(U_{21}-U_{12}),\\  Z&=\frac{1}{2i}(U_{11}-U_{22}). \end{align*}

In our particular case we get:

(15)   \begin{align*} W&=\frac{1+\sqrt{3}}{2\sqrt{2}}\cos t/2,\\  X&=\frac{\sqrt{3}-1}{2\sqrt{2}}\sin t/2,\\  Y&=\frac{\sqrt{3}-1}{2\sqrt{2}}\cos t/2,\\  Z&=\frac{1+\sqrt{3}}{2\sqrt{2}}\sin t/2. \end{align*}

We then get rather awfully looking the parametric formulas for the stereographic projection:

(16)   \begin{align*} x(t)&=-\frac{\sqrt{2} \left(\sqrt{3}-1\right) \sin \left(\frac{t}{2}\right)}{\left(\sqrt{2}+\sqrt{6}\right) \cos \left(\frac{t}{2}\right)-4},\\ y(t)&=-\frac{\sqrt{2} \left(\sqrt{3}-1\right) \sin \left(\frac{t}{2}\right)}{\left(\sqrt{2}+\sqrt{6}\right) \cos \left(\frac{t}{2}\right)-4},\\ z(t)&=-\frac{\sqrt{2} \left(\sqrt{3}-1\right) \sin \left(\frac{t}{2}\right)}{\left(\sqrt{2}+\sqrt{6}\right) \cos \left(\frac{t}{2}\right)-4}. \end{align*}

The result (running t from 0 to r\pi) is not very interesting – just a circle:

It is better than the straight line with a non tilted top, but not a big deal. But the next thing we want to do is to rotate the tilt axis in the (x,y) plane, that is about the z-axis of the laboratory. If we want to rotate by an angle \psi, we should replace V(\phi) by U(\psi)V(\phi)U(\psi)^* – according to Eq. (9). That is we should look at the trajectory:


Now, U^*(\psi)U(t)=U(t-\psi). Since anyway we are are going to run with t through the whole (0,4\pi) interval, we can as well use U(t) instead of U(t-\psi). Therefore we can draw


We can calculate W,X,Y,Z and then x,y,z. Here are the results from my Mathematica code:

(17)   \begin{align*} W&=\cos \left(\frac{\phi }{2}\right) \cos \left(\frac{t+\psi }{2}\right),\\ X&=-\sin \left(\frac{\phi }{2}\right) \sin \left(\frac{\psi -t}{2}\right),\\ Y&=\sin \left(\frac{\phi }{2}\right) \cos \left(\frac{\psi -t}{2}\right),\\ Z&=\sin \left(\frac{\phi }{2}\right) \cos \left(\frac{\psi -t}{2}\right). \end{align*}

(18)   \begin{align*} x(\psi,\phi,t)&=\frac{\sin \left(\frac{\phi }{2}\right) \sin \left(\frac{\psi -t}{2}\right)}{\cos \left(\frac{\phi }{2}\right) \cos \left(\frac{t+\psi }{2}\right)-1},\\ y(\psi,\phi,t)&=\frac{\sin \left(\frac{\phi }{2}\right) \cos \left(\frac{\psi -t}{2}\right)}{1-\cos \left(\frac{\phi }{2}\right) \cos \left(\frac{t+\psi }{2}\right)},\\ z(\psi,\phi,t)&=\frac{\cos \left(\frac{\phi }{2}\right) \sin \left(\frac{t+\psi }{2}\right)}{1-\cos \left(\frac{\phi }{2}\right) \cos \left(\frac{t+\psi }{2}\right)}. \end{align*}

And here is the resulting of plotting 30 such circles, each for \phi=\pi/6, with \psi increasing from 0 t0 180 degrees, every 6 degrees:

It looks like half of a torus. With \phi=\pi/3, we get another, larger torus …
In the next posts we will disturb the motion of the top to see what kind of trajectories we will be getting then.

Mathematics and sex

Nowadays mathematics is almost everywhere. To know some cool mathematics (or mathematician) is nowadays sexy.

Therefore I am helping you here, on this blog, to have better relations in your lives. Just looking at the equations can already make you feel better. When you look for a little at Euler’s equations, you will sleep better. That is why today I recall some of these pretty formulas from the dynamics of spinning bodies.

Repetitio est mater studiorum – repetition is the mother of study/learning, that is what they say. So let us repeat, recollect, before traveling another road. If we climb high enough, then going down will be easy.

It is not that climbing up is difficult. It is not. We can take stairs. But it takeS time. And that is what we are going to do now: take time, do it easy way. Enjoying happy sliding down afterwards.

The two essential ingredients af the asymmetric top spinning in zero gravity are:

  1. Euler’s equations
  2. Attitude matrix equations

Euler equations, written in terms of the angular momentum vector \mathbf{L}, with components (L_1,L_2,L_3) are

(1)   \begin{eqnarray*} \frac{dL_1(t)}{dt}&=&(\frac{1}{I_3}-\frac{1}{I_2})L_2(t)L_3(t),\\ \frac{dL_2(t)}{dt}&=&(\frac{1}{I_1}-\frac{1}{I_3})L_3(t)L_1(t),\\ \frac{dL_3(t)}{dt}&=&(\frac{1}{I_2}-\frac{1}{I_1})L_1(t)L_2(t). \end{eqnarray*}

I_1,I_2,I_3 are principal moments of inertia. We assume that our spinning body is asymmetric, therefore I_1,I_2,I_3 are all different from each other, and we will order them as I_1<I_2<I_3. Perhaps it is worthwhile to mention that physics of ordinary materials requires that all three numbers are strictly positive, and that I_1+I_2\geq I_3.  I have my pet rigid body that is essentially flat, with I_1=1,I_2=2,I_3=3. It looks as on the picture below

It is kind of simple but it shows all the nontrivial behavior that a rigid body can have.

Remark: As noted by Bjab in the discussion of this post, the sentence above is incorrect. At the other end of the spectrum of rigid body shapes there are one-dimensional objects. These are flying rods with I_1=0,I_2=I_3\neq 0. They deserve a separate study.

But that is not important now.  What is important is that the components (L_1,L_2,L_3) are the components of the angular momentum vector with respect to the noninertial frame that rotates with the body. Its components with respect to the inertial laboratory frame are constant in time. This is “conservation of angular momentum”- one of the most fundamental “laws” of classical mechanics. Why such a law holds, or at least approximately holds, in our Universe is not known. Physicists and philosophers are debating about “the origin of inertia”.  Some say that “because of space”, some say “because of distant parallelism”, some other say “because of distant matter”. We will not worry about these problems now. We have our laboratory inertial frame, we have frame rotating with the body, and we have rotation matrix Q(t) that maps coordinates in the rotating frame to coordinates in the laboratory frame.

Recall from Angular momentum

Let \mathbf{E}_1(t), \mathbf{E}_2(t), \mathbf{E}_3(t) be an orthonormal frame corotating with the body, and aligned with its principal axes, and let \mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3 be an inertial laboratory frame, both centered at the center of mass of the body. The two frames are related by time-dependent orthogonal matrix A(t)

(2)   \begin{equation*} \mathbf{E}_i=A_{ji}\mathbf{e}_j,\quad X_i= A_{ij}x_j.\end{equation*}

The inverse of A is denoted by Q

(3)   \begin{equation*} Q= A^{-1}=A^t,\end{equation*}

and it is often called the attitude matrix. For a rotating body, if X_i are coordinates of a fixed point in the body, then its coordinates in the laboratory system change in time:

(4)   \begin{equation*} {\bf x}(t)=Q(t){\bf X}.\end{equation*}

Matrix Q(t) satisfies differential equation that is essentially nothing more than the definition of the angular velocity vector:

(5)   \begin{equation*} Q(t)^{-1}\frac{dQ(t)}{dt}=W(\boldsymbol{\Omega}(t)),\end{equation*}


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

where for any vector \vec{v}

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

We have

(8)   \begin{equation*}L_1=I_1\,\Omega_1,\quad L_2=I_2\,\Omega_2,\quad L_3=I_3\,\Omega_3.\end{equation*}

A couple of comments are due at this point. First of all a careful Reader will notice that in Angular momentum we had

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

while now we have

(10)   \begin{equation*}\frac{dQ(t)}{dt}=Q(t)W(\boldsymbol{\Omega}(t)).\end{equation*}

What is going on?
Several things at once. First we have the map \vec{v}\mapsto W(\vec{v}) that associates a matrix to every vector. In fact it associates an antisymmetric matrix. The association has the property that is easy to verify: for every vector \vec{w} we have


Acting with the matrix W(\vec{v}) is the same as taking cross product with vector \vec{v}. That such an association should exist should be not a surprise. Taking cross-product with a fixed vector is a linear operation, and every linear operation on vectors is implemented by a matrix. However many students who learn about cross products,

do not learn about their relation to antisymmetric matrices and bivectors. And many students that are learning about matrices, operations with them, eigenvalues and eigenvectors, do not learn about cross-product of vectors. Perhaps because the cross-product is particular to 3D?
Anyway, we have this association, and this association has a very nice property of “covariance”: for every orthogonal matrix R of determinant one we have:

(11)   \begin{equation*}R\,W(\vec{v})\,R^t = W(R\vec{v}).\end{equation*}

That is a very important and very useful property. It can be equivalently expressed as the “covariance of the cross-product”

(12)   \begin{equation*}(R\vec{w})\times (R\vec{v})=R(\vec{w}\times\vec{v}).\end{equation*}

But equivalent expression does not replace the proof, and to prove it needs a bunch of simple but lengthy algebraic calculations and taking into account the definitions of cross-product and determinant. It can be easily done with any computer algebra software. I skip it now.
The second thing is the difference between \boldsymbol{\omega} that represents the angular momentum vector in the laboratory frame, and \boldsymbol{\Omega} that represents the same vector in the body frame. Our matrix Q, by definition, connects the two representations (see Eq. (4) above, but now we skip the time-dependence):


Now, going back to Eq. (9), we have (again we skip the time dependence)


which solves our puzzle.
: Above we have used \boldsymbol{\Omega} and \boldsymbol{\omega} to distinguish between numerical representation of the same vector with respect to two different frames. But when there is no possibility of a confusion, we can denote a vector by any convenient letter whatsoever.
Let us see how the above recollection of facts can be applied. In Towards the road less traveled with spin there was the following statement:

The vertical z-axis is the natural axis to try to spin the thing. Imagine our top is floating in space, in zero gravity. We take the z-axis between our fingers, and spin the device. If our hand is not shaking too much, our top will nicely spin about the z-axis. This is the most stable axis for spinning.

The corresponding solution of Euler’s equations is \vec{\omega}=(0,0,\omega_3), where \omega_3 is a constant. The solution of the attitude matrix equation is Q(t)=\exp(W(t\vec{\omega})).

I promised to answer the question why is it so.

In general, when we have matrix differential equation of the type:


where X is a constant (i.e. time independent) matrix, we can verify that M(t)=\exp(tX) is a solution. In fact, it is the unique solution with the initial data M(0)=I. To prove it, one would have to expand e^{tX} into power series, take care about the convergence of the infinite series (no care needed, it is absolutely convergent), and then differentiate term by term. Let us take it for granted that is the case. In the quote from the other post above I used the letters \vec{\omega}=(0,0,\omega_3), instead of “more adequate”, \vec{\Omega}=(0,0,\Omega_3), but, as in the remark above, we are supposed to adjust the meaning to the context. For the matrix X in the attitude equation we have W(\vec{\omega}), so the solution with the property W(0)=I is Q(t)=\exp(tW(\vec{\omega})), which, because \vec{v}\mapsto W(\vec{v}) is a linear map, is the same as Q(t)=\exp(W(t\vec{\omega})).
In fact here we have the opportunity to make another application of today’s recollections. Suppose that our spinning top is completely symmetric – a perfect ball spinning one of its axes. Perfect ball means I_1=I_2=I_3. If so, the right side of the Euler’s equations is automatically zero. Therefore any constant vector \boldsymbol{\Omega} is a solution. Let us choose \boldsymbol{\Omega} a unit vector, so that the absolute value of the angular velocity is 1. That means our top makes a complete rotation about the axis along the vector \boldsymbol{\Omega} every 2\pi units of time. As noticed above the solution of the attitude equation is \exp(tW(\omega})). Therefore \exp(tW(\boldsymbol{\omega})) is the rotation matrix that describes the rotation about the axis in the direction \boldsymbol{\omega} by an angle t. We have already used this fact before, but now it was a good place for recollecting.
After all these recollections we are now ready to look again at the roads that are less traveled.

Of course it is not Albert Einstein who is the author, but, in the age of fake news, who cares?