Day: February 2, 2017

Solving Euler’s equations

Posted on by arkajad

We are going to solve now Euler’s equations. In order to solve them in an intelligent way we need some preparation.

In “Angular momentum” we have derived Euler’s equations that govern rotations of a free rigid body about its center of mass:

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

Here \Omega_1,\Omega_2,\Omega_3 are the components of the angular velocity vector in the body frame, while I_1,I_2,I_3 are the three principal moments of inertia.

Euler’s equations follow from conservation of the angular momentum vector with respect to space, as seen from the inertial laboratory frame. This is how we derived them in the previous post. The angular momentum vector in the laboratory frame is Q\mathbf{L}, where \mathbf{L} has components (I_1\Omega_1,I_2\Omega_2,I_3\Omega_3). Although we do not know Q (it needs to be found through solving differential equations), we know that it is an orthogonal matrix. And orthogonal matrix preserves the length of the vector. Since the components of the vector Q\mathbf{L} are constants, also its length is a constant. But the length of Q\mathbf{L} is the same as the length of \mathbf{L}, and the length squared L^2 of Q\mathbf{L} is

(1)   \begin{equation*} L^2=I_1^2\Omega_1^2+I_2^2\Omega_2^2+I_3^2\Omega_3^2.\end{equation*}

There is also another important constant of motion – kinetic energy. We denote by 2E the doubled kinetic energy:

(2)   \begin{equation*}2E=I_1\Omega_1^2+I_2\Omega_2^2+I_3\Omega_3^2.\end{equation*}

To see that indeed 2E is a constant of motion, we differentiate 2E with respect to t:

    \[\frac{d}{dt}\,2E=2(I_1\Omega_1\dot{\Omega_1}+I_2\Omega_2\dot{\Omega_2}+I_3\Omega_3\dot{\Omega_3}).\]

Now, we use the Euler’s equations and, surprise-surprise, everything on the right hand side cancels out:

(3)   \begin{multline*}\frac{d}{dt} 2E=2(I1\Omega_1(I_2-I_3)\Omega_2\Omega_3+\\ I_2\Omega_2(I_3-I_1)\Omega_1\Omega_3+\\ I_3\Omega_3(I_1-I_2)\Omega_1\Omega_2)=\\ 2(I_1I_2-I_1I_3+I_2I_3-I_2I_1+I_3I_1-I_3I_2)\Omega_1\Omega_2\Omega_3=0. \end{multline*}

Thus we get \frac{d}{dt}\,2E=0.

So, we have identified two important constants of motion. Whenever a rigid body rotates, it will change its orientation in space, its rotation axis may change its direction, but the absolute value L of the angular momentum given by Eq. (1) is constant in time. Kinetic energy is also constant in time. Yet these two constants are not completely independent. There are restrictions on the possible values of 2E if L^2 is given. To see this let us order I_1,I_2,I_3 so that

(4)   \begin{equation*}I_1\leq I_2\leq I_3.\end{equation*}

We then have:

    \[ I_12E=I_1(I_1\Omega_1^2+I_2\Omega_2^2+I_3\Omega_3^2)=I_1^2\Omega_1^2+I_1I_2\Omega_2^2+I_1I_3\Omega_3^2\leq L^2.\]

For the same reason we have

    \[ I_32E\geq L^2.\]

Let us introduce the constant d defined by

(5)   \begin{equation*}d=\frac{2E}{L^2}.\end{equation*}

Then we have

(6)   \begin{equation*}\frac{1}{I_3}\leq d\leq\frac{1}{I_1}.\end{equation*}

We can now move to solving the Euler’s equations (Euler) – at the top. To this end we compare them with the properties of the derivatives of Jacobi elliptic functions \sn,\cn,\dn that we have derived in “Derivatives of Jacobi elliptic am, sn, cn, dn“:

(7)   \begin{eqnarray*}\frac{d}{du}\mathrm{sn}(u,m)&=&\mathrm{cn}(u,m)\,\mathrm{dn}(u,m),\\ \frac{d}{du}\mathrm{cn}(u,m)&=&-\mathrm{sn}(u,m)\,\mathrm{dn}(u,m),\\ \frac{d}{du}\mathrm{dn}(u,m)&=&-m\mathrm{sn}(u,m)\,\mathrm{cn}(u,m).\end{eqnarray*}

According to our convention (4) on the right hand side of Euler’s equation we have signs (-,+,-). Our guess is that \Omega_2 is proportional to \sn. We do not know whether to assume \Omega_1 proportional to \cn or to \dn? Let us try \Omega_1 proportional \cn and \Omega_3 proportional to \dn. So, here is our “Ansatz”:

(8)   \begin{eqnarray*} \Omega_1(t)&=&A_1\cn(Bt,m),\\ \Omega_2(t)&=&A_2\sn(Bt,m),\\ \Omega_3(t)&=&A_3\dn(Bt,m).\\ \end{eqnarray*}

Question is: can we find A_1,A_2,A_3,B,m such that the Euler’s equation are satisfied automatically, given L and d?

The answer has been found long ago. We can see a solution in the 1894 book by Alfred Cardew Dixon, “The elementary properties of the elliptic functions” . But we want to do it in a modern way. Here, for instance, is REDUCE code that provides us with a solution:

%We derive a solution of Euler’s equations
o1(t):=A1*cn(B*t);
o2(t):=A2*sn(B*t);
o3(t):=A3*dn(B*t);
%To solve Euler’s equation we must have e1=e2=e3=0 where
e1:=(I1*DF(o1(t),t)-(I2-I3)*o2(t)*o3(t))/(dn(b*t)*sn(b*t));
e2:=(I2*DF(o2(t),t)-(I3-I1)*o3(t)*o1(t))/(cn(b*t)*dn(b*t));
e3:=(I3*DF(o3(t),t)-(I1-I2)*o1(t)*o2(t))/(sn(b*t)*cn(b*t));
% Here are properties of sn, cn, dn
For all v let DF(sn(v),v)=cn(v)*dn(v);
For all v let DF(cn(v),v)=-sn(v)*dn(v);
For all v let DF(dn(v),v)=-m*sn(v)*cn(v);
For all v let cn(v)**2=1-sn(v)**2;
For all v let dn(v)**2=1-m*sn(v)**2;
% en2 is twice kinetic energy:
en2:=I1*o1(t)^2+I2*o2(t)^2+I3*o3(t)^2;
% ll2 is square of angular momentum:
ll2:=I1^2*o1(t)^2+I2^2*o2(t)^2+I3^2*o3(t)^2;
OFF NAT;
OUT “H:\Reduce MyFiles\eu1.txt”;
e1;
e2;
e3;
% we use the three Euler’s equations in order to express A1,A2,A3 in terms of i1,i2,i3,b,m
xx:=solve({e1=0,e2=0,e3=0},{A1,A2,A3});
% We get four solutions, with different signs, we choose one with all signs plus:
PART(xx,3);
a1:= part(part(part(xx,3),1),2);
a2:= part(part(part(xx,3),2),2);
a3:= part(part(part(xx,3),3),2);
en2;
ll2;
% now we call angular momentum squared l^2 and we introduce d=en2/ll2.
% we then find m,b expressed through l and d
% again there are two solutions, with different signs. We choose both plus signs:
xxx:=solve({ll2=l^2,en2=l^2*d},{b,m});
b:=part(part(part(xxx,1),1),2);
m:=part(part(part(xxx,1),2),2);
% now we find final expressions for a1,a2,a3
a1:=a1;
a2:=a2;
a3:=a3;
SHUT “H:\Reduce MyFiles\eu1.txt”;
END;

The results spilled out by REDUCE are not yet nicely arranged. After the arrangement, taking into account our conventions, the answer takes the following form:

(9)   \begin{eqnarray*} A_1&=&L \sqrt{\frac{dI_3-1}{I_1(I_3-I_1)}},\\ A_2&=&L\sqrt{\frac{dI_3-1}{I_2(I_3-I_2)}},\\ A_3&=&L \sqrt{\frac{1-dI_1}{I_3(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*}

Now I should verify that it is indeed a solution. I was trying to do it using REDUCE, but I did not succeed. So, I will return to this problem in my next post (and I will also try to fix typos in this post). Hopefully I will be able to find a way of verifying for sure the above statements.

Update: The following REDUCE code verifies the solution on random data. On output it spits out three zeros.:
a1:=l*sqrt((d*i3-1)/(i1*(i3-i1)));
a2:=l*sqrt((d*i3-1)/(I2*(i3-i2)));
a3:=l*sqrt((1-d*i1)/(i3*(i3-i1)));
b:=l*sqrt((1-d*i1)*(I3-i2)/(i1*i2*i3));
m:=(d*i3-1)*(i2-i1)/((1-d*i1)*(i3-i2));
o1(t):=A1*cn(B*t);
o2(t):=A2*sn(B*t);
o3(t):=A3*dn(B*t);
For all v let DF(sn(v),v)=cn(v)*dn(v);
For all v let DF(cn(v),v)=-sn(v)*dn(v);
For all v let DF(dn(v),v)=-m*sn(v)*cn(v);
For all v let cn(v)**2=1-sn(v)**2;
For all v let dn(v)**2=1-m*sn(v)**2;
e1:=(I1*DF(o1(t),t)-(I2-I3)*o2(t)*o3(t))/(dn(b*t)*sn(b*t));
e2:=(I2*DF(o2(t),t)-(I3-I1)*o3(t)*o1(t))/(cn(b*t)*dn(b*t));
e3:=(I3*DF(o3(t),t)-(I1-I2)*o1(t)*o2(t))/(sn(b*t)*cn(b*t));
for all a,b let sqrt(a)*sqrt(b)=sqrt(a*b);
off rounded;
i1:=random(1000)/1000;
i2:=1.0+random(1000)/1000;
i3:=2.0+random(1000)/1000;
l:=random(10)/2;
d:=1/i3+(1/i1-1/i3)/2;

OFF NAT;
OUT “H:\Reduce MyFiles\eu3.txt”;
e1:=e1;
e2:=e2;
e3:=e3;
SHUT “H:\Reduce MyFiles\eu3.txt”;
END;

Update 2:

Both REDUCE and Mathematica are able to do the task in a symbolic mode, but not without human help. These programs have very smart simplification algorithms, but sometimes our eye see that two expressions are the same, though they have different forms, while the programs do not want to see it. Then we need to help the program. For instance with Mathematica, it did not want to see that a particular expression is zero, but it recognized that its square is zero! Like here:


Here is the Mathematica notebook verifying Euler’s equations.
So, I consider that the job is done, and that we do have a solution of Euler’s equations.