Category: Elliptic functions

Attitude matrix and quaternion path for m>1

Posted on by arkajad

I noticed that somehow I did not finish with the case m>1. So today, without further ado, I am posting the algorithm.
We have the body with I_1<I_2<I_3, and we are considering the case with d>1/I_2, where d is, as always, the ratio of the doubled kinetic energy to angular momentum squared.
Then we define

(1)   \begin{eqnarray*} A_1&=&\sqrt{\frac{I_1(dI_3-1)}{I_3-I_1}},\\ A_2&=&\sqrt{\frac{I_2 (1-dI_1)}{I_2-I_1}},\\ A_3&=&\sqrt{\frac{I_3 (1-dI_1)}{I_3-I_1}},\\ B&=&\sqrt{\frac{(dI_3-1)(I_2-I_1)}{I_1I_2I_3}},\\ \mu&=&\frac{1}{m}=\frac{(1-dI_1)(I_3-I_2)}{(dI_3-1)(I_2-I_1)}. \end{eqnarray*}

(2)   \begin{eqnarray*} L_1(t)&=&A_1\,\dn(Bt,\mu),\\ L_2(t)&=&A_2\,\sn(Bt,\mu),\\ L_3(t)&=&A_3\,\cn(Bt,\mu). \end{eqnarray*}

(3)   \begin{equation*} \alpha=\frac{I_3-I_1}{\sqrt{\frac{I_1(dI_3-1)(I_2-I_1)I_3}{I_2}}}, \end{equation*}

(4)   \begin{equation*} \nu=\frac{I_3-dI_1I_3}{I_1-dI_1I_3}. \end{equation*}

(5)   \begin{equation*} \psi(t)=\frac{t}I_3-\arctan\left((A_2/A_1)\mathrm{sd}(Bt,\mu) \right)+\alpha \Pi(\nu,\am(Bt,\mu),\mu),\ \end{equation*}

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

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

(8)   \begin{equation*} Q(t)=Q_2(t)Q_1(t). \end{equation*}

(9)   \begin{eqnarray*} q_0(t)&=&\sqrt{\frac{1+L_3(t)}{2}}\cos\frac{\psi(t)}{2},\\ q_1(t)&=&\frac{1}{\sqrt{2(1+L_3(t))}}\left(L_2(t)\cos\frac{\psi(t)}{2}+L_1(t)\sin\frac{\psi(t)}{2}\right),\\ q_2(t)&=&\frac{1}{\sqrt{2(1+L_3(t))}}\left(L_2(t)\sin\frac{\psi(t)}{2}-L_1(t)\cos\frac{\psi(t)}{2}\right),\\ q_3(t)&=&\sqrt{\frac{1+L_3(t)}{2}}\,\sin\frac{\psi(t)}{2},\\ q(t)&=&(q_0(t),q_1(t),q_2(t),q_3(t))=q_0(t)+\mathbf{i}\,q_1(t)+\mathbf{j}\,q_1(t)+\mathbf{k}\,q_3(t). \end{eqnarray*}

I use the above formulas to draw a stereographic projection of one particular path. So, I take I_1=1,I_2=2,I_3=3,d=0.5000001, and do the parametric plot of the curve \mathbf{r}(t) in \mathbf{R}^3,
with

    \[\mathbf{r}(t)=\left(\frac{q_1(t)}{1-q_0(t)},\frac{q_2(t)}{1-q_0(t)},\frac{q_3(t)}{1-q_0(t)}\right).\]

I show below two plots. One with t\in(-1000,1000), and one with t\in(-10000,10000). For this selected value of d, the time between consecutive flips, given by the formula

(10)   \begin{equation*}\tau =4\mathrm{EllipticK}(\mu)/B= 116.472.\end{equation*}

So, for t\in(-10000,10000) we have somewhat less than 200 flips, and the lines are getting rather densely packed in certain regions.
Notice that is just one geodesic line, geometrically speaking the straightest possible line in the geometry determined by the inertial properties of the body.

It is this “geometry”that will become the main subject of the future notes.

Geodesic line for t between -1000 and 1000

Geodesic line for t between -10000 and 10000

How can such a line be “straight”? Well, it is….

Attitude matrix for m<1

Posted on by arkajad

The last post was about The final answer for the Universe in which m>1. Yet, as I see it today, it was neither complete nor final. It was a very bad approximation. It may be also that our universe is also a very bad approximation to the one that has been intended. But that is another story. I think I will return to this problem in my blogging, but now I will concentrate on the algorithm for the rotation matrix. At the same time I will make a comment how the algorithm from the last post, for the quaternions, should be improved. The fact is that I was making wrong assumptions in my mind. Making wrong assumptions, without being conscious of the fact that one is making such assumptions, that is very dangerous. Lot of bad things happens around us for this reason….

Plane crash
The result of wrong assumptions

So, here is the code:

The case m<1.
The solution \mathbf{L}(t) of the Euler’s equations has two trajectories. For one with L_3(t)>0 we take

(1)   \begin{eqnarray*} A_1&=&\sqrt{\frac{I_1 (d I_3-1)}{I_3-I_1}},\\ A_2&=&\sqrt{\frac{I_2 (d I_3-1)}{I_3-I_2}},\\ A_3&=&\sqrt{\frac{I_3 (1-d I_1)}{I_3-I_1}},\\ \end{eqnarray*}

while for the other one, with L_3(t)<0 we take

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

Let

(3)   \begin{eqnarray*} B&=&\sqrt{\frac{(1-d I_1) (I_3-I_2)}{I_1 I_2 I_3}},\\ m&=&\frac{(d I_3-1) (I_2-I_1)}{(1-d I_1) (I_3-I_2)}. \end{eqnarray*}

The solution \mathbf{L}(t) of the Euler’s equations is given by

(4)   \begin{eqnarray*} L_1(t)&=&A_1\, \cn (Bt,m),\\ L_2(t)&=&A_2\, \sn (Bt,m),\\ L_3(t)&=&A_3\, \dn (Bt,m), \end{eqnarray*}

With constants \alpha and \nu defined as

(5)   \begin{equation*} \alpha=\frac{I_3-I_1}{\sqrt{\frac{I_1(1-d I_1)(I_3-I_2)I_3}{I_2}}}, \end{equation*}

(6)   \begin{equation*} \nu=\frac{I_1-dI_1I_3}{I_3-dI_1I_3} \end{equation*}

we set the phase variable \psi(t) as

(7)   \begin{equation*} \psi(t)=\frac{t}I_1+\arctan\left((A_2/A_3)\mathrm{sd}(Bt,m) \right)-\alpha \Pi(\nu,\am(Bt,m),m),\ \end{equation*}

where the Jacobi function \mathrm{sd} is defined as \mathrm{sd}(u,m)=\sn(u,m)/\dn(u,m).
Let

(8)   \begin{equation*} \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*}

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

Then the attitude matrix Q(t) is given by

(10)   \begin{equation*}Q(t)=Q_2(t)Q_1(t).\end{equation*}

If you compare this with the formulas I gave yesterday, you will find the following important difference: now I am playing with signs in the constants A_1,A_3 rather than with the signs in formulas for L_1(t),L_3(t). It took me whole day to realize that this way is a better way. I was not realizing before that there is one sign in the formula for \psi(t) that depends on the trajectory we are taking. The constant A_2/A_3 does this job now.

In a couple of days I will also change the formulas in the quaternionic algorithm to make them universal as well.

Now, I check the modified example3. The modified initial data in init_example3.dat file are as follows:
1.0000000000D0
10.0000000000D0
-0.544332842491675D0
0.729131780907662D0
-0.414811526666455D0
1.000000000000000D0
0.000000000000000D0
0.000000000000000D0
0.000000000000000D0
1.000000000000000D0
0.000000000000000D0
0.000000000000000D0
0.000000000000000D0
1.000000000000000D0
1.000000000000000D0
1.012686988782515D0
3.306237422473038D0
5
They are the same as in the quaternionic case discussed yesterday. The initial attitude matrix is the identity matrix. We get the same t_0. I use Mathematica to calculate the transposed final Q

Transpose[Inverse[Q[t0]].Q[t0 + 10.0]]
{{0.0656754, 0.550118, -0.832501}, {0.995487, 0.0211483,
0.0925082}, {0.0684964, -0.834819, -0.546246}}

We compare it with the Fortran output from 903 code:

0.0656754297
0.55011776
-0.832500564
0.995487311
0.0211482921
0.0925081752
0.0684963551
-0.834819262
-0.546246326

And it seems that our code is doing the same job as the Fortran code from the giants!
Yet, as I said at the beginning, it took me the whole days to figure it out why for some initial data I was in agreement, but for another initial data there was a disagreement.

There is one little plus of the code above. The coefficient in front of \arctan is fixed, and equal 1.! Zanna and Celledoni in their paper and in their code have complicated formulas for this coefficient. It would take some work to verify that their complicated formula simplifies to 1. I hope I am right here….

The final answer for the Universe in which m<1

Posted on by arkajad

This is not an independent post or note or whatever. It is a continuation. In fact it is the continuation of Approaching the ultimate answer for m<1. We have approached enough. We are ready.

There was a moment’s expectant pause while panels slowly came to life on the front of the console. Lights flashed on and off experimentally and settled down into a businesslike pattern. A soft low hum came from the communication channel.
“Good morning,” said Deep Thought at last.
“Er … good morning, O Deep Thought,” said Loonquawl nervously, “do you have … er, that is …”
“An answer for you?” interrupted Deep Thought majestically. “Yes. I have.”
The two men shivered with expectancy. Their waiting had not been in vain.
“There really is one?” breathed Phouchg.
“There really is one,” confirmed Deep Thought.

Douglas Adams, The Hitchhiker’s Guide to the Galaxy

We did not get to the answer yesterday, but today is the day. We will be using the formulas from the last post. Here they are again, repeated for your convenience:

The case m<1.

(1)   \begin{eqnarray*} A_1&=&\sqrt{\frac{I_1 (d I_3-1)}{I_3-I_1}},\\ A_2&=&\sqrt{\frac{I_2 (d I_3-1)}{I_3-I_2}},\\ A_3&=&\sqrt{\frac{I_3 (1-d I_1)}{I_3-I_1}},\\ B&=&\sqrt{\frac{(1-d I_1) (I_3-I_2)}{I_1 I_2 I_3}},\\ m&=&\frac{(d I_3-1) (I_2-I_1)}{(1-d I_1) (I_3-I_2)}. \end{eqnarray*}

The solution \mathbf{L}(t) of the Euler’s equations has two trajectories: one with L_3(t)>0 is given by

(2)   \begin{eqnarray*} L_1(t)&=&A_1\, \cn (Bt,m),\\ L_2(t)&=&A_2\, \sn (Bt,m),\\ L_3(t)&=&A_3\, \dn (Bt,m), \end{eqnarray*}

while the other one, with L_3(t)<0 is given by

(3)   \begin{eqnarray*} L_1(t)&=&-\,A_1\, \cn (Bt,m),\\ L_2(t)&=&A_2\, \sn (Bt,m),\\ L_3(t)&=&-\,A_3\, \dn (Bt,m). \end{eqnarray*}

With constants \alpha and \nu defined as

(4)   \begin{equation*} \alpha=\frac{I_3-I_1}{\sqrt{\frac{I_1(1-d I_1)(I_3-I_2)I_3}{I_2}}}, \end{equation*}

(5)   \begin{equation*} \nu=\frac{I_1-dI_1I_3}{I_3-dI_1I_3} \end{equation*}

we set the phase variable \psi(t) as

(6)   \begin{equation*} \psi(t)=\frac{t}I_1-\arctan\left((A_2/A_3)\mathrm{sd}(Bt,m) \right)-\alpha \Pi(\nu,\am(Bt,m),m),\ \end{equation*}

where the Jacobi function \mathrm{sd} is defined as \mathrm{sd}(u,m)=\sn(u,m)/\dn(u,m).
Then the quaternionic attitude solution is given by q(t)=W(t)+\mathbf{i}X(t)+\mathbf{j}Y(t)+\mathbf{k}Z(t), with

(7)   \begin{eqnarray*} W(t)&=&\cos \frac{\psi(t)}{2}\sqrt{\frac{1+L_1(t)}{2}},\\ X(t)&=&\sin \frac{\psi(t)}{2}\sqrt{\frac{1+L_1(t)}{2}},\\ Y(t)&=&\frac{L_3(t)\cos \frac{\psi(t)}{2}+L_2(t)\sin \frac{\psi(t)}{2}}{\sqrt{2(1+L_1(t))}},\\ Z(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 consider Example Four with

(8)   \begin{eqnarray*} I_1&=&1.0,\\ I_2&=&1.012686988782515,\\ I_3&=&3.306237422473038 \end{eqnarray*}

(9)   \begin{eqnarray*} m_1(0)&=&-0.544332842491675,\\ m_2(0)&=&0.729131780907662,\\ m_3(0)&=&-0.414811526666455. \end{eqnarray*}

m_3(0) is negative, the therefore we take option as in Eq. 3). We calculate t_0 that reproduces the initial angular momenta. The method is known from the previous posts:

(10)   \begin{equation*}am=\arctan(-m_1(0)/A_1,m_2(0)/A_2)=0.925158,\end{equation*}

(11)   \begin{equation*}t_0=\mathrm{EllipticF}(am,m)=3.17243.\end{equation*}

So we have t_0. At t=t_0 our angular momentum is the same as angular momentum in example 4 at t=0. We calculate q_0=q(t_0) from Eq'(7)

    \[q_0=W(t_0)+X(t_0)\,\mathbf{i}+Y(t_0)\,\mathbf{j}+Z(t_0)\,\mathbf{k},\]

(12)   \begin{eqnarray*} W(t_0)&=&0.435964,\\ X(t_0)&=&0.194343,\\ Y(t_0)&=&-0.0858987,\\ Z(t_0)&=-&0.874521 \end{eqnarray*}

thus

(13)   \begin{equation*}q_0=0.435964+0.194343\,\mathbf{i}-0.0858987\,\mathbf{j}-0.874521\,\mathbf{k}.\end{equation*}

We verify that we have indeed a unit quaternion:

    \[||q_0||^2=W(t_0)^2+X(t_0)^2+Y(t_0)^2+Z(t_0)^2=1.0.\]

Thus the inverse quaternion is the same as conjugated one:

    \[q_0^{-1}=0.435964-0.194343\,\mathbf{i}+0.0858987\,\mathbf{j}+0.874521\,\mathbf{k}.\]

In the Example 4 the initial quaternion, at t=0 is 1. Therefore in order for our solution to reproduce the initial data of the example we set

(14)   \begin{equation*}\tilde{q}(t)=q_0^{-1}\,q(t+t_0).\end{equation*}

The final time in Example 4 is t=10. Therefore we calculate

(15)   \begin{equation*}\tilde{q}(10)=q_0^{-1}q(10+t_0),\end{equation*}

where the multiplication is the quaternion multiplication.
Calculated with Mathematica the answer is

(16)   \begin{eqnarray*} W&=&=-0.36762,\\ X&=&-0.630629,\\ Y&=&-0.612723,\\ Z&=&0.302874. \end{eqnarray*}

The result from the Fortran code in the file out_example4.dat is

-0.3676198430772359
-0.6306293413288832
-0.6127232632258010
0.3028737154869889


It seems that therefore that we have obtained the ultimate answer. At least for quaternions. We need to get it also for rotation matrices. So the saga will continue.