### New wine in new bottles

In my recent post I was presenting new formulas for the attitude matrix, different from those of a month ago. And I did not explain in detail how I got them, and why they look like being in disagreement with the previous ones. This caused a legitimate concern from Bjab, when he noticed the apparent disagreement. And indeed putting new wine into old bottles would lead to a trouble.

Mar 2:21 No man also seweth a piece of new cloth on an old garment: else the new piece that filled it up taketh away from the old, and the rent is made worse.
Mar 2:22 And no man putteth new wine into old bottles: else the new wine doth burst the bottles, and the wine is spilled, and the bottles will be marred: but new wine must be put into new bottles.

What we did is this: we poured new wine into new bottles. And here is how it is done.

We consider free rigid body with normalized (unit length) angular momentum vector, and the doubled kinetic energy equal to Euler’s equations can be then solved explicitly and we choose the following solution:

(1)

(2)

where
We will need the angular velocity vector with components

(3)

We want to solve now the differential equation for the attitude matrix

(4)

where

(5)

The strategy is this: We seek the solution of the form where is a certain matrix cleverly constructed out of and is a simple rotation matrix by some angle
There are two methods: old one, and new one.

### The old method.

In the old method, as used for instance in Taming the T-handle continued (we use now instead of there), we set:

(6)

(7)

(8)

Then we solve for :

(9)

(10)

(11)

(12)

where

(13)

### Comparing old and new wines

In both cases we get a solution of the attitude differential equation (4). But these are different solutions. If is a solution, and if is a fixed rotation matrix of determinant one, then is another solution. This is the same body rotating, but observed from two different laboratory frames.
Let and be the two solutions, one obtained with the old, the other one with the new method.
With some little algebra it can be verified that where

(14)

.
Therefore for all

(15)

I verified it numerically. The above comparison was for rotation matrices. For quaternions it is similar, though the quaternion case deserves a comment. I will give it in another note.

Update: Mathematica notebook verifying Eq. (15) numerically is available here.

### Attitude matrix for m<1

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….

So, here is the code:

The case
The solution of the Euler’s equations has two trajectories. For one with we take

(1)

while for the other one, with we take

(2)

Let

(3)

The solution of the Euler’s equations is given by

(4)

With constants and defined as

(5)

(6)

we set the phase variable as

(7)

where the Jacobi function is defined as
Let

(8)

(9)

Then the attitude matrix is given by

(10)

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 rather than with the signs in formulas for 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 that depends on the trajectory we are taking. The constant 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 I use Mathematica to calculate the transposed final

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 is fixed, and equal ! 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 I hope I am right here….

### The final answer for the Universe in which m<1

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

(1)

The solution of the Euler’s equations has two trajectories: one with is given by

(2)

while the other one, with is given by

(3)

With constants and defined as

(4)

(5)

we set the phase variable as

(6)

where the Jacobi function is defined as
Then the quaternionic attitude solution is given by with

(7)

We consider Example Four with

(8)

(9)

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

(10)

(11)

So we have At our angular momentum is the same as angular momentum in example 4 at We calculate from Eq'(7)

(12)

thus

(13)

We verify that we have indeed a unit quaternion:

Thus the inverse quaternion is the same as conjugated one:

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

(14)

The final time in Example 4 is Therefore we calculate

(15)

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

(16)

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.