数値積分に用いた運動方程式は、
\begin{align*}
\big( I_X + Ma^2 \big) \left({ d\omega_X\over d t}
\textcolor{red}{+\alpha|\dot\psi|\omega_X}\right)
& =
- \big( I_Z +Ma^2\big) \omega_Y\omega_Z {-}I_X \omega_Y^2 \tan {\theta}
{+} Ma^2 {g\over a} \sin {\theta}
\\%--------------
I_X \left( { d\omega_Y\over d t}
\textcolor{red}{+\alpha|\dot\psi|\omega_Y}\right)
& =
{+}I_X\omega_X \omega_Y\tan {\theta} + I_Z\omega_Z \omega_X
\\%--------------
\big( I_Z + Ma^2 \big)\left( { d\omega_Z\over d t}
\textcolor{red}{+\alpha|\dot\psi|\omega_Z}\right)
& =
+ Ma^2 \omega_X \omega_Y
\\%------------------------------
{d\theta\over dt} & = \omega_X
\\
{d\phi\over dt} & = {1\over\cos\theta}\,\omega_Y
\\
{d\psi\over dt} & = \omega_Z + \tan\theta\,\omega_Y
\end{align*}
で、これを4次のルンゲ・クッタ法で数値積分した。
この導出はノートにあります。
赤字で示した項は、現象論的に導入した散逸項で、散逸のない場合は
$\alpha=0$。
