数値的安定性

力学的安定性に加えて、私たちは、数値的安定性についても考慮する必要があります。

これは、私たちのケースでは特に真実です。なぜならば、今は、単純なオイラー法を用いて、回転を表現する四元数の積分を行っているからです。

        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();

図は、四元数qの成分の時間的発展を示しています。初期条件は：

    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;

数値的には、以下のようになります：

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

四元数qは、２秒間で１回転しかしないことに注意して下さい。これに対して、物体に固定した単位Xベクトルは、１秒で１回転します。

さて、４秒間ではなく40秒間のシミュレーションをやってみましょう。

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

まあ、それらしいと言えますかね。