Category: Elliptic functions

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.

Standing on the shoulders of giants – Reboot

Posted on by arkajad

Standing on the shoulders of giants may be good, even very-very good. But at the same time it can be very dangerous. First of all it is a dangerous balancing act. It is easy to fall.

But what if the giant gets a hiccup? Or stumble? You may die. Or someone can die because of aftershocks.

In Standing on the shoulders of giants I quoted two papers that I used in working out my computer simulations of Dzhanibekov’s effect. These were

Ramses van Zon, Jeremy Schofield, “Numerical implementation of the exact dynamics of free rigid bodies“, J. Comput. Phys. 225, 145-164 (2007)
and
Celledoni, Elena; Zanna, Antonella.E Celledoni, F Fassò, N Säfström, A Zanna, The exact computation of the free rigid body motion and its use in splitting methods, SIAM Journal on Scientific Computing 30 (4), 2084-2112 (2008)

The algorithms were working smoothly. But recently I decided to use quaternions instead of rotation matrices. And in recent post Introducing geodesics, I have presented a little animation:

It looks impressive, but if we pay attention to details, we see something strange:

We have a nice curly trajectory, but there are also strange straight line spikes. On animation they are, perhaps, harmless. But if something like this happens with the software controlling the flight of airplanes, space rockets or satellites – people may die. There is a BUG in the algorithm. Yes, I did stay on shoulders of giants, but perhaps they were not the right giants for my aims? And indeed this happens. The method I used is not appropriate for working with quaternions.

That is why I have to do the reboot, and look for another giant. Here is my giant:
Antonella Zanna, professor of Mathematics, University of Bergen, specialty Numerical Analysis/Geometric Integration.

I am going to use the following 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.

The paper is only seven years old. And it has all what we need. As in Standing on the shoulders of giants we will split the attitude matrix Q into a product: Q=Q_1Q_0, but this time we will do it somewhat smarter. Celledoni and Zanna denote the angular momentum vector in the body frame using letter \mathbf{m}=(m_1,m_2,m_3), and they assume that \mathbf{m} is normalized: m=\sqrt{\mathbf{m}^2}=1. They assume, as we do, I_1<I_2<I_3, and use T to denote the kinetic energy. Let us first look at the following part of their paper (which is available here):

Though it is not very important, I must say that it is very nice of the two authors that they acknowledge Kosenko. I was trying to study Kosenko’s paper, but I gave up not being able to understand it. Then wen notice that the case considered is that of d>1/I_2 using the notation in my posts. That means “high energy regime”.
The angular momentum \mathbf{m} is constant (fixed) in the laboratory frame owing to the conservation of angular momentum. Suppose we orient the laboratory frame so that its z axis is oriented along the angular momentum vector. Then, in the laboratory frame angular momentum has coordinates (0,0,1). Suppose we construct a matrix Q_0 that transforms \mathbf{m}=(m_1,m_2,m_3) into (0,0,1). The attitude matrix Q transforms vector components in the body frame into their components in the laboratory frame. In particular it transforms \mathbf{m}=(m_1,m_2,m_3) into (0,0,1). Therefore,if we write Q=Q_1Q_0, the matrix Q_0 must transform (0,0,1) into itself. Therefore Q_0 must be a rotation about the third axis, that is it is of the form:

(1)   \begin{equation*}Q_1= \begin{bmatrix}\cos\psi&-\sin\psi&0\\\sin\psi&\cos\psi&0\\0&0&1\end{bmatrix}.\end{equation*}

All of that we have discussed already in Standing on the shoulders of giants. But now we are going to choose Q_0 differently.

So suppose we have vector \mathbf{m} of unit length. What would be the most natural way of constructing a rotation matrix that rotates \mathbf{m} into (0,0,1)? Simple. Project \mathbf{m} vertically onto (x,y) plane. Draw a line \mathbf{n} on (x,y) plane that is perpendicular to the projection. Rotate about this line by the angle \theta between \mathbf{m} and \mathbf{k}=(0,0,1).
Let us now do the calculations. Vector \mathbf{m} has components (m_1,m_2,m_3). Its projection on (x,y) plane has components (m_1,m_2,0). Orthogonal vector in (x,y) plane has components (m_2,-m_1,0) (to check orthogonality calculate scalar product), So the normalized vector \mathbf{n} has components

    \[\mathbf{n}=(\frac{m_2}{m_p}, -\frac{m_1}{m_p},0),\]

where

    \[m_p=\sqrt{m_1^2+m_2^2}.\]

The cosinus of the angle \theta is m_3. Its sinus is \sqrt{1-m_3^2}=m_p.
In Spin – we know that we do not know we have derived a general formula for a rotation about an angle \theta about axis defined by a unit vector \mathbf{n}:

    \[R=I+\sin\theta\, W(\vec{n})+(1-\cos\theta)\,W(\vec{n})^2,\]

where for any vector \mathbf{v}

(2)   \begin{equation*}W=W({\vec{v})=\begin{pmatrix}0& -v_3&v_2\\v_3&0&-v_1\\-v_2&v_1&0.\end{pmatrix}\end{equation*}

Applying the formula to our case, simple algebra leads to:

(3)   \begin{equation*}Q_0=\begin{bmatrix}(m_2^2 + m_1^2 m_3)/m_p^2& (m_1 m_2 (-1 + m_3))/m_p^2& -m_1\\ (m_1 m_2 (-1 + m_3))/m_p^2& ( m_1^2 + m_2^2 m_3)/m_p^2& -m_2\\m_1& m_2& m_3\end{bmatrix}.\end{equation*}

We can easily check that indeed Q_0 is an orthogonal matrix with determinant 1, Q0 acting on (m_1,m_2,m_3) vector gives (0,0,1). One can also check that the vector (-m_2,m_1,0) is invariant under Q_0 – as it should be as it is on the rotation axis. REDUCE program that checks it all is here.

At this point we have to make a break. It is not yet clear what is the relation of my Q_0 to the paper quoted, what is the difference between this and the old version, and why is this version better?

We will continue climbing on the shoulders of giants in the following notes. But, please, remember, Standing on shoulders of giants (like “Nobel Prize Winner” Obama), is risky:

Photo of Kellyanne Conway kneeling on Oval Office couch sparks…

Elliptic Pi in Mathematica and Maple

Posted on by arkajad

We use EllipticPi when we write exact solutions of rotation of a free asymmetric top. While solving Euler’s equations for angular velocity or angular momentum in the body frame we need Jacobi elliptic functions \cn,\sn,\dn, solving the differential equation for the attitude matrix involves EllipticPi function. As I have explained it in Taming the T-handle continued we need the integral

(1)   \begin{equation*}\psi(t)=c_1 t+c_2\int_0^t \frac{1}{1+c_3\,\sn^2(Bs,m)}\,ds.\end{equation*}

In Mathematica this is easily implemented as

(2)   \begin{equation*}\psi(t)=c_1 t+\frac{c_2}{B}\,\Pi(-c_3;\am(Bt,m)|m).\end{equation*}

However, as pointed out by Rowan in a comment to Taming the T-handle continued , the same formula does not work with Maple.

While the documentations of both Mathematica and Maple contain links to Abramowitz and Stegun Handbook of Mathematical Functions, they use different definitions. Here is what concerns us, from p. 590 of the 10th printing:

What we need is 17.2.16, while Maple is using 17.2.14. To convert we need to set x=\sn u, but such a conversion is possible only in the domain where \sn can be inverted. We can do it easily for sufficiently small values of u, but not necessarily for values that contain several quarter-periods.

The following Maple procedure does the job:


epi := proc (t::float, nu::float, k::float) local t2, n, dt, ep0, res; ep0 := EllipticPi(nu, k); t2 := EllipticK(k); n := floor(t/t2); dt := t-t2*n; if type(n, even) then res := Re(n*ep0+EllipticPi(JacobiSN(dt, k), nu, k)) else res := Re((1+n)*ep0-EllipticPi(JacobiSN(t2-dt, k), nu, k)) end if end proc

HAs an example here is the Maple plot for nu=-3, k=0.9:
plot(('epi')(t, -3.0, .9), t = -20 .. 20)

And here is the corresponding Mathematica plot:

The function epi(t,nu, k) defined above for Maple gives now the same result as EllipticPi(nu,JacobiAM(t,k^2),k^2) in Mathematica.