剛体の自由回転：数学的定式化
座標系
慣性主軸を座標軸とする剛体固定座標系を $B$ 系として座標軸を $XYZ$ 軸とし、 実験室座標系（慣性系）を $S$ 系として座標軸を $xyz$ 軸とする。
慣性楕円体：
\[
I_X X^2 + I_Y Y^2 + I_Z Z^2 = 1
\]
ただし、
$I_X$, $I_Y$, $I_Z$ はそれぞれ慣性主軸 $X$, $Y$, $Z$ 軸回りの慣性モーメント。
不変平面： \[ \vec r\cdot\vec L = \sqrt{2 K}; \qquad \vec r := (x, y, z), \quad \mbox{$K$ は回転運動エネルギー} \]
シミュレーションで解いた方程式
運動方程式： 剛体に固定された慣性主軸 $XYZ$ に対する 角速度ベクトル $\vec\omega$ の成分の従う運動方程式を Euler方程式という。 外力のない場合には角運動量ベクトル $\vec L$ は一定で、 \[ {d\vec L\over dt} = \left({\partial\vec L\over\partial t}\right)_B + \vec\omega\times\vec L = 0. \] ただし、 \[ \left({\partial\vec L\over\partial t}\right)_B := {d L_X\over dt}\hat e_X + {d L_Y\over dt}\hat e_Y + {d L_Z\over dt}\hat e_Z \] で、$\vec L$ の $B$ 系における見かけの時間変化を表す。 これより、Euler方程式は \begin{equation} \left\{\begin{array}{rl}\displaystyle {d\omega_X\over dt} &\displaystyle ={I_Y-I_Z\over I_X}\,\omega_Y\,\omega_Z, \\ \displaystyle {d\omega_Y\over dt} &\displaystyle ={I_Z-I_X\over I_Y}\,\omega_Z\,\omega_X, \\ \displaystyle {d\omega_Z\over dt} &\displaystyle ={I_X-I_Y\over I_Z}\,\omega_X\,\omega_Y. \end{array}\right. \end{equation} で与えられる。これは、すでに述べたように非線形の方程式なので、 一般の運動は基本モードの重ね合わせで表せない。回転を表す4元数： 上のオイラー方定式は剛体に固定された $B$ 系で記述されているので、 実験室系で見た運動を表すには解を回転して実験室系（ $S$ 系）に 座標変換しなければならない。
剛体固定系と実験室系の間の回転座標変換を表すのには、 オイラー角を用いる方法と四元数を用いる方法がある。 オイラー角を用いる方法のほうが直感的に理解しやすいが、 オイラー角と角速度ベクトルとの間の変換式に特異点があり、 その近傍で数値解法が破綻する。 そのため、数値計算には四元数を用いる方がよい。
ベクトル $\vec a=(a_x, a_y, a_z)$ を $\hat n$ 軸の回りに $\theta$ 回転すると ベクトル $\vec b=(b_x, b_y, b_z)$ となる場合、 その回転を表す四元数 $R_q$ は \[ b_q = R_q\, a_q R_q^{-1}; \qquad R_q = \exp\left[ {\theta\over 2}\,\hat n_q\right] = \cos{\theta\over 2} + \hat n_q \sin{\theta\over 2} \] で与えられる。 ただし、$\hat n$ は単位ベクトルで、 \[ a_q := a_x \hat i + a_y \hat j + a_z \hat k, \quad b_q := b_x \hat i + b_y \hat j + b_z \hat k, \quad n_q := n_x \hat i + n_y \hat j + n_z \hat k. \] ここで、$\hat i$, $\hat j$, $\hat k $ は四元数の虚数単位元で $\hat i^2 = \hat j^2 =\hat k^2 = \hat i\hat j\hat k =-1 $ を満たす。
剛体固定系から実験室系への変換： あるベクトル $\vec a$ の剛体固定系の $XYZ$ 成分を 実験室系の $xyz$ 成分に変換する4元数を $R_q$ とする。 すなわち、 \[ \vec a =: a_x \vec e_x + a_y \vec e_y + a_z \vec e_z =: a_X \vec e_X + a_Y \vec e_Y + a_Z \vec e_Z \] のとき、４元数 $a^S_q$ および $a^B_q$ を \[ a^S_q := a_x \hat i + a_y \hat j + a_z \hat k, \quad a^B_q := a_X \hat i + a_Y \hat j + a_Z \hat k \] と定義して、 \[ a^S_q = R_q a^B_q R_q^{-1} \] とする。
回転角速度ベクトル $ \vec\omega $ の $ B $ 系での $ XYZ $ 成分が与えられている場合、 $R_q$ は \begin{equation} {dR_q\over dt} = {1\over 2}\, R_q \omega^B_q; \qquad \omega^B_q:= \omega_X\hat i + \omega_Y\hat j + \omega_Z\hat k % \tag{A} \end{equation} から求められる。 この式を成分で表すと、 \begin{equation} {d\over dt} \left(\begin{array}{c} R_W \\ R_X \\ R_Y \\ R_Z \end{array}\right) ={1\over 2} \left(\begin{array}{cccc} 0, & -\omega_X, & -\omega_Y, & -\omega_Z \\ \omega_X, & 0, & \omega_Z, & -\omega_Y \\ \omega_Y, & -\omega_Z, & 0, & \omega_X \\ \omega_Z, & \omega_Y, & -\omega_X, & 0 \end{array}\right) \left(\begin{array}{c} R_W \\ R_X \\ R_Y \\ R_Z \end{array}\right) \end{equation} と与えられる。ただし、 \[ R_q := R_W + R_X \hat i + R_Y\hat j + R_Z\hat k \] である。
右辺の $(\omega_X,\omega_Y, \omega_Z)$ に オイラー方程式で求めたものを用いて微分方程式を数値積分し、 各時刻における $R_q$ を求めて、実験室系における物体の方向を計算する。
式(A)の導出：
$\vec r$ を剛体固定系 $B$ に固定された質点の位置ベクトルとし、 時刻ゼロで $\vec r=\vec r^0$ とする。剛体固定系 $B$ ( $XYZ$ 軸)が 慣性系 $S$ ( $xyz$ 軸)に対して角速度 $\vec\omega$ で回転しているとすると、 質点の速度 $\vec v$ は \[ \vec v = {d\vec r\over dt} = \vec\omega\times\vec r \tag{1} \] で与えられる。一方、 \[ r^S_q := x\hat i + y\hat j + z\hat k, \qquad r^0_q := x^0\hat i + y^0\hat j + z^0\hat k \] とすると、 \[ r^S_q = R_q r^0_q R_q^{-1} \] なので、 \begin{eqnarray} v^S_q & := & {dr^S_q\over dt} = {dR_q\over dt}r^0_q R_q^{-1} + R_q r^0_q {d R_q^{-1}\over dt} \\ & = & \Omega_q r^S_q - r^S_q\Omega_q \tag{2} \end{eqnarray} と表される。ここで、 \[ \Omega_q:= {dR_q\over dt}R_q^{-1} = - R_q{dR_q^{-1}\over dt} \tag{3} \] とする。
式(1)を用いて、式(2)の両辺を比較することにより、 \[ \Omega_q = {1\over 2}\,\omega_q^S \] とすればよいことが分かる。 これより式(3)から式(A) \[ {dR_q\over dt} = {1\over 2}\,\omega_q^S R_q = {1\over 2}\, R_q \omega_q^B \] が得られる。 ここで、最後の等式では、 \[ \omega_q^S = R_q \omega_q^B R_q^{-1} \] を用いた。■