Ham Radio Blog
画像をクリックすると、リアルタイム・アニメーションが見れます。
もう１つの単純な場合は、I1=I2≠I3であるときです。またもや、Landau and Lifshitzから。
これは、MIT OpenCourseWare Courseからの素敵な図です。
シミュレーションのパラメータは：
Ibody.set(1,0,0, 0,1,0, 0,0,2); L.set(00, 01, 2.0*Math.PI);
従って、数式(36.6)のωの値はπに、また、A (=√ Ω1^2 + Ω2^2)の値は1.0になります。グラフのwx, wy, wzの曲線を見て下さい。
幾つかの特別な場合について見てみましょう。再び、Landau and Lifshitzから：
明らかに、最も単純な場合は、I1=I2=I3であるときです。この場合、Omegaは定数ベクトルになります。ちょうど、図に示されているように。用いられたパラメータは：
Ibody.set(1,0,0, 0,1,0, 0,0,1); L.set(1, 2, 3);
もし、
L.set(1, 1, 1);
とすれば、以下のようになります。
２つの図で、Y軸のスケールが大幅に異なることに注意して下さい。
画像をクリックすると、リアルタイム・アニメーションが見れます。
パラメータは、以下のとおりです。
Ibody.set(1,0,0, 0,1,0, 0,0,2); L.set(2.0*Math.PI*0.0,2.0*Math.PI*0.0,2.0*Math.PI*5.0); var dt = 1.0/1000.0;
物体は毎秒2.5回だけ回転すべきで、シミュレーションでは正にそうなっています。グラフは、単位四元数の成分を示していますが、これは毎秒1.25回だけ回転します。
リアルタイム・アニメーションでは、物体はゆっくりと回転していますが、これはフレームレートが30fpsしかないからです。
どちらかをクリックしてみて下さい。非安定な場合、もしくは、安定な場合。
ソースコードは誤りを含んでいるかもしれませんが、まあ、基本的にはOKなのでしょう。私は、matrix4_from_matrix3()とかmatrix3_from_matrix4()と言った関数を自分で書く必要があるのか良く分かりません。
疑問であれば、下記のドキュメントを調べてみて下さい。
https://threejs.org/docs/api/math/Vector3.html
https://threejs.org/docs/api/math/Quaternion.html
https://threejs.org/docs/api/math/Matrix3.html
https://threejs.org/docs/api/math/Matrix4.html
<html>
<head>
<title>My first Three.js app</title>
<style>
body { margin: 0; }
canvas { width: 100%; height: 100% }
</style>
</head>
<body>
<script src="https://ajax.googleapis.com/ajax/libs/threejs/r76/three.min.js"></script>
<script>
function matrix4_from_matrix3 (m4, m3) {
var se = m3.elements;
var te = m4.elements;
te[0] = se[0]; te[4] = se[3]; te[ 8] = se[6]; te[12] = 0;
te[1] = se[1]; te[5] = se[4]; te[ 9] = se[7]; te[13] = 0;
te[2] = se[2]; te[6] = se[5]; te[10] = se[8]; te[14] = 0;
te[3] = 0; te[7] = 0; te[11] = 0; te[15] = 0;
}
function matrix3_from_matrix4 (m3, m4) {
var se = m4.elements;
var te = m3.elements;
te[0] = se[0]; te[3] = se[4]; te[6] = se[ 8];
te[1] = se[1]; te[4] = se[5]; te[7] = se[ 9];
te[2] = se[2]; te[5] = se[6]; te[8] = se[10];
}
var scene = new THREE.Scene();
var camera = new THREE.PerspectiveCamera(7, window.innerWidth/window.innerHeight, 0.1, 1000 );
var renderer = new THREE.WebGLRenderer();
renderer.setSize( window.innerWidth, window.innerHeight );
document.body.appendChild( renderer.domElement );
var geometry = new THREE.BoxGeometry( 4, 9, 1 );
var material2= new THREE.MeshLambertMaterial( { color: 0xf01010 } );
var cube= new THREE.Mesh( geometry, material2 );
cube.position.x = 0.0;
scene.add( cube );
var geometry0= new THREE.BoxGeometry( 1000, 1000, 1);
var material0 = new THREE.MeshBasicMaterial( { color: 0x101040 } );
var cube0 = new THREE.Mesh( geometry0, material0 );
cube0.position.z = -10;
scene.add( cube0 );
camera.position.x = 100;
camera.position.y = 100;
camera.position.z = 250;
camera.lookAt(cube.position);
scene.add( camera );
var light = new THREE.DirectionalLight(0xffffff);
light.position.set(250, 250, 250);
scene.add( light );
var light0 = new THREE.AmbientLight( 0x404040 );
scene.add( light0 );
var Ibody = new THREE.Matrix3();
var IbodyInv = new THREE.Matrix3();
var IbodyInv4 = new THREE.Matrix4();
Ibody.set(82,0,0, 0,17,0, 0,0,97);
IbodyInv.getInverse(Ibody);
matrix4_from_matrix3 (IbodyInv4, IbodyInv);
var q = new THREE.Quaternion();
var qq = new THREE.Quaternion();
var qdot = new THREE.Quaternion();
var axis = new THREE.Vector3(0,1,0).normalize();
var angle = Math.PI * 0.0;
q.setFromAxisAngle(axis,angle);
var x, y, z, w;
var L = new THREE.Vector3();
var L4 = new THREE.Vector4();
L.set(13.0*2.0*Math.PI*82.0,2.0*Math.PI*0.1,2.0*Math.PI*0.1);
L4.set(L.x, L.y, L.z, 0);
var omega = new THREE.Vector3();
var R = new THREE.Matrix3();
var RInv = new THREE.Matrix3();
var IInv = new THREE.Matrix3();
var R4 = new THREE.Matrix4();
var RInv4 = new THREE.Matrix4();
var IInv4 = new THREE.Matrix4();
var dt = 0.001;
var render = function () {
requestAnimationFrame( render );
x = q.x; y = q.y; z = q.z; w = q.w;
R.set(1.0-2.0*y*y-2.0*z*z, 2.0*x*y+2.0*w*z, 2.0*x*z-2.0*w*y,
2.0*x*y-2.0*w*z, 1.0-2.0*x*x-2.0*z*z, 2.0*y*z+2.0*w*x,
2.0*x*z+2.0*w*y, 2.0*y*z-2.0*w*x, 1.0-2.0*x*x-2.0*y*y);
RInv.getInverse(R);
matrix4_from_matrix3 (R4 , R );
matrix4_from_matrix3 (RInv4, RInv);
// IInv = R * IbodyInv * RInv;
IInv4.multiplyMatrices(R4, IbodyInv4);
IInv4.multiply (RInv4);
matrix3_from_matrix4 (IInv, IInv4);
// omega = IInv * L;
omega.copy(L);
omega.applyMatrix3(IInv);
qq.set(omega.x/2.0, omega.y/2.0, omega.z/2.0, 0.0);
// qdot = qq * q;
qdot.multiplyQuaternions(qq, q);
// q += qdot * dt;
q.set(q.x+qdot.x*dt, q.y+qdot.y*dt, q.z+qdot.z*dt, q.w+qdot.w*dt);
q.normalize();
cube.rotation.setFromQuaternion(q);
renderer.render(scene, camera);
}
render();
</script>
</body>
</html>
画像をクリックすると、リアルタイム・アニメーションが見れます。
あなたが、宇宙空間で３つのモノリスを見つけたとしましょう。それぞれは、３つの主慣性軸のいずれかの回りに回転しています。あなたは、どのモノリスが安定であるために能動的に制御されているかが解りますか？
力学的安定性に加えて、私たちは、数値的安定性についても考慮する必要があります。
これは、私たちのケースでは特に真実です。なぜならば、今は、単純なオイラー法を用いて、回転を表現する四元数の積分を行っているからです。
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) {
まあ、それらしいと言えますかね。
あなたのモノリスが宙返りをするかどうかを、色んな状況で調べてみましょう。
図は、物体に固定されたVector3d::UnitZ()の先端を示しています。あなたは、ベクトルが常に上を向いていて、回転が安定であることが見えるでしょう。
初期条件は：
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;
Ix < Iy < Izなので、ｚ軸回りの回転は安定であることに注意して下さい。良い視覚的効果を与えるために、物体はｘ軸の周りで少しだけ回転されています。角運動量Lが与えられ、このベクトルは外力やトルクが無いので変化しません。
Ibody << 3.0,0.0,0.0, 0.0,5.0,0.0, 0.0,0.0,4.0;
今度は、Ix < Iz < Iyにしたので、ｚ軸回りの回転は非安定になりました。
あなたは木製軌道上のTMA-2を覚えていますか？辺の長さの比が1 : 4 : 9の直方体の主慣性モーメントの比は1^2+4^2:4^2+9^2:9^2+1^2、すなわち、17:97:82になります。
主慣性モーメントの値が全て異なる物体は、非対称コマと呼ばれます。テニスラケット定理、もしくは、Landau and Lifshitzに書かれているように、慣性モーメントがI2である軸の周りの回転は安定ではないのに対して、I1もしくはI3周りの回転は安定です。ここで、I1 > I2 > I3と仮定していますが。
ところで、以前の記事で私が、Ibody=(1,2,4)としたのは、間違いでした。なぜなら、このようなことは決して起きないからです。再び、Landau and Lifshitzから引用すると：
それらしく見えるでしょうか？19行のIbodyと26行のLの値に注意して下さい。
#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
今度は、IbodyとLの両方を少し変えました。
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;
解は安定なようです。
Eigenは、四元数の操作も含む線型代数のC++テンプレートライブラリです。
#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
単位四元数は、３次元空間内での物体の向きと回転とを表現するのに回転行列よりも優れた方法を提供します。