Category: Euler’s equations

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…

Mathematics and sex

Posted on by arkajad

Nowadays mathematics is almost everywhere. To know some cool mathematics (or mathematician) is nowadays sexy.

Therefore I am helping you here, on this blog, to have better relations in your lives. Just looking at the equations can already make you feel better. When you look for a little at Euler’s equations, you will sleep better. That is why today I recall some of these pretty formulas from the dynamics of spinning bodies.

Repetitio est mater studiorum – repetition is the mother of study/learning, that is what they say. So let us repeat, recollect, before traveling another road. If we climb high enough, then going down will be easy.

It is not that climbing up is difficult. It is not. We can take stairs. But it takeS time. And that is what we are going to do now: take time, do it easy way. Enjoying happy sliding down afterwards.

The two essential ingredients af the asymmetric top spinning in zero gravity are:

  1. Euler’s equations
  2. Attitude matrix equations

Euler equations, written in terms of the angular momentum vector \mathbf{L}, with components (L_1,L_2,L_3) are

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

I_1,I_2,I_3 are principal moments of inertia. We assume that our spinning body is asymmetric, therefore I_1,I_2,I_3 are all different from each other, and we will order them as I_1<I_2<I_3. Perhaps it is worthwhile to mention that physics of ordinary materials requires that all three numbers are strictly positive, and that I_1+I_2\geq I_3.  I have my pet rigid body that is essentially flat, with I_1=1,I_2=2,I_3=3. It looks as on the picture below

It is kind of simple but it shows all the nontrivial behavior that a rigid body can have.

Remark: As noted by Bjab in the discussion of this post, the sentence above is incorrect. At the other end of the spectrum of rigid body shapes there are one-dimensional objects. These are flying rods with I_1=0,I_2=I_3\neq 0. They deserve a separate study.

But that is not important now.  What is important is that the components (L_1,L_2,L_3) are the components of the angular momentum vector with respect to the noninertial frame that rotates with the body. Its components with respect to the inertial laboratory frame are constant in time. This is “conservation of angular momentum”- one of the most fundamental “laws” of classical mechanics. Why such a law holds, or at least approximately holds, in our Universe is not known. Physicists and philosophers are debating about “the origin of inertia”.  Some say that “because of space”, some say “because of distant parallelism”, some other say “because of distant matter”. We will not worry about these problems now. We have our laboratory inertial frame, we have frame rotating with the body, and we have rotation matrix Q(t) that maps coordinates in the rotating frame to coordinates in the laboratory frame.

Recall from Angular momentum

Let \mathbf{E}_1(t), \mathbf{E}_2(t), \mathbf{E}_3(t) be an orthonormal frame corotating with the body, and aligned with its principal axes, and let \mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3 be an inertial laboratory frame, both centered at the center of mass of the body. The two frames are related by time-dependent orthogonal matrix A(t)

(2)   \begin{equation*} \mathbf{E}_i=A_{ji}\mathbf{e}_j,\quad X_i= A_{ij}x_j.\end{equation*}

The inverse of A is denoted by Q

(3)   \begin{equation*} Q= A^{-1}=A^t,\end{equation*}

and it is often called the attitude matrix. For a rotating body, if X_i are coordinates of a fixed point in the body, then its coordinates in the laboratory system change in time:

(4)   \begin{equation*} {\bf x}(t)=Q(t){\bf X}.\end{equation*}

Matrix Q(t) satisfies differential equation that is essentially nothing more than the definition of the angular velocity vector:

(5)   \begin{equation*} Q(t)^{-1}\frac{dQ(t)}{dt}=W(\boldsymbol{\Omega}(t)),\end{equation*}

thus

(6)   \begin{equation*} \frac{dQ(t)}{dt}=Q(t)W(\boldsymbol{\Omega}(t)),\end{equation*}

where for any vector \vec{v}

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

We have

(8)   \begin{equation*}L_1=I_1\,\Omega_1,\quad L_2=I_2\,\Omega_2,\quad L_3=I_3\,\Omega_3.\end{equation*}

A couple of comments are due at this point. First of all a careful Reader will notice that in Angular momentum we had

(9)   \begin{equation*}\frac{dQ(t)}{dt}=W(\boldsymbol{\omega}(t))Q(t), \end{equation*}

while now we have

(10)   \begin{equation*}\frac{dQ(t)}{dt}=Q(t)W(\boldsymbol{\Omega}(t)).\end{equation*}

What is going on?
Several things at once. First we have the map \vec{v}\mapsto W(\vec{v}) that associates a matrix to every vector. In fact it associates an antisymmetric matrix. The association has the property that is easy to verify: for every vector \vec{w} we have

    \[W(\vec{v})\,\vec{w}=\vec{v}\times\vec{w}.\]

Acting with the matrix W(\vec{v}) is the same as taking cross product with vector \vec{v}. That such an association should exist should be not a surprise. Taking cross-product with a fixed vector is a linear operation, and every linear operation on vectors is implemented by a matrix. However many students who learn about cross products,

do not learn about their relation to antisymmetric matrices and bivectors. And many students that are learning about matrices, operations with them, eigenvalues and eigenvectors, do not learn about cross-product of vectors. Perhaps because the cross-product is particular to 3D?
Anyway, we have this association, and this association has a very nice property of “covariance”: for every orthogonal matrix R of determinant one we have:

(11)   \begin{equation*}R\,W(\vec{v})\,R^t = W(R\vec{v}).\end{equation*}

That is a very important and very useful property. It can be equivalently expressed as the “covariance of the cross-product”

(12)   \begin{equation*}(R\vec{w})\times (R\vec{v})=R(\vec{w}\times\vec{v}).\end{equation*}

But equivalent expression does not replace the proof, and to prove it needs a bunch of simple but lengthy algebraic calculations and taking into account the definitions of cross-product and determinant. It can be easily done with any computer algebra software. I skip it now.
The second thing is the difference between \boldsymbol{\omega} that represents the angular momentum vector in the laboratory frame, and \boldsymbol{\Omega} that represents the same vector in the body frame. Our matrix Q, by definition, connects the two representations (see Eq. (4) above, but now we skip the time-dependence):

    \[\boldsymbol{\omega}=Q\,\boldsymbol{\Omega}.\]

Now, going back to Eq. (9), we have (again we skip the time dependence)

    \[W(\boldsymbol{\omega})Q=W(Q\boldsymbol{\Omega})Q=Q\,W(\boldsymbol{\Omega})Q^t\,Q=Q\,W(\boldsymbol{\Omega}),\]

which solves our puzzle.
Remark: Above we have used \boldsymbol{\Omega} and \boldsymbol{\omega} to distinguish between numerical representation of the same vector with respect to two different frames. But when there is no possibility of a confusion, we can denote a vector by any convenient letter whatsoever.
Let us see how the above recollection of facts can be applied. In Towards the road less traveled with spin there was the following statement:

The vertical z-axis is the natural axis to try to spin the thing. Imagine our top is floating in space, in zero gravity. We take the z-axis between our fingers, and spin the device. If our hand is not shaking too much, our top will nicely spin about the z-axis. This is the most stable axis for spinning.

The corresponding solution of Euler’s equations is \vec{\omega}=(0,0,\omega_3), where \omega_3 is a constant. The solution of the attitude matrix equation is Q(t)=\exp(W(t\vec{\omega})).

I promised to answer the question why is it so.

In general, when we have matrix differential equation of the type:

    \[\frac{dM(t)}{dt}=M(t)X,\]

where X is a constant (i.e. time independent) matrix, we can verify that M(t)=\exp(tX) is a solution. In fact, it is the unique solution with the initial data M(0)=I. To prove it, one would have to expand e^{tX} into power series, take care about the convergence of the infinite series (no care needed, it is absolutely convergent), and then differentiate term by term. Let us take it for granted that is the case. In the quote from the other post above I used the letters \vec{\omega}=(0,0,\omega_3), instead of “more adequate”, \vec{\Omega}=(0,0,\Omega_3), but, as in the remark above, we are supposed to adjust the meaning to the context. For the matrix X in the attitude equation we have W(\vec{\omega}), so the solution with the property W(0)=I is Q(t)=\exp(tW(\vec{\omega})), which, because \vec{v}\mapsto W(\vec{v}) is a linear map, is the same as Q(t)=\exp(W(t\vec{\omega})).
In fact here we have the opportunity to make another application of today’s recollections. Suppose that our spinning top is completely symmetric – a perfect ball spinning one of its axes. Perfect ball means I_1=I_2=I_3. If so, the right side of the Euler’s equations is automatically zero. Therefore any constant vector \boldsymbol{\Omega} is a solution. Let us choose \boldsymbol{\Omega} a unit vector, so that the absolute value of the angular velocity is 1. That means our top makes a complete rotation about the axis along the vector \boldsymbol{\Omega} every 2\pi units of time. As noticed above the solution of the attitude equation is \exp(tW(\omega})). Therefore \exp(tW(\boldsymbol{\omega})) is the rotation matrix that describes the rotation about the axis in the direction \boldsymbol{\omega} by an angle t. We have already used this fact before, but now it was a good place for recollecting.
After all these recollections we are now ready to look again at the roads that are less traveled.

Of course it is not Albert Einstein who is the author, but, in the age of fake news, who cares?

More is different

Posted on by arkajad

We all know that more is different. Don’t we? But listen to Philip Anderson, who in 1977 was awarded the Nobel Prize in Physics for his investigations into the electronic structure of magnetic and disordered systems, which allowed for the development of electronic switching and memory devices in computers. He ends his paper “More is different” in Science journal as follows

The arrogance of the particle physicist and his intensive research may be behind us (the discoverer of the positron said “the rest is chemistry”), but we have yet to recover from that of some molecular biologists, who seem determined to try to reduce everything about the human organism to “only” chemistry, from the common cold and all mental disease to the religious instinct. Surely there are more levels of organization between human ethology and DNA than there are between DNA and quantum electrodynamics, and each level can require a whole new conceptual structure.
In closing, I offer two examples from economics of what I hope to have said. Marx said that quantitative differences become qualitative ones, but a dialogue in Paris in the 1920’s sums it up even more clearly:

FITZGERALD: The rich are different from us.
HEMINGWAY: Yes, they have more money.

Things are not simple, they are complex, sometimes very hard to understand. We need skills and tricks, intelligence and luck, we need to find our creative muse, which is not an easy task, and even with all that, understanding complexity is often a very complex task.

Now, here we are, with our asymmetric spinning top – the hardest concept in all physics – and there are questions that we need to answer. Like for instance the question asked by Bjab in a comment to the last post Crack in the Cosmic egg:

“But another question is interesting?
Can white point move along a curve in the northern hemisphere in the opposite direction?”

To which I replied (without any explanation whatsoever)

“Interesting question. Yes and no. But why yes and why no? Not so trivial.”

We need to answer those questions that can be answered. Then we will see in full light those questions that are still without answers. So, let’s do it. More is different, so now I go for more.

  1. Angular velocity vector \vector{\Omega} of our free asymmetric top, in the body frame, satisfies Euler’s equations

    (Euler)   \begin{eqnarray*}I_1\dot{\Omega_1}&=&(I_2-I_3)\Omega_2\Omega_3,\\ I_2\dot{\Omega_2}&=&(I_3-I_1)\Omega_3\Omega_1,\\ I_3\dot{\Omega_3}&=&(I_1-I_2)\Omega_1\Omega_2. \end{eqnarray*}

    From the very structure of Euler’s equations we see that if the three components of \vector{\Omega}, (\Omega_1(t),\Omega_2(t),\Omega_3(t)), satisfy Euler’s equations, and we change the sign of any two of them, then we will have another solution of Euler’s equations. The same concerns angular momentum components L_i=I_i\Omega_i,\,(i=1,2,3). That is because either both of these changed components are on the right, then the product of two minuses cancels out giving plus, or one is on the left and one is on the right, and again the changes cancel out.
    Now, changing two components of a vector is equivalent to 180 degrees rotation about the third axis.

  2. Let us have a look at the first animation from the previous post:

    Here d=0.4999 and the code for the solution has these lines:

    A1 = Sqrt[I1*(d*I3 – 1)/(I3 – I1)];
    A2 = Sqrt[I2*(d*I3 – 1)/(I3 – I2)];
    A3 = Sqrt[I3*(1 – d*I1)/(I3 – I1)];
    m = ((d*I3 – 1)*(I2 – I1))/((1 – d*I1)*(I3 – I2))
    m1[t_] = A1*JacobiCN[B*t, m];
    m2[t_] = A2*JacobiSN[B*t, m];
    m3[t_] = A3*JacobiDN[B*t, m];

    The vertical axis on the animation is the third axis. A3>0, and the function \dn is always non-negative. That is why the point moves in the northern hemisphere. If we change the signs of the first two components of the angular momentum vector, that is equaivalent to the rotation by 180 degrees around the vertical axis. The motion will be the same, but the point will be on the opposite side of the curve.

  3. Euler’s equations concern \vector{\Omega}(t), whose components are functions of time, but the equations do not contain time t explicitly. Therefore if \vector{\Omega}(t) is a solution, then \vector{\Omega}(t+t_0) is another solution, for any t_0. We can move the origin of the time axis to any point in time. Let us have a look at the table from Periods of Jacobi elliptic functions – Part 1:

    We see that time translation by 2K is the same as changing the signs of the first two components. And this is the same as shifting the point on the trajectory by half of its period, which is 4K.

  4. Changing the sign of components (1,3) or (2,3) takes us to the trajectory on the southern hemisphere. The function \dn is non-negative, so -\dn is non-positive.
  5. Let us look now at the second animation, for d=0.5001.

    Now the parameter m is larger than 1: m>1. One has to be careful now, because the definitions of the Jacobi functions for m>1 are tricky (see The case of inverted modulus – Treading on Tiger’s tail and Jacobi elliptic cn and dn):

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

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

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

    Now it is the first component (that contains \cn, that is \dn) that is always non-negative. That’s why the trajectory is on the right. If we change the signs of 1 and 2 (or 1 and 3) components, we can move to the trajectory on the left etc.

  6. But we never go on the same trajectory in the opposite direction. Why? This morning I received an email from one of the readers. He notices the fact that the angular velocity and angular momentum vectors are, in fact, pseudovectors, and they are “abstract constructs describing reality.” And indeed. Angular momentum, by its very nature is a “bivector”, not a vector. We represent it as vector using a certain convention. And we draw it using certain conventions. In Angular momentum we have seen this piece:

    … matrix W can be written as

    W = W({\bf v})=\begin{pmatrix}0& -v_3&v_2\\v_3&0&-v_1\\-v_2&v_1&0.\end{pmatrix}

    .

    Why this way? Why not opposite signs? And then, why to draw in a right-handed coordinate system? We can’t represent a bivector as a vector without using some conventions about what is right and what is left. And in more than 3 dimensions vectors and bivectors belong to completely different species!

    I know it is getting complex…. But so is life.