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.

13 thoughts on “The final answer for the Universe in which m<1”

  1. Something fishy

    Yes, something fishy is going on.
    The new algorithm, presented here, though it seems to agree with 903 code, it is not in agreement with my old algoritm for quaternions. For large t they produce certainly different numbers.
    I have no idea why.
    To debug the universe I will have ro start with rotation matrices rather than quaternions, which are somewhat exotic species.
    So the next post will be about the new algorithm with rotation matrices, and then I will start the debugging hearing in the Congress.

  2. Maple’s EllipticF is it defined differently to Mathematica?
    {\it mambiguous} \left( {\it mambiguous} \left( `F \left(z ,k \right) \mbox {{\tt ```}}:=\mbox {{\tt `'`}} \int _{0}^{z }\frac{1}{\sqrt{- \alpha \mbox {{\tt `1`}} ^{2}+1}~\sqrt{-k ^{2}~\alpha \mbox {{\tt `1`} } ^{2}+1}}d \alpha \mbox {{\tt `1`}} \esapos , \right) \right)
    from Spinning Gulliver among giants you calculate 14.9607 is got 5.975-12.471i.

    Today a got a real number but now yours.
    Edit:-
    Latex got messed up a bit but it shows the function. Basicially F(z,k) goes to infinity just after 5 pie /16 and doesn’t become real again.

    1. Did a bit of research found this inMaple help

      “It is worth noting the difference between the Legendre normal form of the Incomplete Elliptic integral of the first kind (see A&S 17.2.7), in Maple represented by EllipticF(z,k) but for the splitting of the square root in the denominator of the integrand (see definition lines above), and the normal trigonometric form of this elliptic integral (see A&S 17.2.6), in Maple represented by the InverseJacobiAM function
      InverseJacobiAM(phi,k);
      InverseJacobiAM(phi, k)”
      Gives both answers for t0, 14.9607 and todays 3.17243

  3. Problem solved. New problem:-
    psi(t) is probably causing an issue.
    Can you give the value for psi(t0) along with the arctan and EllipitcPi functions in it please?
    Must sleep now.

    1. The following gives the same answer as Mathematica

      I1 := 1.0; I2 := 1.012686988782515; I3 := 3.306237422473038; d := .8733147046233154; m := (I3*d-1)*(I2-I1)/((-I1*d+1)*(I3-I2)); k := sqrt(m);
      A2 := sqrt(I2*(I3*d-1)/(I3-I2)); A3 := sqrt(I3*(-I1*d+1)/(I3-I1)); alpha := (I3-I1)/sqrt(I1*(-I1*d+1)*(I3-I2)*I3/I2); nu := (-I1*I3*d+I1)/(-I1*I3*d+I3); B := sqrt((-I1*d+1)*(I3-I2)/(I1*I2*I3)); t0 := 3.17243;
      psi(t):=t/(I1)-alpha*epi(B*t,nu,k)-arctan((A2)/(A3)JacobiSD(B t,k));

      psi(t0);
      0.838675316

  5. Can you explain this formula please?

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

    Don’t understan it? But then I would need you to explain why we need a spacial t_0. What is wrong with 0 or why would zero be wrong?

    1. A general solution of equations of motion involves integration constants that would allow accomodating arbitrary initial conditions. We do not have such constants, so we accomodate initial data in a different way. First we need to have correct angular momenta. For this we are shifting our solution by t0. Then we need to take into account initial q. For this we shift the shifted solution by multiplication from the left by q_0^{-1}.
      If we would set

          \[\tilde{q}(t)=q(0)^{-1}q(t),\]

      then we would have \tilde{q}(0)=1 satisfied, but that would be a solution with wrong initial angular momenta (that is, it would correspond to a different initial values of \dot{q}(0)).

      If the above is not clear, please let me know, and I will will elaborate more.

      1. Edit- copy paste of letex script hasn’t worked.Removing $ signs worked but doesnt directly display equation.
        My basic question concerned http://quicklatex.com/cache3/ce/ql_39e56426e65caff39b00d673d49df4ce_l3.png
        That deals with quaternion “logic”. So far I only have a rudimentry level of knowlege of them. It is a tough situation even to figure out the revelant question to ask. I don’t want to waste you time.
        So 2 questions,
        1) If it is not to time consuming to explain,I would appreciate you pointing me to a reference that explains why you need to multiply by q(0)^-1 intoq(t+t0).
        2) What does the ~ over q(t) actually indicate/mean/represent.

Leave a Reply