Month: November 2016

Apart from the mechanical stability, we must also be concerned with numerical stability.

This is especially true in our case, because, right now, a simple Euler method is employed in integrating the quaternion that represent the rotation of the body.

        Quaterniond qq(0.0, omega.x()/2.0, omega.y()/2.0, omega.z()/2.0);
        qdot =  qq * q;
        q.w()   += qdot.w()   * dt;
        q.vec() += qdot.vec() * dt;
        q.normalize();

The figure shows the time evolution of the four components of the quaternion, q with the initial conditions of:

    Ibody << 1.0,0.0,0.0, 0.0,1.0,0.0, 0.0,0.0,1.0;
    Quaterniond q(AngleAxisd(M_PI*0.0, Vector3d::UnitX()));
    q.normalize();
    L << 0.00, 0.00, 2.0*M_PI;

Numerically, they are:

double dt = 0.001;
for(double t=0; t<4.0; t+= dt) {

     1   1         0   0   0
     2   0.999995  0   0   0.00314158
     3   0.99998   0   0   0.00628312

  3999   0.99998   0   0  -0.00632446
  4000   0.999995  0   0  -0.00318292
  4001   1         0   0  -4.13415e-05

Note that the quaternion, q, rotates only once in two seconds, while the unit X vector, fixed to the body, rotates once in very second.

Now, we do 40sec simulation instead of 4sec.

double dt = 0.001;
for(double t=0; t<40.0; t+= dt) {

So can we say that it is reasonably stable?

Let’s see if your monolith tumbles or not in various situations.

The figure shows the tip of the Vector3d::UnitZ() fixed to the body. You can see that the vector is always up and the rotation is stable.

The initial conditions are:

    Ibody << 3.0,0.0,0.0, 0.0,4.0,0.0, 0.0,0.0,5.0;
    Quaterniond q(AngleAxisd(M_PI*0.03, Vector3d::UnitX()));
    L << 0.00, 0.00, 1.00;

Note that Ix < Iy < Iz, and therefore the rotation around the z-axis should be stable. The body is slightly rotated around the x-axis to give better visual impression. The angular momentum, L is given, and this vector does not change because there is no external force nor torque.

    Ibody << 3.0,0.0,0.0, 0.0,5.0,0.0, 0.0,0.0,4.0;

Now, Ix < Iz < Iy, and the rotation around the z-axis becomes unstable.

Do you remember TMA-2 orbiting Jupiter? A rectangular cuboid whose sides extend in the ratio of 1 : 4 : 9 has the principal moments of the ratio 1^2+4^2:4^2+9^2:9^2+1^2, or 17:97:82.

An object whose three principal moments are all different is called an asymmetrical top. According to the tennis racket theorem, or as is described in Landau and Lifshitz, the rotation around axis with moment of inertia I2 is unstable, while the ones with I1 and I3 are both stable, assuming that I1 > I2 > I3.

By the way, I was wrong in writing Ibody=(1,2,4) in my previous article, because this can never happen. Again from Landau and Lifshitz:

Looks reasonable? Note the values of Ibody (line 19) and L (line 26) in the code.

#include <iostream>
#include <Eigen/Dense>
#include <Eigen/Core>
#include <Eigen/Geometry>
using namespace std;
using namespace Eigen;

Matrix3d    Ibody;
Matrix3d    Ibodyinv;
Quaterniond q;
Quaterniond qdot;
Vector3d    L;
Matrix3d    Iinv;
Matrix3d    R;
Vector3d    omega;

int main(int argc, const char * argv[]) {

    Ibody << 1.0,0.0,0.0, 0.0,2.0,0.0, 0.0,0.0,4.0;
    Ibodyinv = Ibody.inverse();
    cout << "Ibody: "    << endl << Ibody    << endl;
    cout << "Ibodyinv: " << endl << Ibodyinv << endl;

    Quaterniond q(AngleAxisd(M_PI*0.0, Vector3d::UnitZ()));
    q.normalize();
    L << 0.00, 0.00, 1.00;
    
    
    double dt = 0.1;
    
    for(double t=0; t<200.0; t+= dt) {
        R     = q.toRotationMatrix();
        Iinv  = R * Ibodyinv * R.transpose();
        omega = Iinv * L;
        
        Quaterniond qq(0.0, omega.x()/2.0, omega.y()/2.0, omega.z()/2.0);
        qdot =  qq * q;
        
        cout << "t    : " << t << endl;
        cout << "q    : " << q.w() << " " << q.vec().transpose() << endl;
        cout << "L    : " << L.transpose() << endl;
        cout << "R    : " << endl << R << endl;
        cout << "Iinv : " << endl << Iinv << endl;
        cout << "omega: " << omega.transpose() << endl;
        cout << "qdot : " << qdot.w() << " " << qdot.vec().transpose() << endl;

        q.w()   += qdot.w()   * dt;
        q.vec() += qdot.vec() * dt;
        q.normalize();
    }
    
    return 0;
}

% ./MyEigenTest1 | grep "q    :" | colrm 1 6 | cat -n > t
% ./MyEigenTest1 | grep "omega:" | colrm 1 6 | cat -n > u
% gnuplot
gnuplot> plot "t" using 1:2, "t" using 1:5, "u" using 1:2, "u" using 1:3, "u" using 1:4

Now both Ibody and L are modified slightly.

    Ibody << 1.0,0.0,0.0, 0.0,4.0,0.0, 0.0,0.0,2.0;
    L << 0.01, 0.00, 1.00;

The solution seems to be stable.

Eigen is a C++ template library for linear algebra including the manipulation of quaternion.

#include <iostream>
#include <Eigen/Dense>
#include <Eigen/Core>
#include <Eigen/Geometry>
using namespace std;
using namespace Eigen;

int main(int argc, const char * argv[]) {
    Matrix3f A, B;
    A << 1,0,0, 0,2,0, 0,0,4;
    B = A.inverse();
    cout << A << endl;
    cout << B << endl;
    
    Quaternionf q1(AngleAxisf(M_PI/3.0f, Vector3f::UnitZ()));
    cout << q1.w() << ", " << q1.x() << ", " << q1.y() << ", " << q1.z() << endl;
    return 0;
}

1 0 0
0 2 0
0 0 4
   1    0    0
   0  0.5    0
   0    0 0.25
0.866025, 0, 0, 0.5

Unit quaternions offer you a better way to represent orientations and rotations of objects in three dimensions than rotation matrices do.