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
q(R)=

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

and

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

OPERATOR R;
MATRIX RR(3,3);
FOR I:=1:3 DO
FOR J:=1:3 DO
RR(I,J):=R(I,J);
R(1,1):=W^2+X^2-Y^2-Z^2;
R(1,2):=2*(X*Y-W*Z);
R(1,3):=2*(W*Y+X*Z);
R(2,1):=2*(W*Z+X*Y);
R(2,2):=W^2-X^2+Y^2-Z^2;
R(2,3):=2*(Y*Z-W*X);
R(3,1):=2*(X*Z-W*Y);
R(3,2):=2*(W*X+Y*Z);
R(3,3):=W^2-X^2-Y^2+Z^2;
LET X**2=1-Y^2-Z^2-W^2;
%LET ABS(W)=W;
s:=2*SQRT(TRACE(RR)+1);
q0:=s/4;
q1:=(R(3,2)-R(2,3))/s;
q2:=(R(1,3)-R(3,1))/s;
q3:=(R(2,1)-R(1,2))/s;
END;

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

    \[\vec{k}=(0,0,1),\]

    \[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.

Leave a Reply