The case of m>1

When I was in the middle of writing this post, question came in a comment from Bjab

New Note after posting everything here: As has just been pointed to me by Ronan, I have already covered this issue in Infeld, Einstein and blogging . While, there is no harm in covering the same issue twice, the fact that I have forgotten tells something about the workings of my memory. But collective memory (of all readers) seems to be the answer.

Ark,
where can I find, in your posts, formulas for A1, A2, A3, B, m for m>1 ?

(For m<1 and m=1 they are in posts:"Taming the T-handle" and "Taming the T-handle – continued"

I checked, and indeed I could not find the place where I have treated this problem with clarity. It is my ambition to write extremely clear, better in this respect than anyone else in the whole universe, but here we have a hole. Since this post was anyway supposed to be exactly about the case of m>1, I will start with describing the case in details (and, as we know, the devil is hiding just there).

To begin with the beginning: when considering a free rigid body rotating about its center of mass, the most important parameter characterizing the rotation is the parameter d defined as the ratio of the doubled kinetic energy (that is constant during the motion) to the square L^2 of the angular momentum vector (that is also constant during the motion). In order to simplify the notation let us assume that L is normalized so that L=1. Then

(1)   \begin{equation*}d=\frac{L_1^2}{I_1}+\frac{L_2^2}{I_2}+\frac{L_3^2}{I_3},\end{equation*}

where I_1<I_2<I_3 are the principal moments of inertia, and L_1,L_2,L_3 are the components of the angular momentum vector with respect to the frame attached to the body. It follows almost directly from the definitions that the value of d is in the interval 1/I_3\leq d\leq 1/I_1. The case of d=1/I_2 is very special, there is only one flip of the body – it is “the archetype of Dzhanibekov effect”.
The first step in finding explicit solution to the equations of motion is to solve Euler’s equations:

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

A solution is given by the following explicit formula adapted from Solving Euler’s equations:

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

In Solving Euler’s equations the solution is given for \Omega_i. To get the solution for L_i we simply multiply the constants A_i from there by I_i to obtain

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

In the formulas (3) Jacobi elliptic functions \cn,\sn,\dn are somewhat formally extended to cover the case of m>1 that corresponds to d>1/I_2. It is safer to transform these formulas using the conversion formulas (see More is different and links from there)

(5)   \begin{equation*}\mathrm{cn}(u,m)=\mathrm{dn}(ku,1/m),\quad m>1\end{equation*}

(6)   \begin{equation*}\mathrm{sn}(u,m)=\frac{1}{k}\mathrm{sn}(ku,1/m)\quad m>1\end{equation*}

(7)   \begin{equation*}\mathrm{dn}(u,m)=\mathrm{cn}(ku,\frac{1}{m}),\quad m>1.\end{equation*}

Let us call 1/m=\mu. Then, with m=k^2, L_1(t), L_2(t),L_3(t) become

(8)   \begin{eqnarray*}  L_1(t)&=&A_1\,\cn(Bt,m)=A_1\,\dn(kBt,\mu),\\ L_2(t)&=&A_2\,\sn(Bt,m)=\frac{1}{k}A_2\,\sn(kBt,\mu),\\ L_3(t)&=&A_3\,dn(Bt,m)=A_3\,\cn(kBt,\mu). \end{eqnarray*}

We see that perhaps it is useful to introduce new constant A_2 that is equal to the old A_2 divided by k. We have

    \[1/k=\sqrt{1/m}=\sqrt{\frac{(1-dI_1)(I_3-I_2)}{(dI_3-1)(I_2-I_1)}}.\]

Therefore, after some little calculation, we get

    \[ A_2/k=\sqrt{\frac{I_2(dI_3-1)}{I_3-I_2}}\sqrt{\frac{(1-dI_1)(I_3-I_2)}{(dI_3-1)(I_2-I_1)}}=\sqrt{\frac{I_2(1-dI_1)}{I_2-I_1}}.\]

This will be our new A_2. It is also useful to introduce new B that equals to k times old B. We have

    \[kB=\sqrt{\frac{(dI_3-1)(I_2-I_1)}{(1-dI_1)(I_3-I_2)}}\sqrt{\frac{(1-dI_1)(I_3-I_2)}{I_1I_2I_3}}=\sqrt{\frac{(dI_3-1)(I_2-I_1)}{I_1I_2I_3}}\]

This will be our new B. Finally, we can now rename old 1/m into m. Altogether we obtain

If d>1/I_2 we define:

(9)   \begin{eqnarray*} m&=&\frac{(1-dI_1)(I_3-I_2)}{(dI_3-1)(I_2-I_1)},\\ 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}}. \end{eqnarray*}

A solution of the Euler’s euqations for this case can be written as

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

I wrote “a solution”. In fact it is almost the most general solution. Every solution can be obtained by shifting the time t\mapsto t-t_0, and, if necessary, putting minus signs in front of two of the three formulas for L_i(t).
That was an extra for this post because of the question from Bjab. Now comes what I was going to post today originally.

This is a continuation from the previous post Spinning Gulliver among giants. Another example from the 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

Elena Celledoni has been employed at the Department of Mathematical Sciences since 2004. She is a professor of mathematics since 2009. She is a member and currently the leader of the group of differential equations and numerical analysis.

I will make it short and sweet, because the reasoning will be like in Spinning Gulliver among giants. At the same time I will stress the differences. This should also help to understand the algorithm I am using for the case d>1/I_2. I am using Mathematica, and in the forthcoming posts I will mention how to adapt it to Maple (which at the same time will answer the comment from Ronan).

I am taking now “example0” from Fortran code examples that accompany the paper by Celledoni and Zanna. These examples can be downloaded in a zip file using the link to the code here. (though I got the file directly from one of the Authors – A.Z.).

The description of example0 in a Readme file: “driver_example0.f (example of FRB integration in a single step from 0 to Tfin). To be linked with frb_step.f”.

The initial data for this example are (I am decoding here their meaning)

Time to run:
1.0000000000000000e-01
Moments of inertia:
1.0000000000000000e+00
1.6487857827119290e+00
1.9720127096641928e+00
Initial angular momentum
-0.709894965287627e+00
-0.685144717153487e+00
0.163174308075589e+00
Initial attitude matrix (columnwise|)
-0.139359936806125e+00
-0.414251479478192e+00
-0.899430108326112e+00
-0.726901415681411e+00
-0.574008408466476e+00
0.376999574124619e+00
-0.672453076350874e+00
0.706335655874678e+00
-0.221126211350745e+00

At first I was playing with these data. But then, reading Readme.txt, I noticed that: “driver_example2.f (Like driver_example0.f, but uses quaternions in place for rotations). To be linked with quat_step.f“. So I checked example2, and I noticed that there the initial attitude matrix is given by the unit quaternion 1. That is the initial attitude matrix is the identity matrix as it was in the example we have already discussed in Spinning Gulliver. Therefore I decided to change the example and to use the identity matrix here as well. So I will be discussing the slightly modified version of the original example.

To Mathematica I enter these data:

Tfin = 0.1;
I1 = 1;
I2 = 1.6487857827119290;
I3 = 1.9720127096641928;
L10 = -0.709894965287627;
L20 = -0.685144717153487;
L30 = 0.163174308075589;

I checked that the initial angular momentum vector has indeed length 1, as it is supposed to have. We now calculate d

    \[d =\frac{L10^2}{I_1}+\frac{L20^2}{I_2}+\frac{L30^2}{I_3}=0.802161.\]

Since 1/I_2= 0.606507, we have the case when d>1/I_2.

This is more difficult case. The modulus k (and m=k^2) would come greater than 1. While Mathematica has elliptic functions with m>1, other software may have them only for m<1. We were discussing before how to convert from m to 1/m. Here I will be using the results of our previous discussions. So, m below will be already inverted, therefore smaller than 1. Here are the formulas for our case:

A1 = Sqrt[I1*(d*I3 – 1)/(I3 – I1)]
A2 = Sqrt[I2*(1 – d*I1)/(I2 – I1)]
A3 = Sqrt[I3*(1 – d*I1)/(I3 – I1)]
B = Sqrt[((d*I3 – 1)*(I2 – I1))/(I1*I2*I3)]
m = ((1 – d*I1)*(I3 – I2))/((d*I3 – 1)*(I2 – I1))

We get

A1 = 0.773709
A2 = 0.709067
A3 = 0.633541
B = 0.340743
m = 0.169391

In our simulations for this case in many previous posts we were using:

L1 = A1*JacobiDN[B*t,m]
L2 = A2*JacobiSN[B*t,m]
L3 = A3*JacobiCN[B*t,m]

But now we see that L10 is negative. This is a different Poinsot orbit. Therefore we have to change the sign of L1. We know that when we want to have Euler’s equations still satisfied, we have to change also one other sign. Let us choose L2 ( I could instead change L3, then t_0 would be different, but all the conclusions would be the same). Thus what we will use is

L1[t_] = -Re[A1*JacobiDN[B*t, m]];
L2[t_] = -Re[A2*JacobiSN[B*t, m]];
L3[t_] = Re[A3*JacobiCN[B*t, m]];

We have put minus sign in the formula for L1[t] and L2[t].

The formula for \psi is

    \[\psi (\text{t$\_$})=\frac{(\text{I3}-\text{I1}) \Pi \left(-\frac{\text{I3} (1-d \text{I1})}{\text{I1} (d \text{I3}-1)};\left.\text{am}\left(\left.t \sqrt{\frac{(\text{I2}-\text{I1}) (d \text{I3}-1)}{\text{I1} \text{I2} \text{I3}}}\right|m\right)\right|m\right)}{\sqrt{\frac{\text{I1} \text{I3} (d \text{I3}-1) (\text{I2}-\text{I1})}{\text{I2}}}}+\frac{t}{\text{I3}};\]

Note: I have not yet discussed the formula for \psi for d>1/I_2. It needs to be discussed. But we will not be using it in this post anyway.

The rest goes like in the previous post. As before, we compute the amplitude

Thus, using Mathematica ArcTan function:

    \[\am(B t_0,m)=ArcTan[L30/A3,-L20)/A2].\]

Notice that here I had to put minus sign in front of L20. That is because we have changed the sign in the formula for L2.
I use Mathematica and get

    \[\am(B t_0,m)=1.3103.\]

Then

    \[t_0=EllipticF[2.25039,m]/B=3.98931.\]

The output file from the Fortan program gives both angular momentum and the attitude matrix at the final time. Let us first check angular momentum. Notice that even if I have changed the initial attitude matrix, this change does not affect the angular momentum. Therefore I can read the final angular momentum from the original output file. Here it is:

-0.708844791922432
-0.690514805418774
0.143973485273913

We go to Wolfram Alpha, enter there our data for L1(t_0+0.1):

0.773709 JacobiDN[0.340743*(3.98931+0.1,0.169391)

and obtain the result: 0.7088446163317107…
Comparing with 0.70884479192243 the result is not too bad.

The Reader can check L2, and L3 (I have already checked with Mathematica).

As this post became longer than I indented, and as there are things that I have to ponder about before next post, let us call it a day!

24 thoughts on “The case of m>1

  1. I compared the aboue to what I coded the other night. A1,A2/k ,A3 B.k m=1/k^2 agree.
    L1,L2,L3 all comeout positive but numerically agree whic the figures you give.

    EllifticF is defined as ${\it EllipticF} \left( z,k \right) =\int_{0}^{z}\!{\frac {1}{\sqrt {-{
    {\it \_}\alpha \mbox {{\tt `1`}}}^{2}+1}\sqrt {-{k}^{2}{{\it \_}
    \alpha \mbox {{\tt `1`}}}^{2}+1}}}\,{\rm d}{\it \_}\alpha \mbox {{\tt $ which must be different to mathematica. No real values after about 1.6 on the x axis.
    Looking at your formula in Infield… post I deduced epic uses
    ${\it epi} \left( Bt,-{\frac {{{\it A3}}^{2}}{{{\it A1}}^{2}}},{k}^{-1}
    \right) $
    asssuming I have epi correct my psi is
    %{\frac {t}{{\it I3}}}+{({\it I3}-{\it I1}){\it epi} \left( Bt,-{\frac
    {{{\it A3}}^{2}}{{{\it A1}}^{2}}},{k}^{-1} \right) \left( \sqrt{{
    \frac {{\it I1}\,{\it I3}\, \left( {\it I2}-{\it I1} \right) \left( d
    {\it I3}-1 \right) }{{\it I2}}}} \right) ^{-1}}$

    There is a typo in am(B*t_0,m)=Acrtan( error in here)

    I don’t get the having to have a special value for t_o. I ran it from t=-50..50 and t=0..100. Doesn’t mean whatI ran is correct though.
    Anyhow greal work!
    I actually ran into the problem of issues with results printed in papers too. Was trying to model the forward kinematics of a Stewart Gough Platform (Hexapod). Some of the papers would have completly incorrect tables of results printed in them. So I know aboult the trials and stresses that can cause.

  2. I compared the above to what I coded the other night. A1,A2/k ,A3 B.k m=1/k^2 agree.
    L1,L2,L3 all comeout positive but numerically agree whic the figures you give.

    EllifticF is defined as {\it EllipticF} \left( z,k \right) =\int_{0}^{z}\!{\frac {1}{\sqrt {-{ {\it \_}\alpha \mbox {{\tt `1`}}}^{2}+1}\sqrt {-{k}^{2}{{\it \_} \alpha \mbox {{\tt `1`}}}^{2}+1}}}\,{\rm d}{\it \_}\alpha \mbox {{\tt  `1`}}
    No real values after about 1.6 on the x axis.
    Looking at your formula in Infield… post I deduced “epi” uses
    {\it epi} \left( Bt,-{\frac {{{\it A3}}^{2}}{{{\it A1}}^{2}}},{k}^{-1}  \right)
    asssuming I have epi correct my psi is
    {\frac {t}{{\it I3}}}+{({\it I3}-{\it I1}){\it epi} \left( Bt,-{\frac  {{{\it A3}}^{2}}{{{\it A1}}^{2}}},{k}^{-1} \right)  \left(  \sqrt{{ \frac {{\it I1}\,{\it I3}\, \left( {\it I2}-{\it I1} \right)  \left( d {\it I3}-1 \right) }{{\it I2}}}} \right) ^{-1}}
    There is a typo in am(B*t_0,m)=Acrtan( error in here)

    I don’t get the having to have a special value for t_o. I ran it from t=-50..50 and t=0..100. Doesn’t mean whatI ran is correct though.
    Anyhow greal work!
    I actually ran into the problem of issues with results printed in papers too. Was trying to model the forward kinematics of a Stewart Gough Platform (Hexapod). Some of the papers would have completly incorrect tables of results printed in them. So I know aboult the trials and stresses that can cause.

      1. Sorry if the communication was abrupt.

        you have am(B*t_0,m)=ArcTan[L30/A3,-L20)/A2]

        there is “,” and “)” should it be
        am(B*t_0,m)=ArcTan[L30/A3-L20/A2]

        Do you type in the LaTeX code yourself? I get Maple to produce mine, otherwise I would post unintelligable garbage.

        1. It is OK. From Maple help:

          For real arguments x, y, the two-argument function arctan(y, x), computes the principal value of the argument of the complex number x+I*y, so -Pi < arctan(y, x) <= Pi. In Mathematica the first argument is x the second y.

  3. ” the fact that I have forgotten tells something about the workings of my memory.”

    Well, I asked the question:
    “Where can I find, in your posts, formulas for A1, A2, A3, B, m for m>1 ?”
    so I also was not quite sure if you managed the case of m>1. I was looking for it but I couldn’t find it, so I asked. I was not expecting you would write the whole new post on it but that you would show me the previous post on that subject.
    Anyway this new post helped me to clarify why A2 and B are so different for m>1 and m<1 so thanks a lot.

  4. Another thing which you wrote in this post was not quite enough recognised by me in former posts, that is:
    “solution can be obtained […] putting minus signs in front of two of the three formulas for L_i(t).”

    I have a feeling that this can be the answer for my former question:
    “From symmetry there should be fifth term:
    \pi/3 or \frac{4}{3}\pi ?”

    So may be you would make another update to your post:”Never twice the same” taking this second solution into account.

  5. You wrote in that post:
    \Omega_1(t)=(A_1/I_1)/\cosh(Bt),\Omega_3(t)=(A_3/I_3)/\cosh(Bt).

    \[A_1/I_1=1/2,\, A_3/I_3=\frac{1}{2\sqrt{3}}.
    Which gaved angle of \pi/3.

    But there is another solution:
    \[A_1/I_1=-1/2,\, A_3/I_3=-\frac{1} {2\sqrt{3}}.
    which gives angle \frac{4}{3}\pi.

    So the attacking side of a T-handle can be red or with the other solution can be the other one (white?).
    (Notice that these colours are as on flag of Poland.)

    1. It is simply changing the direction of the vector. In fact we want to know the angle between a vector and the x axis (line). It should always be less than pi. If it is not, we subtract pi.

      1. “In fact we want to know the angle between a vector and the x axis (line).”

        You want what you want.

        But I want to know the angle between a velocity vector and the x axis (not simply a line but directed line). This gives me the opportunity to distinguish the sides of T-handle. I don’t want to be colorblinded. That gives me also the knowledge by looking at simulations of movement of colored T-handle whether m>1 or m<1.

          1. As ever. Left is the other than right. Positive charge is the other than negative. Imaginary unit i is the other than minus i.

            What is important: after painting the plane of the T-handle it is painted – the sides are easly distinguishable. Even if not painted the sides of T – handle are distinguishable (but harder).

            So if there was not the second solution (4pi/3) there wouldn’t be symetry.

        1. But I prefer to say it this way: you can always paint in such a way that the angle is pi/3. I call it: painting method adjusted to the actual situation.

          In general I advocate flexibility.

          1. So you will have to repaint the T-handle every each flip when m>1. A waste of time and paint.

          2. The angle is constant in time when m is exactly 1 but then there is no next flip so “you will have to repaint …”.

            I mentioned the case m>1.
            What if m>1 (even slightly)? I think that in that case the angle is near pi/3 and then, after a flip, near 4pi/3.

          3. That is a new interesting problem. Will look at it.
            OK. I am looking. The angle with x-axis is always between -pi/2 and pi/2 for a given painting.

            But probably we are not really interested in this angle. I believe we should be interested in the angle with z-axis?

          4. “The angle with x-axis is always between -pi/2 and pi/2 for a given painting.”

            It is in contrast to what I wrote at 14:08

          5. Omega1 is (A1/I1) dn, so it has a constant sign.
            But it does not tell us anything about the angle with (x,y) plane which oscillates like Omega3 that is proportional to cn.

          6. If angle of velocity vector with x-axis is pi/3 or 4pi/3 so the angle with z-axis is accordingly pi/6 or 7pi/6.

Leave a Reply