はじめに

「テニスラケットの定理」って知ってますか？

テニスラケットの定理とは、3軸の慣性モーメントが異なる物体の回転運動を考えた際に、第二慣性主軸（慣性モーメントが大中小となるように並べたときの二番目の慣性主軸）周りの運動は不安定になるという定理のことです。

この定理による不思議な運動はジャニベコフ効果とも呼ばれ、「NASA Dzhanibekov effect」で検索すると無重力環境でT字形の物体（T形レンチか何か？）を回している動画が見つかると思います。見たことがない方はまず調べてみて下さい。

動画の中で起きている不思議な運動は空力的な効果（乱流とか）とは関係なく発生するので、本記事では剛体の回転運動シミュレーションを行い、テニスラケットの定理（とそれにより発生する不思議な運動）を確認します。

シミュレーションで回す物体は何でも良いのですが、折角なのでNASAの動画のものと同じ形状にします。適当に寸法を決めて以下のようにモデリングしました。

シミュレーションにおいては、この物体を空気抵抗や重力といった外力が全く働かない状況で回転させます。

剛体の回転運動方程式

回転運動に限らず、物体の運動はダイナミクス（Dynamics）とキネマティクス（Kinematics）に分けて考えることができます。一般に、ダイナミクスは外力と速度の関係を表し、キネマティクスは速度と位置の関係を表します。両者の関係については文献[1]に以下のような図が掲載されています。

＜DynamicsとKinematicsの関係（文献[1]中の図6を参考に作成）＞

図中の上段が並進運動の場合、下段が回転運動の場合のダイナミクスとキネマティクスの関係を表しています。今回は回転運動のみのシミュレーションなので直接関係があるのは下段だけですね。

具体的な式展開を行う前に、慣性空間に固定された慣性座標系($X\text{-}Y\text{-}Z$)と原点が物体の重心位置に固定され物体と共に回転する物体座標系($x\text{-}y\text{-}z$)を考えます。

＜慣性座標系と物体座標系＞

また、式中の変数がどちらの座標系上の値であるかを上付きの添字で表します。例えば慣性座標系上の角速度なら$\boldsymbol{\omega}^\mathrm{i}$とし、物体座標系上の角速度なら$\boldsymbol{\omega}^\mathrm{b}$とします。

キネマティクス

説明の都合上キネマティクスの式を先に示します。 回転運動の場合には姿勢表現の方法が複数あり、今回行うような物体がくるくる回るシミュレーションで特異点があると非常に厄介なので姿勢表現には四元数を用います。 本記事中においては、四元数を表す変数は上にチルダを付けて表記します。

物体の姿勢を表す四元数を$\tilde{q} \in \mathbb{H}$とすると、$\tilde{q}$の時間微分は以下のようになります[2, 3]。

$$\frac{d\tilde{q}}{dt} = \frac{1}{2} \tilde{q} \boldsymbol{\omega}^\mathrm{b} \tag{1}$$

ここで$\boldsymbol{\omega}^\mathrm{b} \mathrm{[rad/s]} \in \mathbb{R}^3$は物体座標系上における角速度です。最終的にはこの微分方程式を解くことで物体の回転運動を計算します。

ダイナミクス

※ 追記(2024-05-05)：ダイナミクスの式が間違っていたので修正しました。式(2)から式(3)の導出についてはページ下部の補足をご覧ください。

トルクを $\boldsymbol{\tau} \mathrm{[N \, m]} \in \mathbb{R}^3$ 、角運動量を $\boldsymbol{L} \mathrm{[N \, m \, s]} \in \mathbb{R}^3$ とすると剛体のダイナミクスは以下の式で表されます[4]。

$$\begin{align*} \boldsymbol{\tau}^\mathrm{i} &= \frac{d \boldsymbol{L}^\mathrm{i}}{dt} \tag{2} \\ \boldsymbol{\tau}^\mathrm{b} &= \frac{d \boldsymbol{L}^\mathrm{b}}{dt} + \boldsymbol{\omega}^\mathrm{b} \times \boldsymbol{L}^\mathrm{b} \tag{3} \end{align*}$$

これが一般に「オイラーの運動方程式」と呼ばれる回転の運動方程式です（$\times$は外積）。ここで、物体の慣性テンソルを$I$とすると物体座標系上における角速度と角運動量の間には

$$\boldsymbol{L}^\mathrm{b} = I \boldsymbol{\omega}^\mathrm{b} \tag{4}$$

の関係が成り立ち、この関係と式(3)より以下の式が得られます。

$$\frac{d \boldsymbol{\omega}^\mathrm{b}}{dt} = I^{-1} ( \boldsymbol{\tau}^\mathrm{b} – \boldsymbol{\omega}^\mathrm{b} \times I \boldsymbol{\omega}^\mathrm{b} ) \tag{5}$$

上式は$\boldsymbol{\omega}^\mathrm{b}$に関する微分方程式ですから、これを数値的に解けばキネマティクスに必要な角速度が得られます。 しかし、この式を用いて実際にシミュレーションを行うと（特にオイラー法では）徐々に角運動量が変化してしまいます。 数値計算的な誤差は仕方ないとはいえ、外力が働かない場合には慣性座標系上における角運動量$\boldsymbol{L}^\mathrm{i}$は一定となるはずなのでこの結果は望ましくありません。

そこで、式(2)と式(4)より角運動量を陽に扱います。 具体的には、式(2)より慣性座標系上で$\boldsymbol{L}^\mathrm{i}$を更新し、次に式(4)を用いて角速度を求めます。 ただし、式(4)における角運動量は物体座標系上の値ですから、四元数を用いて物体座標系上の値に変換します。 慣性座標系から物体座標系への座標変換は共役四元数$\tilde{q}^*$を用いて

$$\boldsymbol{L}^\mathrm{b} = \tilde{q}^* \boldsymbol{L}^\mathrm{i} \tilde{q} \tag{6}$$

と書くことができるので[5]、式(4)より物体座標系上における角速度が以下のように求まります。

$$\begin{align*} \boldsymbol{\omega}^\mathrm{b} &= I^{-1} \boldsymbol{L}^\mathrm{b} \tag{7} \\ &= I^{-1} ( \tilde{q}^* \boldsymbol{L}^\mathrm{i} \tilde{q} ) \tag{8} \end{align*}$$

これにより、角運動量保存則を成立させた状態で角速度を求めることができるようになりました。

数値解法

ここまでに示した微分方程式を数値解法で解いていきます。

まずは角運動量を更新するために式(2)を解く必要がありますが、今回は外力が加わらないので$\boldsymbol{L}^\mathrm{i}$は常に一定となります。 つまり、初期時刻$t = 0$において慣性座標系と物体座標系が一致しているものとすれば

$$\boldsymbol{L}_{(0)}^\mathrm{i} = I \boldsymbol{\omega}_{(0)}^\mathrm{b} \tag{9}$$

となるので、シミュレーション開始時に$\boldsymbol{L}_{(0)}^\mathrm{i}$を計算しておくだけで済みます。仮に外力が加わる場合でも台形公式などを使って数値積分すれば計算できます。

$\boldsymbol{L}^\mathrm{i}$が得られれば式(8)より$\boldsymbol{\omega}^\mathrm{b}$も簡単に計算できます。 $\boldsymbol{\omega}^\mathrm{b}$が得られたらあとは式(1)を解くだけですが、式(8)には姿勢$\tilde{q}$が含まれているので数値解法で解くために式(8)を式(1)に代入します。

$$\begin{align*} \frac{d\tilde{q}}{dt} &= \frac{1}{2} \tilde{q} \{ I^{-1} ( \tilde{q}^* \boldsymbol{L}^\mathrm{i} \tilde{q} ) \} \tag{10} \\ &= f(I, \boldsymbol{L}^\mathrm{i}, \tilde{q}) \tag{11} \end{align*}$$

これで$\tilde q$に関する微分方程式が得られたので、上式をルンゲ・クッタ法[6]で解くことで物体の回転運動を計算することができます。

ルンゲ・クッタ法では、現時刻における物体の姿勢を$\tilde{q}_{(t)}$、刻み幅$h$だけ進んだ時刻の姿勢を$\tilde{q}_{(t + h)}$とすると、更新のための変数

$$\begin{align*} \tilde{f_1} &= f(I, \boldsymbol{L}^\mathrm{i}, \tilde{q}_{(t)}) \tag{12} \\ \tilde{f_2} &= f(I, \boldsymbol{L}^\mathrm{i}, \tilde{q}_{(t)} + h \tilde{f_1} / 2) \tag{13} \\ \tilde{f_3} &= f(I, \boldsymbol{L}^\mathrm{i}, \tilde{q}_{(t)} + h \tilde{f_2} / 2) \tag{14} \\ \tilde{f_4} &= f(I, \boldsymbol{L}^\mathrm{i}, \tilde{q}_{(t)} + h \tilde{f_3}) \tag{15} \end{align*}$$

より、各計算ステップにおける更新式が以下のように得られます。

$$\tilde{q}_{(t + h)} = \tilde{q}_{(t)} + \frac{h}{6} \left(\tilde{f_1} + 2 \tilde{f_2} + 2 \tilde{f_3} + \tilde{f_4} \right) \tag{16}$$

ただし、初期時刻における姿勢は$\tilde{q}_{(0)} = 1$（恒等四元数）とします。

実装

行列演算や四元数演算がやりやすいのでシミュレーション部分はRustで実装します。

[package]
name = "tennis_racket_theorem"
version = "0.1.0"
edition = "2021"

[dependencies]
quaternion-core = "0.3.4"
nalgebra = "0.31"

use std::fs;
use std::io::{Write, BufWriter};
use quaternion_core as quat;
use quat::Quaternion;
use nalgebra::{SVector, SMatrix};

const DT: f64 = 0.03125;  // 32Hz
const TIME_RANGE: f64 = 10.0;  // [s]
const SPLIT_NUM: usize = (TIME_RANGE / DT) as usize + 1;
const OUT_FILE_PATH: &'static str = "./result.csv";

type SVector3 = SVector<f64, 3>;
type SMatrix3x3 = SMatrix<f64, 3, 3>;

fn main() {
    // シミュレーション結果の保存先（同一ファイルが存在したら上書き）
    let mut file = BufWriter::new(fs::File::create(OUT_FILE_PATH).unwrap());

    // トルク
    let torque = SVector3::zeros();  // 外力無し
    // 初期角速度[rad/s]
    let gyr_init = SVector3::new(0.01, 8.0, 0.01);
    // 慣性テンソル（慣性モーメントの単位は[kg m^2]）
    let inertia_tensor = 1e-6 * SMatrix3x3::from_diagonal(&SVector3::new(62.2, 171.5, 210.5));
    // 慣性座標系上における角運動量
    let mut inertial_angular_momentum = inertia_tensor * gyr_init;
    // 四元数
    let mut q: Quaternion<f64> = (1.0, [0.0; 3]);

    let mut pre_gyr = gyr_init;
    for i in 0..SPLIT_NUM {
        // 表示用
        let body_angular_momentum: SVector3 = quat::frame_rotation(q, inertial_angular_momentum.into()).into();
        let body_gyr = inertia_tensor.try_inverse().unwrap() * body_angular_momentum;
        let alpha = (body_gyr - pre_gyr) / DT;  // 角加速度
        pre_gyr = body_gyr;

        // ダイナミクス
        inertial_angular_momentum = inertial_angular_momentum + DT * torque;

        // キネマティクス
        q = rk4(inertia_tensor, inertial_angular_momentum, q);
        q = quat::normalize(q);

        // シミュレーション結果保存
        let ypr = quat::to_euler_angles(quat::RotationType::Intrinsic, quat::RotationSequence::ZYX, q);
        file.write(format!(
            "{:.4},{:.4},{:.4},{:.4},{:.4},{:.4},{:.4},{:.4},{:.4},{:.4},{:.7},{:.7},{:.7}\n",
            DT * i as f64, ypr[0], ypr[1], ypr[2],
            body_gyr[0], body_gyr[1], body_gyr[2],
            alpha[0], alpha[1], alpha[2],
            body_angular_momentum[0], body_angular_momentum[1], body_angular_momentum[2]
        ).as_bytes()).unwrap();
    }

    println!("Result has been saved to {}", OUT_FILE_PATH);
}

/// 4次のルンゲ・クッタ法
fn rk4(inertia_tensor: SMatrix3x3, inertial_angular_momentum: SVector3, q: Quaternion<f64>) -> Quaternion<f64> {
    let f1 = dq_dt(inertia_tensor, inertial_angular_momentum, q);
    let f2 = dq_dt(inertia_tensor, inertial_angular_momentum, quat::scale_add(DT*0.5, f1, q));
    let f3 = dq_dt(inertia_tensor, inertial_angular_momentum, quat::scale_add(DT*0.5, f2, q));
    let f4 = dq_dt(inertia_tensor, inertial_angular_momentum, quat::scale_add(DT, f3, q));
    let tmp1 = quat::scale_add(2.0, f2, f1);
    let tmp2 = quat::scale_add(2.0, f3, f4);
    quat::scale_add(DT/6.0, quat::add(tmp1, tmp2), q)
}

/// 四元数の時間微分
fn dq_dt(inertia_tensor: SMatrix3x3, inertial_angular_momentum: SVector3, q: Quaternion<f64>) -> Quaternion<f64> {
    let body_angular_momentum: SVector3 = quat::frame_rotation(q, inertial_angular_momentum.into()).into();
    let body_gyr = inertia_tensor.try_inverse().unwrap() * body_angular_momentum;
    quat::scale(0.5, quat::mul(q, (0.0, body_gyr.into())))
}

四元数の演算にはquaternion-coreクレートを使用しています。 結構使いやすいと思うのでぜひ使って下さい（自分が作ったクレートなので使いやすいと感じるのは当然かもしれませんが…）。

シンタックスハイライトがRustに対応していなくて見にくいのでGithubにも上げています。

https://github.com/HamaguRe/tennis_racket_theorem

シミュレーション結果のグラフを表示するPythonスクリプトとアニメーション表示用のProcessingコードもリンク先で見れます。

シミュレーション結果

初期時刻における角速度を$\boldsymbol{\omega}_{(0)}^\mathrm{b} = [0.01, 8.0, 0.01]$としてシミュレーションを行いました。$y$軸以外の角速度を$0$にしてしまうと全く遷移しなくなるので小さい値を設定しています。

シミュレーション結果はCSVファイルに保存され、グラフにすると以下のようになります。 グラフだけ見てもよくわからないかもしれませんが、物体座標系の$y$軸周りの角速度（上から2段目）を見ると何か規則的な動きをしていることが確認できると思います。

＜シミュレーション結果（角速度・角加速度・角運動量は全て物体座標系上の値）＞

0秒から10秒の範囲をアニメーション表示すると以下のようになります。

Wow! NASAの動画と同じような挙動になっていますね！

おわりに

誰か宇宙に行く機会があったら本当にシミュレーション通りの挙動を示すか確認してみて下さい。

参考文献

1. 山口 功，木田 隆，岡本 修，狼 嘉彰，“クォータニオンとオイラー角によるキネマティックス表現の比較について”，航空宇宙技術研究所資料，NAL TM-636，航空宇宙技術研究所，1991．
2. 堀 源一郎，『ハミルトンと四元数：人・数の体系・応用』，海鳴社，2007．
3. 矢田部 学，“クォータニオン計算便利ノート”，MSS技報，Vol.18，pp.29-34，三菱スペース・ソフトウェア（現 三菱電機ソフトウェア），2007．
4. ファインマン，レイトン，サンズ，『ファインマン物理学 I』（坪井忠二訳），岩波書店，1967．
5. 宇宙電波実験室，“四元数まとめ資料（2022年8月27日版）”，2022．
6. 長沢 工，桧山澄子，『パソコンで見る天体の動き』，地人書館，1992．

補足ー四元数を用いた式(3)の導出ー

慣性座標系に対する物体座標系の姿勢$\tilde q$が与えられた場合、両座標系におけるトルク・角運動量・角速度の関係は以下の式で表されます。

$$\begin{align*} \boldsymbol{\tau}^\mathrm{b} &= \tilde{q}^* \boldsymbol{\tau}^\mathrm{i} \tilde{q} \\ \boldsymbol{L}^\mathrm{b} &= \tilde{q}^* \boldsymbol{L}^\mathrm{i} \tilde{q} \\ \boldsymbol{\omega}^\mathrm{b} &= \tilde{q}^* \boldsymbol{\omega}^\mathrm{i} \tilde{q} \end{align*}$$

上記の関係を式(2)に代入して展開すると

$$\begin{align*} \tilde{q} \boldsymbol{\tau}^\mathrm{b} \tilde{q}^* &= \frac{d}{dt} (\tilde{q} \boldsymbol{L}^\mathrm{b} \tilde{q}^*) \\ &= \tilde{q}^\prime (\boldsymbol{L}^\mathrm{b} \tilde{q}^*) + \tilde{q} (\boldsymbol{L}^\mathrm{b} \tilde{q}^*)^\prime \\ &= \tilde{q}^\prime (\boldsymbol{L}^\mathrm{b} \tilde{q}^*) + \tilde{q} ({\boldsymbol{L}^\mathrm{b}}^\prime \tilde{q}^* + \boldsymbol{L}^\mathrm{b} {\tilde{q}^*}^\prime) \\ &= \frac{1}{2} \tilde{q} \boldsymbol{\omega}^\mathrm{b} \boldsymbol{L}^\mathrm{b} \tilde{q}^* + \tilde{q} \frac{d \boldsymbol{L}^\mathrm{b}}{dt} \tilde{q}^* + \tilde{q} \boldsymbol{L}^\mathrm{b} \left( \frac{1}{2} \tilde{q} \boldsymbol{\omega}^\mathrm{b} \right)^* \quad \left( \because \frac{d \tilde{q}^*}{dt} = \left( \frac{d \tilde{q}}{dt} \right)^* \right) \\ &= \frac{1}{2} \tilde{q} \boldsymbol{\omega}^\mathrm{b} \boldsymbol{L}^\mathrm{b} \tilde{q}^* + \tilde{q} \frac{d \boldsymbol{L}^\mathrm{b}}{dt} \tilde{q}^* + \frac{1}{2} \tilde{q} \boldsymbol{L}^\mathrm{b} {\boldsymbol{\omega}^\mathrm{b}}^* \tilde{q}^* \\ &= \frac{1}{2} \tilde{q} \boldsymbol{\omega}^\mathrm{b} \boldsymbol{L}^\mathrm{b} \tilde{q}^* + \tilde{q} \frac{d \boldsymbol{L}^\mathrm{b}}{dt} \tilde{q}^* - \frac{1}{2} \tilde{q} \boldsymbol{L}^\mathrm{b} \boldsymbol{\omega}^\mathrm{b} \tilde{q}^* \qquad (\because {\boldsymbol{\omega}^\mathrm{b}}^* = -\boldsymbol{\omega}^\mathrm{b}) \end{align*}$$

となります。ただし、見易さのため変数$x$の時間微分を$x^\prime$と書きました。 続いて、単位四元数では$\tilde{q}^* = \tilde{q}^{-1}$であることを利用して両辺に左右からそれぞれ$\tilde{q}^*$と$\tilde{q}$を掛け整理すると

$$\begin{align*} \boldsymbol{\tau}^\mathrm{b} &= \frac{1}{2} \boldsymbol{\omega}^\mathrm{b} \boldsymbol{L}^\mathrm{b} + \frac{d \boldsymbol{L}^\mathrm{b}}{dt} - \frac{1}{2} \boldsymbol{L}^\mathrm{b} \boldsymbol{\omega}^\mathrm{b} \\ &= \frac{1}{2} \boldsymbol{\omega}^\mathrm{b} \boldsymbol{L}^\mathrm{b} + \frac{d \boldsymbol{L}^\mathrm{b}}{dt} - \frac{1}{2} ({\boldsymbol{\omega}^\mathrm{b}}^* {\boldsymbol{L}^\mathrm{b}}^*)^* \\ &= \frac{1}{2} \boldsymbol{\omega}^\mathrm{b} \boldsymbol{L}^\mathrm{b} + \frac{d \boldsymbol{L}^\mathrm{b}}{dt} - \frac{1}{2} (\boldsymbol{\omega}^\mathrm{b} \boldsymbol{L}^\mathrm{b})^* \\ &= \frac{d \boldsymbol{L}^\mathrm{b}}{dt} + \mathsf{V} (\boldsymbol{\omega}^\mathrm{b} \boldsymbol{L}^\mathrm{b}) \qquad (\because \tilde{a} - \tilde{a}^* = 2\mathsf{V} \tilde{a}) \\ &= \frac{d \boldsymbol{L}^\mathrm{b}}{dt} + \boldsymbol{\omega}^\mathrm{b} \times \boldsymbol{L}^\mathrm{b} \end{align*}$$

となり、式(3)を導出できました。ただし、$\mathsf{V}$は四元数のベクトル部を取り出す演算子です。

簡潔ですね。四元数すごい。因みに、四元数の微分に関しては文献[2]に詳しく書かれています。