Numerical stability

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?