Dzhanibekov Effect | Open System - Ark's blog

Standing on the shoulders of giants – Reboot

Posted on 2017/02/282017/03/17 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…

Introducing geodesics

Posted on 2017/02/262017/02/26 by arkajad

Why are geodesics important? Probably because they are very simple. They are, in a sense, the simplest possible paths. Quoting from Richard Buckminster Fuller Domes and archives, 1960, 1965

Fuller inspired by his observations of nature. The inventor applies the concept of the geodesic line (the shortest line joining two points on a surface) to construct the most balanced, lightweight and resistant structure possible. His domes are a synthesis of all of the inventor’s fundamental precepts, combining a reasoned and aesthetic use of technological progress with a holistic conception of man’s relationship to nature. Such was the reputation of the inventor in the scientific domain that a family of carbon-based molecules with a geodesic structure was named after him: Buckminster fullerenes, later changed to fullerenes. Many of these molecules have played a role in recent nanotechnology discoveries.

Quoting from Wikipedia – Geodesics on an ellipsoid :

Geodesics on an ellipsoid of revolution
There are several ways of defining geodesics (Hilbert & Cohn-Vossen 1952, pp. 220–221). A simple definition is as the shortest path between two points on a surface. However, it is frequently more useful to define them as paths with zero geodesic curvature—i.e., the analogue of straight lines on a curved surface. This definition encompasses geodesics traveling so far across the ellipsoid’s surface (somewhat more than half the circumference) that other distinct routes require less distance. Locally, these geodesics are still identical to the shortest distance between two points.

If the Earth is treated as a sphere, the geodesics are great circles (all of which are closed) and the problems reduce to ones in spherical trigonometry. However, Newton (1687) showed that the effect of the rotation of the Earth results in its resembling a slightly oblate ellipsoid and, in this case, the equator and the meridians are the only closed geodesics. Furthermore, the shortest path between two points on the equator does not necessarily run along the equator. Finally, if the ellipsoid is further perturbed to become a triaxial ellipsoid (with three distinct semi-axes), then only three geodesics are closed and one of these is unstable.
….
On a triaxial ellipsoid, there are only three simple closed geodesics, the three principal sections of the ellipsoid given by X = 0, Y = 0, and Z = 0.

I am not playing with geodesics on three axial ellipsoid. But studying unexpected flips of asymmetric spinning top in zero gravity brings us very close to this circle of ideas. We do have three axial ellipsoid, because we have three different moments of inertia $I_1<I_2<I_3.$ And we do have geodesics, even if we have to wait a little bit to see them nicely introduced in this series of posts. Histories of rotations in the rotation group are geodesics. We have seen a bunch of them in two previous posts. These were closed. Nothing particularly interesting – who did not see a circle? True, we have seen a bunch of circles spanning a torus, but these were generated artificially by a rotating observer.

But no it is time to move beyond the safe mode. If you have a spinning top in free space, and if you struck it, like Peggy Whitson

it will start to nutate. And when all three moments of inertia are different, this nutation is non-periodic. Geodesics have infinite length. They are certainly not the shortest connections between points in the group, in any sense. One geodesic line is wandering through the three-dimensional sphere $S^3$ sometimes almost returning to the starting point, then traveling far away. Strange are these trajectories.

Below is a part of one such trajectory, for $I_1=1,I_2=2,I_3=3$ and $d=0.4.$ I am rotating it in animation, so that you can have a better view of its 3D structure. Of course, as in previous posts, this is stereographic projection from $S^3.$

Geodesic on three sphere — Geodesic of the left-invariant metric in the group $\mathrm{SU}(2)$

Of course I will have to explain the details …

Taming the T-handle continued

Posted on 2017/02/132017/02/23 by arkajad

This is a simple continuation from the last post “Taming the T-handle“. We ended up with the equation

(1) $\begin{equation*}\frac{d\psi(t)}{dt}=c_1+\frac{c_2}{1+c_3\,\sn^2(Bt,m)},\end{equation*}$

where

(2) $\begin{equation*}c_1=\frac{L}{I_3},\quad c_2=L\left(\frac{1}{I_1}-\frac{1}{I_3}\right),\quad c_3=\frac{I_3(I_2-I_1)}{I_1(I_3-I_2)}.\end{equation*}$

To find $\psi(t),$ let us assume that $\psi(0)=0$ , we need to integrate

$\psi(t)=c_1 t+c_2\int_0^t \frac{1}{1+c_3\,\sn^2(Bs,m)}\,ds.$

Setting $Bt=u$ we transform the above into

$\psi(t)=c_1 t+\frac{c_2}{B}\int_0^{Bt} \frac{1}{1+c_3\,\sn^2(u,m)}\,du.$

Searching the net we can find the integral of this type, for instance on the Wolfram’s page about elliptic integrals of the third order we can find:

$\Pi(n;z|m)=\int_0^{F(z|m)}\,\frac{1}{1-n\,\sn(u|m)^2}\, du.$

Then we set $F(z|m)=Bt,$ therefore $z= \am(Bt,m)$ – see Jacobi amplitude- realism or cubism,
so that

$\psi(t)=c_1 t+\frac{c_2}{B}\,\Pi(-c_3;\am(Bt,m)|m).$

In principle we could now substitute the values of $c_1,c_2,c_3$ and consider our task essentially done. But that would be not very prudent! The point is that we have the parameter $m$ and, depending on the case, this parameter can have value: $<1,$ $1,$ and $>1.$ The case $m=1$ is very special and should be treated separately. No elliptic functions are needed in this case. Moreover not every software used for calculations and simulations will know what to do with $m>1,$ and even if it pretends to know, it may happen that in this domain the software is not sufficiently tested and debugged. It is for this reason that we now consider the three cases separately.

The case of $d<1/I_2,$ i.e. $m<1.$

This is the most straightforward case. From Taming the T-handle, Eq. (10), we have

(3) $\begin{equation*}B=L\sqrt{\frac{(1-dI_1)(I_3-I_2)}{I_1I_2I_3}},\end{equation*}$

(4) $\begin{equation*}m=\frac{(dI_3-1)(I_2-I_1)}{(1-dI_1)(I_3-I_2)}.\end{equation*}$

Thus

$\frac{c_2}{B}=(I_3-I_1)\sqrt{\frac{I_2}{I_1I_3(I_3-I_2)(1-dI_1)}},$

and we are done.

The case of $d=1/I_2,$ i.e. $m=1.$

This is a tricky case. For $m=1$ we have (see Jacobi Elliptic Functions, Eqs. (33)-(35), and this comment)

(5) $\begin{eqnarray*} \cn(u,1)&=&\frac{1}{\cosh u},\\ \sn(u,1)&=&\tanh u,\\ \dn(u,1)&=&\frac{1}{\cosh u}. \end{eqnarray*}$

To find $\psi(t)$ we need to integrate

$\psi(t)=c_1 t+c_2\int_0^t \frac{1}{1+c_3\,\tanh^2(Bs,m)}\,ds.$

We can go to Wolfram Alpha, and ask it to integrate it. It will, but we will not like the answer. It looks very-very ugly!

In principle we could try to simplify it, but it is better to make one step back in order to make two steps forward! Let us go back to Eq. (1) and write it as just one fraction:

(6) $\begin{equation*}\frac{d\psi(t)}{dt}=c_1+\frac{c_2}{1+c_3\,\tanh^2 Bt} =\frac{c_1+c_2+c_1c_3\tanh^2 Bt}{1+c_3\tanh^2 Bt}.\end{equation*}$

Now, we write $\tanh^2 u=\sinh^2 u /\cosh^2 u=(\cosh^2u -1)/\cosh^2 u,$ and get

(7) $\begin{equation*}\frac{d\psi(t)}{dt}=\frac{-c_1c_3+(c_1+c_2+c_1c_3)\cosh^2 Bt}{-c_3+(1+c_3)\cosh^2 Bt}.\end{equation*}$

So, we need to integrate the function

(8) $\begin{equation*}\frac{d\psi}{dt}=\frac{\alpha+\beta \cosh^2Bt}{\gamma+\delta\cosh^2 Bt}.\end{equation*}$

We can go again to Wolfram Alpha and what we get this time is so much nicer!

All we have to substitute the constants. It is a mechanical task, so we can use REDUCE, so we use REDUCE or any other software capable of symbolic computation. I used Mathematica. Here is the result, followed by a complete summary of the case:

(9) $\begin{equation*}\psi(t)=\frac{Lt}{I_2}+L\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*}$

(10) $\begin{eqnarray*}A_1&=&L\sqrt{\frac{I_1(I_3-I_2)}{I_2(I_3-I_1)}},\\ A_2&=&L,\\ A_3&=&L\sqrt{\frac{I_3(I_2-I_1)}{I_2(I_3-I_1)}},\\ B&=&\frac{L}{I_2}\sqrt{\frac{(I_2-I_1)(I_3-I_2)}{I_1 I_3}},\end{eqnarray*}$

(11) $\begin{eqnarray*} L_1(t)&=&\frac{A_1}{\cosh Bt},\\\ L_2(t)&=&A_2 \tanh Bt,\\\ L_3(t)&=&\frac{A_3}{\cosh Bt},\end{eqnarray*}$

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

(13) $\begin{equation*}Q_0(t)=\begin{bmatrix}\frac{L_1(t)L_3(t)}{LL_p(t)}&\frac{L_2(t)L_3(t)}{LL_p(t)}&-\frac{L_p(t)}{L}\\-\frac{L_2(t)}{L_p(t)}&\frac{L_1(t)}{L_p(t)}&0\\\frac{L_1(t)}{L}&\frac{L_2(t)}{L}&\frac{L_3(t)}{L}\end{bmatrix},\end{equation*}$

(14) $\begin{equation*}Q_1(t)=\begin{bmatrix}\cos\psi(t)&-\sin\psi(t)&0\\\sin\psi(t)&\cos\psi(t)&0\\0&0&1\end{bmatrix},\end{equation*}$

(15) $\begin{equation*}Q(t)=Q_1(t)Q_0(t),\end{equation*}$

Changing signs of any two of three components $(L_1,L_2,L_3)$ we can obtain altogether four different evolutions of a painted T_handle. But that, and also the case of $m>1$ must wait for the next post.

But here are four pictures of the painted T-handle, rotated according to the above algorithm, at $t=10,$ where the three last cases correspond to flipped sign of $(L_1,L_2,L_3)$ as follows: $(1,2),(2,3),(1,3)$

Update

You can download my experimental Mathematica cdf file that can be run using free CDF Player. It shows the code and the animation for the four cases mentioned above.