Category: Euler’s equations

Approaching the ultimate answer for m<1

Posted on by arkajad

Far out in the uncharted backwaters of the unfashionable end of the Western Spiral arm of the Galaxy lies a small unregarded yellow sun.
Orbiting this at a distance of roughly ninety-eight million miles is an utterly insignificant little blue-green planet whose ape-descended life forms are so amazingly primitive that they still think digital watches are a pretty neat idea.
This planet has—or rather had—a problem, which was this: most of the people living on it were unhappy for pretty much of the time. Many solutions were suggested for this problem, but most of these were largely concerned with the movements of small green pieces of paper, which is odd because on the whole it wasn’t the small green pieces of paper that were unhappy.
And so the problem remained; lots of the people were mean, and most of them were miserable, even the ones with digital watches.

So sets Douglas Adams the scenery in his Hitchhiker’s Guide to the Galaxy. And he is right, neither green pieces of paper nor digital watches will make us happy. But digital computers – that’s another story. Digital computer, with some patience and determination made me happy today (or I made him happy, doesn’t mater which way), and I am going now to share my happiness with you. Because there is no true happiness which is not shared happiness.

And no, I am not going to calculate the answer to the Ultimate Question of Life, the Universe and Everything. Though we will be not very far away from obtaining such an answer. But please, don’t panic!

In the last post I have perversely asked this question:

So, my guess is that the formula (23) in the paper has wrong signs. What to do?
Should I write to the authors?

Now I admit, I was not completely sincere. I was teasing, I was playing with You, the Reader. Because the truth it that months ago I discovered wrong signs in the formula (20) of this paper. So, being young and naive, I wrote to prof. Zanna announcing my discovery. Here is the reply I got:

So, we should not consider the paper as being sacred, typos sometimes happen. But software – that is something different. Software is a serious matter. The software has been thoroughly tested, so if we test our software against something that has been thoroughly tested, that will be even better than through tests of our own code!

In The Hitchhiker’s Guide to the Galaxy the answer came as forty two.

“I checked it very thoroughly,” said the computer, “and that quite definitely is the answer.

We will also check very thoroughly. Our answer will not exactly be forty two. But we will come out with a quaternion, quaternio, quatro – it can explain 4. It remains to explain the last 2.

We will be testing example four from the code 903 – you need to read a couple of previous posts if you do not know what I am talking about. Alternatively, you can download the file 903.zip, unzip it, go to the folder 903\Fortran77\Drivers, and contemplate the following three files there: driver_example4.f, init_example4.dat, out_example4.dat. All three with number four. That is what we will be doing now.

Helsinky OH 52 - Emil ZÁTOPEK maraton
Emil Zátopek was a Czechoslovak long-distance runner best known for winning three gold medals at the 1952 Summer Olympics in Helsinki. He won gold in the 5,000 metres and 10,000 metres runs, but his final medal came when he decided at the last minute to compete in the first marathon of his life.

According to README.txt

driver_example3.f (example for FRB integration with many steps
of lenght h from 0 to Tfin, using
semi-exact methods)
To be linked with frb_step.f
driver_example4.f (Like driver_example3.f but uses quaternions)
To be linked with quat_step.f

Both examples have the body with the same moments of inertia, but initial momenta are different. So it is impossible to compare their outputs. We return to this issue later on, for now let us concentrate on example4.

The moments of inertia can be read from init_example4.dat. They are

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

So it is the same “future physics” body that we have already met in , when discussing example2.

Then comes initial angular momentum (which sometimes I call \mathbf{m} and sometimes \methbf{L}:

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

As in Spinning Gulliver among giants we have

(3)   \begin{equation*} 1/I_2= 0.987472.\end{equation*}

We calculate d

    \[d =\frac{m_1(0)^2}{I_1}+\frac{m_2(0)^2}{I_2}+\frac{m_3(0)^2}{I_3}=0.873315.\]

So d<1/I_2, therefore we have m<1 case.
Now I am presenting the complete m<1 algorithm based on the method introduced in Standing on the shoulders of giants – Reboot.

The case m<1.

(4)   \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

(5)   \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

(6)   \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

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

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

we set the phase variable \psi(t) as

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

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

My original intention was to verify the above with the 903 paper code using Mathematica, and, if there will be differences in the output, to fix the problems.
But now this Sunday is coming to the end, and I am getting tired. And before sleep I want to watch the second half of the movie “Master and Commander”. It fits this post:

Roger Ebert gave the film 4 stars out of 4, saying that “it achieves the epic without losing sight of the human”

The case of m>1

Posted on by arkajad

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!

Spinning Gulliver among giants

Posted on by arkajad

I am following giants. At least I am trying. For Gulliver it was not so funny to be in the land of giants:

Certainly Jonathan Swift was a giant of literature. Another giant, among those that I am trying to learn from, is G. K. Chesterton:

Gilbert Keith Chesterton was a giant in many ways. His six-feet-two inches and corpulent figure made him a physical giant but it was the extent and diversity of his writings that made him a giant in the field of literature and theology. He proved to be a great defender of the Catholic faith and crossed swords with the likes of George Bernard Shaw on theological matters yet they remained friends, and Shaw called him ‘a man of colossal genius’

Pioneer magazine Chesterton – The Genial Giant

Well well, we really should listen to giants:

The human race, to which so many of my readers belong, has been playing at children’s games from the beginning, and will probably do it till the end, which is a nuisance for the few people who grow up. And one of the games to which it is most attached is called “Keep to–morrow dark,” and which is also named (by the rustics in Shropshire, I have no doubt) “Cheat the Prophet.” The players listen very carefully and respectfully to all that the clever men have to say about what is to happen in the next generation. The players then wait until all the clever men are dead, and bury them nicely. They then go and do something else. That is all. For a race of simple tastes, however, it is great fun.
For human beings, being children, have the childish wilfulness and the childish secrecy. And they never have from the beginning of the world done what the wise men have seen to be inevitable. They stoned the false prophets, it is said; but they could have stoned true prophets with a greater and juster enjoyment. Individually, men may present a more or less rational appearance, eating, sleeping, and scheming. But humanity as a whole is changeful, mystical, fickle, delightful. Men are men, but Man is a woman.

Hmmm…. In my last note, Standing on the shoulders of giants – Reboot,
I started the discussion of one paper by two ladies, 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. But I did not connect in any way to the content of their paper. I am presenting here a long series of posts about rotations of free rigid bodies, but it is in the spirit of playing games like children do. While the ladies with their paper are serious. I remember what the Bible says:

Yes, it is time to put away childish things or, better, to learn how to talk like an adult to adults. So, let’s do it.

There is the paper, and there is the code. Both can be downloaded and examined. More important is the code. Paper may contain typos. Papers often do. Paper you read. Code you run. There is the code, and there are examples of its application. This is a serious matter. We should be able to check if are able to reproduce the results. Let’s take example 2. Here are initial data that we will try to understand:
1.0000000000D0
10.0000000000D0
-0.609860759302936D0
0.761660947972381D0
0.218957654801698D0
1.000000000000000D0
0.000000000000000D0
0.000000000000000D0
0.000000000000000D0
1.000000000000000D0
0.000000000000000D0
0.000000000000000D0
0.000000000000000D0
1.000000000000000D0
1.000000000000000D0
1.012686988782515D0
3.306237422473038D0

First we have time step. Then we have time. Time is t=10.
Then we have initial angular momentum (in the body frame)

(1)   \begin{align*}m_1(0)&=-0.609860759302936,\\ m_2(0)&= 0.761660947972381,\\ m_3(0)&=0.218957654801698. \end{align*}

I checked and we get indeed norm 1.
Then we have the initial attitude matrix Q(0)=I – the identity matrix.
The last three lines are the moments of inertia:

(2)   \begin{align*} I_1&=1.000000000000000,\\ I_2&= 1.012686988782515,\\ I_3&= 3.306237422473038. \end{align*}

The fact that I_1+I_2<I_3 tells us that mathematicians (our two ladies) have more imagination than physicists. At present we do not have materials that would provide such data (for positive mass densities we should always have I_1+I_2\geq I_3,) but in the near future physics and technology may easily move to the place where mathematicians are already today. It is reassuring that already today we can predict and model the behavior of such exotic materials.

Let us calculate the parameter d (the ratio of doubled kinetic energy to the square of angular momentum)

    \[d =\frac{m_1(0)^2}{I_1}+\frac{m_2(0)^2}{I_2}+\frac{m_3(0)^2}{I_3}=0.95929.\]

Since 1/I_2= 0.987472, we have the case when d<1/I_2. I will use the Mathematica code from our children games:

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

Using these formulas we get:

(3)   \begin{align*} A_1&=0.97038,\\ A_2&=0.979214,\\ A_3&=0.241582,\\ B&=0.166993,\\ m&=0.29508. \end{align*}

We see that m<1, as we expected.

We will use the Mathematica code from our children games:

Looking at the initial angular momentum data we see that m_3(0) is positive, therefore we can stay with the formulas above, we do not need to change signs of two of the angular momenta inside the code.

Now, what do we do with initial angular momenta? In our code we need to find the time t_0 at which we obtain these values. From that time on we will have to run the code for t=10. So our problem is to solve the following equations:

(4)   \begin{align*} m_1(0)&=A1 \cn(B t_0, m),\\ m_2(0)&=A2 \sn(B t_0, m). \end{align*}

We do not have to worry about m_3(0), because it is already determined by m_1(0) and m_2(0).
We know that

    \[\cn(B t_0, m)=\cos \am(B t_0,m),\]

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

Thus, using Mathematica ArcTan function:

    \[\am(B t_0,m)=ArcTan[m_1(0)/A_1,m_2(0)/A_2].\]

I use Mathematica and get

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

But \am is the inverse function of EllipticF. So I calculate t_0

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

So we know t_0. Knowing t_0 I use the code above to calculate Q(t_0). We get

(5)   \begin{equation*}Q(t_0)=\begin{bmatrix}\ 0.68312 & 0.365156 & 0.632462 \\ 0.401768 & 0.535288 & -0.743 \\ -0.609861 & 0.761661 & 0.218958 \end{bmatrix}.\end{equation*}

In the example from the code we have Q(t_0)=I. therefore to get what they get I need to calculate \tilde{Q}(t)=Q(t_0)^{-1}Q(t) at t=t_0+10.
We use the Mathematica code above to get:

(6)   \begin{equation*}\tilde{Q}(t_0+10)=\begin{bmatrix} 0.751185 & -0.123316 & -0.64847 \\ -0.165911 & -0.98613 & -0.0046633 \\ -0.638901 & 0.111091 & -0.761226\end{bmatrix}\end{equation*}

To make comparison easier with the results obtained by Celledoni and Zanna, let us transpose the matrix, because they write the output columnwise.

(7)   \begin{equation*}\begin{bmatrix} 0.751185 & -0.165911 & -0.638901 \\ -0.123316 & -0.98613 & 0.111091 \\ -0.64847 & -0.0046633 & -0.761226 \end{bmatrix} \end{equation*}

And here is the result from the code of Celledoni and Zanna at the end of out_example3.dat file:

0.7511853766362874
-0.1659110385305551
-0.6389006630310473
-0.1233163682782040
-0.9861296978083955
0.1110913696692821
-0.6484702022780026
-4.6633029105353817E-03
-0.7612257551892809

It seems that we can now speak the language of adults as well. This encouraging and we can continue our adventure.
The Mathematica notebook with these calculations can be downloaded from here.
Being at peace with matrices, we can move to more dangerous quaternions.