プログラミング
古典力学Python

Pythonでハーポルホード曲線を実装する

プログラミング
2026.01.28

概要

本サイトのトップページやロゴアイコンに使用している図形は、物理学の専門用語でハーポルホード（herpolhode）と呼ばれる曲線を正しく描いたものです。
Pythonを用いてハーポルホード曲線を描画するプログラムを自作したので、その実装方法を、完成に至る苦労を交えながら解説します。

実装においては、arXivで見つけた以下の論文で説明されている理論式に従いました。
[2103.16056] Analytical solution of the Euler-Poinsot problem
ただ、論文に記載の数式をそのままプログラムに直すだけでは正しい解は得られないので、注意が必要です。
途中にはJacobiの楕円関数も登場しますので、楕円関数を用いるプログラムの作成の参考にもなると思います。

本記事は実装のみに焦点を当てており、理論の導出は扱いません。
また、以下の記事にてハーポルホードの語源を徹底調査しています。ぜひご覧ください：

【物理学】ハーポルホードの語源に迫る
概要本サイトのサイト名および管理者ペンネームに使われている「ハーポルホード（herpolhode）」は、物理学の専門用語です。しかしながら、身近な場面で応用されるわけでもなく、恩恵が限定的であるためか、物理学界隈でもその存在はあまり知られて...
www.t-herpolhode.com
2026.01.31

ハーポルホード曲線とは

理論は扱わないとはいえ、日常生活ではまず聞かない単語でしょうから、「ハーポルホード曲線」とはそもそも何なのかを簡単に述べることにします。

ハーポルホード高橋
ハーポルホード高橋

逆に日常生活で聞くというかたがいれば、ご一報ください

3次元空間内に物体があるとします。その物体が

  • 有限の質量を持ち、
  • 有限の大きさを持ち、
  • 決して変形しない

とき、剛体と呼びます。

剛体の運動は、重心運動と回転運動に分解することができます。
重心運動は、剛体の質量を $M$、 重心の速度を $\subb{\bm{v}}{G}$、外力を $\bm{F}$ として$$ M\subb{\dot{\bm{v}}}{G} = \bm{F} $$のように見慣れた運動方程式で表現されます。
一方で回転運動は、剛体の重心回りの慣性モーメントを $I$ 、角速度ベクトルを $\bm{\omega}$ 、外力のトルクを $\bm{N}$ として$$ I\dot{\bm{\omega}} = \bm{N} $$という、よく似た運動方程式で記述することができます。

以下では、重力を含む外力が一切無い環境を考えて、重心を固定点とした回転運動のみに注目します。これを剛体の自由回転と言います。
剛体の重心は常に原点に固定します。($\subb{\bm{v}}{G}=\subb{\bm{x}}{G}=0$）

慣性モーメント $I$ は実対称な3×3行列であり、ざっくり言うと数値が大きいほど回転しにくいと解釈できます。剛体の形状が下図のような楕円体のとき、凸の向きに合わせてx,y,z軸をとると、$I$ は対角行列となります。$$ I = \pmatrix{I_x & 0 & 0 \\ 0 & I_y & 0 \\ 0 & 0 & I_z} $$

実はどんな複雑な形の剛体でも、適切に座標変換することで、楕円体と同じ運動に帰着できることが示されます。このときに対応する仮想的な楕円体を慣性楕円体、$I$ が対角になるように上手く選んだ軸を慣性主軸、対角成分 $I_x,\,I_y,\,I_z$ を主慣性モーメントと呼びます。
そういう訳で一般の剛体に適用できる理論なのですが、イメージし辛い場合は、初めから楕円体が回転していると思って読んで頂いて差支えありません。

角速度ベクトル $\bm{\omega}$ についても補足しておきます。
角速度ベクトルを用いると、どの方向の軸の周りに、どのくらいの速さで回転するのかをまとめて表現できます。ベクトルの始端は原点に固定です。
軸と回転の向きは、右ねじの法則で決まります。例えば $\bm{\omega}={}^t(0\,0\,1)$ であれば、z軸を回転の軸として、x軸正→y軸正→x軸負→y軸負→…という向きに角速度 1 [rad/s] で回転する、という具合です。

実験室系（すなわち静止している人やカメラ）から見て、自由回転する慣性楕円体を観察します。
どのように回転するかというと、回転軸が慣性主軸と一致しているときは、直立した独楽こまのように角速度が一定の安定した回転になります。しかし一般には、回転軸の向きや角速度の大きさが変化する複雑な回転をします。つまり、角速度ベクトル $\bm{\omega}$ は時間変化することになります。

エネルギー保存則および角運動量保存則を用いると、

  • 慣性楕円体は、不変平面と呼ばれる空間内に固定された単一面上を滑らずに転がる
  • 慣性楕円体と不変平面の接点は、慣性楕円体と角速度ベクトル $\bm{\omega}$ の交点である
  • 角速度ベクトル $\bm{\omega}$ の終端は、不変平面と平行な単一面上を動く

ことが示されます。この手法をPoinsotの作図法と呼びます。

これでようやく、ハーポルホードの定義に到達します。

定義：
慣性楕円体と不変平面の接点が、不変平面上に描く軌跡をハーポルホードと呼ぶ

ついでに、慣性楕円体と不変平面の接点が慣性楕円体上に描く軌跡はポルホード（polhode）と呼びます。
ポルホードは必ず閉じた曲線であるのに対して、ハーポルホードは一般には閉じない複雑な曲線となります。

なお、剛体の自由回転には「ハーポルホード錐」という概念も登場しますが、「ハーポルホード曲線」とは異なります。

実装方針

プログラムを実装する前に、前提条件や記法を明らかにしておきます。

・実験室系と剛体系

出力したいハーポルホード曲線は、静止した人から見るのと同じ実験室系（慣性系）で観測されるものです。
しかし理論面からすると、剛体上に固定されて剛体とともに回転する剛体系（非慣性系）を考えるほうが取り扱いやすい事情があります。

そのため、剛体系で時間発展を計算した結果を実験室系に移して最終出力を得る、という構成を取ります。
重心運動は考えず回転運動のみに着目するため、原点 O は両者共通の固定点です。

参考論文では、実験室系の座標は大文字の「XYZ」、剛体系の座標は小文字の「xyz」と区別されています。
本記事のプログラム内では、実験室系は添字「lab」、剛体系は添字「rigid」を付けて表すことにします。
1軸, 2軸, 3軸 は、 X軸, Y軸, Z軸 または x軸, y軸, z軸 に適宜読み替えてください。

実験室系と剛体系の間は、Euler角 $(\psi, \theta, \phi)$ で互いに移り変わることができます。これらの角度変数は時間 $t$ に依存して変化します。
実験室系の角運動量ベクトル $\subb{\bm{L}}{lab}$ は、角運動量保存則より定ベクトルとなります。$\subb{\bm{L}}{lab}$ がZ軸方向となるように、実験室系の座標の向きを選ぶものとします。

Eular angle
$\phi,\theta,\psi$ の順に回転させる

・パラメータ

人が与えるパラメータとして、剛体系の主慣性モーメント $(I_x, I_y, I_z)$ と角速度の初期条件 $(\omega_{0x}, \omega_{0y}, \omega_{0z})$ の計6個を選びます。

論文に揃えて、一般性を失うことなく $I_x>I_y>I_z>0$ に限定します。（例外処理が面倒になるだけなので、分母が0になるような特殊値には対応しません）

・出力の定義

本記事では、実験室系の角速度ベクトルが描く軌跡を、最終的な出力のハーポルホード曲線とします。

定義に厳密に従うなら、不変平面上の軌跡がハーポルホード曲線であり、角速度ベクトルの軌跡とは一般には一致しません。
しかし、角速度ベクトルを含む平面は不変平面と平行なので、角速度ベクトルが描く軌跡は定義通りのハーポルホード曲線を拡大縮小したものになります。

今回はハーポルホード曲線の図形を得ることが動機であり、倍率の情報は重要でないため、不変平面との交点を求めるステップは省略させて頂きます。

・プログラムの構成

以上の仮定の下で、プログラムに計算させるべき内容は次の通りです。

  1. 入力から決まる諸定数を算出する
  2. 剛体系における角速度ベクトルの時間発展を求める
  3. Euler角の時間発展を求める
  4. 剛体系の角速度ベクトルをEuler角で実験室系に変換する
  5. 実験室系の角速度ベクトルを出力する

実装内容の解説

前置きが長くなりましたが、いよいよPythonの実装に入ります。

場合分けと定数計算

前述の通り、始めにパラメータ $(I_x, I_y, I_z)$ および $(\omega_{0x}, \omega_{0y}, \omega_{0z})$ を受け取ります。

エネルギー保存則と角運動量保存則から、2つの運動の定数 $$ F = \frac{1}{2}\left( I_x\omega_{x0}^2+I_y\omega_{y0}^2+I_z\omega_{z0}^2 \right) $$ $$ G^2 = I_x^2\omega_{x0}^2+I_y^2\omega_{y0}^2+I_z^2\omega_{z0}^2 $$ が導かれます。
$F$ は回転エネルギー、$G^2$ は角運動量ベクトルの大きさ（ノルム）の2乗を意味します。
論文では $F$ ではなく $T$ と表記されてありますが、時刻を表す $t$ と紛らわしいので $F$ を用いることにします。

パラメータ $I_x, I_y, I_z$ と $F, G^2$ との大小により、以下の2つに場合分けされます。

  • $2TI_x>2TI_y>G^2>2TI_z$ のとき g2case=1
  • $2TI_x>G^2>2TI_y>2TI_z$ のとき g2case=2

ケースごとに用いる数式と解が異なってくるので、ケースを判定後、定数 $P^2, Q^2, R^2, n, k$ を計算します。

g2case=1 のとき
$$ \begin{align*} P^2 &= \frac{G^2 -2FI_z}{I_x^2-I_x I_z} \\ Q^2 &= \frac{G^2 -2FI_z}{I_y^2-I_y I_z} \\ R^2 &= \frac{G^2 -2FI_x}{I_z^2-I_x I_z} \end{align*} $$ $$ n = \sqrt{\frac{(I_y-I_z)(2FI_x-G^2)}{I_x I_y I_z}} $$ $$ k = \sqrt{\frac{I_x-I_y}{I_y-I_z}\frac{G^2-2FI_z}{2FI_x-G^2}} $$

g2case=2 のとき
$$ \begin{align*} P^2 &= \frac{2FI_z-G^2}{I_x I_z-I_x^2} \\ Q^2 &= \frac{2FI_x-G^2}{I_x I_y-I_y^2} \\ R^2 &= \frac{2FI_x-G^2}{I_x I_z-I_z^2} \end{align*} $$ $$ n = \sqrt{\frac{(I_x-I_y)(G^2-2FI_z)}{I_x I_y I_z}} $$ $$ k = \sqrt{\frac{I_y-I_z}{I_x-I_y}\frac{2FI_x-G^2}{G^2-2FI_z}} $$

いずれのケースにおいても、$0<k<1$ が定義より保証されています。
この性質は、後で楕円関数の定義域を考える上で重要になります。

問1

g2case=1 g2case=2 ともに、$0<k<1$

g2case=1 のとき $$ \begin{align*} &(I_y-I_z)(2FI_x-G^2) -(I_x-I_y)(G^2-2FI_z) \\ &= 2F(I_xI_y-\cancel{I_xI_z}+\cancel{I_xI_z}-I_yI_z) + G^2(-\cancel{I_y}+I_z-I_x+\cancel{I_y}) \\ &= 2F(I_x-I_z)I_y -G^2(I_x-I_z) \\ &= (2FI_y-G^2)(I_x-I_z) >0 \end{align*} $$ よって $$ k^2 = \frac{I_x-I_y}{I_y-I_z}\frac{G^2-2FI_z}{2FI_x-G^2} < 1 $$ g2case=2 のとき $$ \begin{align*} &(I_x-I_y)(G^2-2FI_z) -(I_y-I_z)(2FI_x-G^2) \\ &= 2F(-\cancel{I_xI_z}+I_yI_z-I_xI_y+\cancel{I_xI_z}) + G^2(I_x-\cancel{I_y}+\cancel{I_y}-I_z) \\ &= -2F(I_x-I_z)I_y +G^2(I_x-I_z) \\ &= (G^2-2FI_y)(I_x-I_z) >0 \end{align*} $$ よって $$ k^2 = \frac{I_y-I_z}{I_x-I_y}\frac{2FI_x-G^2}{G^2-2FI_z} < 1 $$

剛体系の時間発展

論文によると、運動方程式の解はJacobiの楕円関数 $\mathrm{sn, cn, dn}$ を用いて与えられます。

g2case=1 のとき
$$ \begin{align} \omega_x(t) &= \mathrm{sgn}(\omega_{z0})\sqrt{P^2}\,\mathrm{cn}(nt+\tau, k) \label{eq:1rigid_x} \\ \omega_y(t) &= -\sqrt{Q^2}\,\mathrm{sn}(nt+\tau, k) \label{eq:1rigid_y} \\ \omega_z(t) &= \mathrm{sgn}(\omega_{z0})\sqrt{R^2}\,\mathrm{dn}(nt+\tau, k) \label{eq:1rigid_z} \end{align} $$

g2case=2 のとき
$$ \begin{align} \omega_x(t) &= \mathrm{sgn}(\omega_{x0})\sqrt{P^2}\,\mathrm{dn}(nt+\tau, k) \label{eq:2rigid_x} \\ \omega_y(t) &= -\sqrt{Q^2}\,\mathrm{sn}(nt+\tau, k) \label{eq:2rigid_y} \\ \omega_z(t) &= \mathrm{sgn}(\omega_{x0})\sqrt{R^2}\,\mathrm{cn}(nt+\tau, k) \label{eq:2rigid_z} \end{align} $$

変数 $\tau$ はいずれのケースも共通で、以下の積分で定義されます。
$\omega_{x0}\omega_{z0}\geq0$ のとき $$ \begin{equation} \tau = \int_0^{\tau_0}\!\frac{\mathrm{d}u}{\sqrt{(1-u^2)(1-k^2u^2)}} \;,\;\; \tau_0 = -\frac{\omega_{y0}}{\sqrt{Q^2}} \label{eq:tau_def1} \end{equation} $$ $\omega_{x0}\omega_{z0}< 0$ のとき $$ \begin{equation} \tau = \frac{K(k)}{2} -\int_0^{\tau_0}\!\frac{\mathrm{d}u}{\sqrt{(1-u^2)(1-k^2u^2)}} \;, \label{eq:tau_def2} \end{equation} $$ $$ \begin{equation} \tau_0 = -\frac{\omega_{y0}}{\sqrt{Q^2}} \;,\;\; K(k) = 4\int_0^{\frac{\pi}{2}}\!\frac{\mathrm{d}u}{\sqrt{1-k^2\sin^2u}} \label{eq:K_def} \end{equation} $$

Jacobiの楕円関数

式\eqref{eq:1rigid_x}～\eqref{eq:2rigid_z}に用いられているJacobiの楕円関数ですが、定義が流派によって異なるようなので、一度定義と性質を整理してから使うことにします。
ここではJacobiによる積分を用いた定義式を紹介します。

定義：
$$ \begin{equation} x = \int_0^\varphi\!\frac{\mathrm{d}\theta}{\sqrt{1-m\sin^2\theta}} = \int_0^{\sin\varphi}\!\frac{\mathrm{d}u}{\sqrt{(1-u^2)(1-mu^2)}} \label{eq:Jacobi_xdef} \end{equation} $$ に対して $$ \mathrm{sn}(x;m) = \sin\varphi $$ $$ \mathrm{cn}(x;m) = \cos\varphi $$ $$ \mathrm{dn}(x;m) = \sqrt{1-m\sin^2\varphi} $$ を満たす関数 $\mathrm{sn, cn, dn}$ を、Jacobiの楕円関数とする

$m\;(0\leq m\leq1)$ は母数と呼ばれます。
式\eqref{eq:Jacobi_xdef}の2つ目の等号は、$u=\sin\theta$ に置換積分することで示せます。

この定義を採用したとき、値域と周期は以下のようになります。

値域周期
$-1\leq\mathrm{sn}(x;m)\leq1$$4K(m)$
$-1\leq\mathrm{cn}(x;m)\leq1$$4K(m)$
$\sqrt{1-m}\leq\mathrm{dn}(x;m)\leq1$$2K(m)$

ただし $$ \begin{equation} K(m) = \int_0^{\frac{\pi}{2}}\!\frac{\mathrm{d}\theta}{\sqrt{1-m\sin^2\theta}} \label{eq:ellipconp} \end{equation} $$

特に $m=0$ のとき $\mathrm{sn}(x;0)=\sin(x),\,\mathrm{cn}(x;0)=\cos(x),\,$$\mathrm{dn}(x;0)=1,\,K(0)=\tfrac{\pi}{2}$ となることから、楕円関数は三角関数を拡張したものであると理解できます。

勘のいい読者は既にお気付きの通り、ハーポルホード曲線の論文は上記とは異なる定義で表記されています。論文の $k^2$ を 式\eqref{eq:Jacobi_xdef}の母数 $m$ と対応させる必要があります。
根拠としては、式\eqref{eq:tau_def1}～\eqref{eq:K_def}の積分を式\eqref{eq:Jacobi_xdef}, \eqref{eq:ellipconp}と整合させること、そして後述する保存量の計算結果を時間発展で不変にすることが理由に挙げられます。

一方、PythonでJacobiの楕円関数は、scipyのspecial.ellipjで実装することができます。
scipy.special.ellipj — SciPy v1.17.0 Manual

Python
from scipy.special import ellipj
sn, cn, dn, ph = ellipj(x, m)

今回は使用しませんが、変数 ph には式\eqref{eq:Jacobi_xdef}の $\varphi$ が格納されます。
$m=0.4$ として得られた数値をプロットすると、以下のようになります。

visualization of ellipj

$\mathrm{dn}$ の下限値が $\sqrt{1-m}\approx0.775$ であることから、scipyの母数 m は式\eqref{eq:Jacobi_xdef}の $m$ に対応することが分かります。

第1種不完全楕円積分

次は、式\eqref{eq:tau_def1},\eqref{eq:tau_def2}に共通して含まれている積分です。$$ \int_0^{\tau_0}\!\frac{\mathrm{d}u}{\sqrt{(1-u^2)(1-k^2u^2)}} $$

解析的に積分実行することはできず、第1種不完全楕円積分と名前が付いています。
式\eqref{eq:Jacobi_xdef}にも全く同じ形が登場しており、Jacobiの楕円関数は第1種不完全楕円積分の逆関数という関係性になっています。

Pythonで第1種不完全楕円積分は、scipyのspecial.ellipkincで実装することができます。
scipy.special.ellipkinc — SciPy v1.17.0 Manual

Python
from scipy.special import ellipkinc
y = ellipkinc(phi, m)

scipy公式ドキュメントには $$ \begin{equation} K(\phi,m) = \int_0^\phi\left[1-m\sin(t)^2\right]^{-1/2}dt \label{eq:ellipkinc_scipy} \end{equation} $$ という定義式が書かれてあり、やはり母数は上記の $m$ のほうで定義されています。

論文の表記をscipyの定義に合わせるため、式\eqref{eq:Jacobi_xdef}の式変形を適用し、次のように書き直します。 $$ \begin{equation} \int_0^{\tau_0}\!\frac{\mathrm{d}u}{\sqrt{(1-u^2)(1-k^2u^2)}} = \int_0^{\mathrm{Arcsin}(\tau_0)}\!\frac{\mathrm{d}\theta}{\sqrt{1-k^2\sin^2\theta}} \label{eq:ellipkinc_k2} \end{equation} $$
よって、scipyの引数 phi, m には、論文の $\mathrm{Arcsin}(\tau_0),\,k^2$ を用いればよいことが分かります。

引数の定義域の確認をします。scipyのellipkinc関数は $0\leq m\leq1$ においてのみ定義されており、それ以外の範囲を入力するとnanになります。
今回の場合、$k$ の定義から $0<k<1$ が保証されていることを前節で示したのでした。そのため、問題なく $m=k^2$ を用いることができます。

第1種完全楕円積分

式\eqref{eq:K_def}で係数 $4$ を無視した積分 $$ \int_0^{\frac{\pi}{2}}\!\frac{\mathrm{d}u}{\sqrt{1-k^2\sin^2u}} $$ は、式\eqref{eq:Jacobi_xdef}で $\varphi=\tfrac{\pi}{2}$ を代入した特別な場合に相当します。こちらは第1種完全楕円積分と呼ばれています。
係数の付け方が異なりますが、楕円関数の周期 $K$ として、式\eqref{eq:ellipconp}にも同じものが登場しています。

実装においては、第1種不完全楕円積分の関数で phi を指定しても良いのですが、scipyにはspecial.ellipkという専用の関数が用意されているので、そちらを使うことにします。
scipy公式ドキュメントには、ellipkincの定義式\eqref{eq:ellipkinc_scipy}で $\phi=\tfrac{\pi}{2}$ としたものが書かれてあります。
scipy.special.ellipk — SciPy v1.17.0 Manual

Python
from scipy.special import ellipk
y = ellipk(m)

scipyのellipk関数を可視化してみると、下図のようになります。

visualization of ellipk

第1種完全楕円積分は、特殊値として $$ K(m=0) = \frac{\pi}{2} \approx 1.57080$$ $$ K(m=\tfrac{1}{2}) = \frac{\Gamma(\tfrac{1}{4})^2}{4\sqrt{\pi}} \approx 1.85407\cdots $$ が知られています。
上図のグラフとも矛盾していないことが確かめられます。

Euler角の時間発展

実験室系と剛体系を結びつけるEuler角は、時刻 $t$ の関数として、以下のように与えられます。 $$ \theta(t) = \mathrm{Arccos}\left(\frac{I_z\omega_z(t)}{\sqrt{G^2}}\right) $$ $$ \phi(t) = \mathrm{Arctan2}(I_x\omega_x(t),\,I_y\omega_y(t)) $$ $$ \dv{\psi(t)}{t} = \sqrt{G^2}\frac{I_x\omega_x^2(t)+I_y\omega_y^2(t)}{I_x^2\omega_x^2(t)+I_y^2\omega_y^2(t)} \;,\; \psi(0) = 0 $$

論文には楕円関数を用いて $\psi(t)$ を表した式も書かれてますが、今回は時間微分の表式を数値的に $t$ 積分して求めることにします。この方法なら、ケース別の場合分けも不要になります。 $$ \psi(t+\Delta t) \simeq \psi(t)+\dv{\psi(t)}{t}\Delta t $$

逆三角関数

逆三角関数は、上記のEuler角のほかに $\tau_0$ に対する式\eqref{eq:ellipkinc_k2}にも登場します。

ハーポルホード高橋
ハーポルホード高橋

私は普段、主値を限定しないものを arc○○、主値を限定して一価関数にしたものを Arc○○ と書くことにしています。
マイナス1乗の書き方は、逆数と紛らわしいので好きではありません。

Pythonで逆三角関数は、numpyで実装することができます。
numpy.arcsin — NumPy v2.4 Manual
numpy.arccos — NumPy v2.4 Manual
numpy.arctan — NumPy v2.4 Manual
numpy.arctan2 — NumPy v2.4 Manual

Python
import numpy as np
f1 = np.arcsin(x)
f2 = np.arccos(x)
f3 = np.arctan(x)
f4 = np.arctan2(x, y)

値域（戻り値）の範囲をまとめると、以下のようになります。すべてラジアンです。

関数名値域
np.arcsin(x)$[-\pi/2,\,\pi/2]$
np.arccos(x)$[0,\,\pi]$
np.arctan(x)$[-\pi/2,\,\pi/2]$
np.arctan2(x,y)$[-\pi,\,\pi]$

Arctan2は2つの引数を持ち、原点から $(x,y)$ を見たときの偏角を返します。
それに対してArctanはtanの逆関数であり、傾き $y/x$ に対応する角度を返します。第1・第3象限および第2・第4象限を区別しません。
今回の状況設定では、$\omega_x(t), \omega_y(t)$ の符号の組み合わせを考慮する必要があるため、ArctanではなくArctan2を使用します。

補足として、numpyにはnp.asin, np.atan2 等の関数も別名で用意されており、同じ結果を得ることができます。
また、標準ライブラリのmath.asin, math.atan2 等も同じ定義ですが、pd.Series型を引数に入れることができるのはnumpyのみです。

剛体系を実験室系に変換

剛体系で観測されたベクトルと実験室系で観測されたベクトルは、時間に依存する変換行列 $R(t)$ とその逆行列 $R^{-1}(t)$ を用いて互いに変換されます。
$$ \subb{\bm{v}}{rigid}(t) = R(t) \subb{\bm{v}}{lab}(t) \;\Longleftrightarrow\; \subb{\bm{v}}{lab}(t) = R^{-1}(t) \subb{\bm{v}}{rigid}(t) $$ 具体的には、角速度ベクトル $\subb{\bm{\omega}}{rigid}$ と角運動量ベクトル $\subb{\bm{L}}{rigid} = I\subb{\bm{\omega}}{rigid}$ を、剛体系から実験室系へ変換することになります。

オイラー角 $(\psi(t), \theta(t), \phi(t))$ の定義より、$R(t)$ は以下で与えられます。 $$ \begin{align*} R(t) &= R_\psi R_\theta R_\phi \\ &= \begin{pmatrix} \cos\psi & \sin\psi & 0 \\ -\sin\psi & \cos\psi & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos\theta & \sin\theta \\ 0 & -\sin\theta & \cos\theta \end{pmatrix} \begin{pmatrix} \cos\phi & \sin\phi & 0 \\ -\sin\phi & \cos\phi & 0 \\ 0 & 0 & 1 \end{pmatrix} \\ &= \begin{pmatrix} \cos\psi\cos\phi-\sin\psi\cos\theta\sin\phi & \sin\psi\cos\phi+\cos\psi\cos\theta\sin\phi & \sin\theta\sin\phi \\ -\cos\psi\sin\phi-\sin\psi\cos\theta\cos\phi & -\sin\psi\sin\phi+\cos\psi\cos\theta\cos\phi & \sin\theta\cos\phi \\ \sin\psi\sin\theta & -\cos\psi\sin\theta & \cos\theta \end{pmatrix} \end{align*} $$

問2

$R(t)$ は直交行列、すなわち $R^{-1}(t)={}^t\!R(t)$

$R_\psi,\,R_\theta,\,R_\phi$ はそれぞれ直交行列であるから、$R_\bullet^{-1} = {}^t\!R_\bullet$ を満たす。 $$ \begin{align*} R^{-1} &= (R_\psi R_\theta R_\phi)^{-1} \\ &= R_\phi^{-1}\, R_\theta^{-1}\, R_\psi^{-1} \\ &= {}^t\!R_\phi \,{}^t\!R_\theta \,{}^t\!R_\psi \\ &= {}^t(R_\psi R_\theta R_\phi) \\ &= {}^t\!R \end{align*} $$

この性質から、逆変換 $R^{-1}(t)$ の表式は以下のようになります。 $$ R^{-1}(t) = \begin{pmatrix} \cos\psi\cos\phi-\sin\psi\cos\theta\sin\phi & -\cos\psi\sin\phi-\sin\psi\cos\theta\cos\phi & \sin\psi\sin\theta \\ \sin\psi\cos\phi+\cos\psi\cos\theta\sin\phi & -\sin\psi\sin\phi+\cos\psi\cos\theta\cos\phi & -\cos\psi\sin\theta \\ \sin\theta\sin\phi & \sin\theta\cos\phi & \cos\theta \end{pmatrix} $$

グラフ出力

以上の計算を、pandasを用いて実行します。
時間を微小な時間幅 dt のステップに区切って、行方向に時刻（＝ステップ数）を並べ、列方向に各物理量を追加していきます。計算過程は、以下のようなpd.DataFrame型変数にまとめることにします。

example of dataframe

実験室系の角速度ベクトル $\subb{\bm{\omega}}{lab}$ のXY座標が、ハーポルホード曲線となります。
w1_lab, w2_lab の2列を取り出し、matplotlibで2次元平面にプロットします。
修飾前のコードを以下に示します。

Python
import pandas as pd
import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.plot(df['w1_lab'], df['w2_lab'])

ax.set_xlabel('x_lab')
ax.set_ylabel('y_lab')
ax.set_aspect('equal', adjustable='box')  # 軸のアスペクト比を揃える
ax.grid(False)

plt.show()

また、検算時には $\subb{\bm{\omega}}{lab}$ が本当に平面上を動いているか・そもそも物理的に有効な解になっているのかの確認をすることは重要です。
理論からの帰結で、以下の変数は保存量、つまり時間に依らない定数となるべきです。DataFrameの該当する列を見て、計算結果が妥当かのチェックに使うことができます。

  • 回転エネルギー $\;F=\tfrac{1}{2}\left( I_x\omega_{x}^2+I_y\omega_{y}^2+I_z\omega_{z}^2 \right)$
  • 角運動量の大きさ $\;G^2 = I_x^2\omega_{x}^2+I_y^2\omega_{y}^2+I_z^2\omega_{z}^2$
  • 実験室系の角運動量ベクトル $\;\subb{\bm{L}}{lab} = R^{-1}I\subb{\bm{\omega}}{rigid}$
  • 実験室系の角速度ベクトルのZ成分 $\;\subb{\omega}{Z} = \left( R^{-1}\subb{\bm{\omega}}{rigid} \right)_Z$

完成したコード

最終的なコードを公開します。
論文に記載がある計算例については、結果が一致することを確認できました。

Python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.special import ellipj
from scipy.special import ellipk
from scipy.special import ellipkinc

# df作成
def calc_w(I1, I2, I3, W10, W20, W30, dt=0.01, t_max=50):

    # 定数を計算
    f = (I1*W10**2 + I2*W20**2 + I3*W30**2) * 0.5
    g2 = I1**2*W10**2 + I2**2*W20**2 + I3**2*W30**2
    g = np.sqrt(g2)
    print('f =', f)
    print('g2 =', g2)
    print('2fI1, 2fI2, 2fI3 =', 2*f*I1, ',', 2*f*I2, ',', 2*f*I3)

    if I1>I2 and I2>I3 and I3>0:
        if 2*f*I2>g2 and g2>2*f*I3:
            g2_case = 1
            print('g2_case =', g2_case)
        elif 2*f*I1>g2 and g2>2*f*I2:
            g2_case = 2
            print('g2_case =', g2_case)
        else:
            print('Error: g2_case undefined')
            return None
    else:
        print('Error: invarid I1, I2, I3')
        return None
    

    # 剛体系の角速度ベクトルの時間発展を計算

    # =========================
    # g2_case==1: 2fI1 > 2fI2 > g2 > 2fI3
    # =========================
    if g2_case==1:
        p2 = (g2 - 2*f*I3)/(I1*I1 - I1*I3)
        q2 = (g2 - 2*f*I3)/(I2*I2 - I2*I3)
        r2 = (g2 - 2*f*I1)/(I3*I3 - I1*I3)
        n = np.sqrt((I2 - I3)*(2*f*I1 - g2)/(I1*I2*I3))
        k = np.sqrt(((I1 - I2)/(I2 - I3)) * ((g2 - 2*f*I3)/(2*f*I1 - g2)))
        print('k =', k)
        tau0 = -W20 / np.sqrt(q2)
        
        if W10*W30>=0:
            tau = ellipkinc(np.arcsin(tau0), k**2)
        else:
            tau = 2*ellipk(k**2) - ellipkinc(np.arcsin(tau0), k**2)
        
        t = np.arange(0.0, t_max, dt)
        s = n * t + tau

        sn, cn, dn, am = ellipj(s, k**2)

        w1 = np.sign(W30) * np.sqrt(p2) * cn
        w2 = -np.sqrt(q2) * sn
        w3 = np.sign(W30) * np.sqrt(r2) * dn


    # =========================
    # g2_case==2: 2fI1 > g2 > 2fI2 > 2fI3
    # =========================
    elif g2_case==2:
        p2 = (2*f*I3 - g2)/(I1*I3 - I1*I1)
        q2 = (2*f*I1 - g2)/(I1*I2 - I2*I2)
        r2 = (2*f*I1 - g2)/(I1*I3 - I3*I3)
        n = np.sqrt((I1 - I2)*(g2 - 2*f*I3)/(I1*I2*I3))
        k = np.sqrt(((I2 - I3)/(I1 - I2)) * ((g2 - 2*f*I1)/(2*f*I3 - g2)))
        print('k =', k)
        tau0 = -W20 / np.sqrt(q2)
        
        if W10*W30>=0:
            tau = ellipkinc(np.arcsin(tau0), k**2)
        else:
            tau = 2*ellipk(k**2) - ellipkinc(np.arcsin(tau0), k**2)
        
        t = np.arange(0.0, t_max, dt)
        s = n * t + tau

        sn, cn, dn, am = ellipj(s, k**2)
        
        w1 = np.sign(W10) * np.sqrt(p2) * dn
        w2 = -np.sqrt(q2) * sn
        w3 = np.sign(W10) * np.sqrt(r2) * cn


    df = pd.DataFrame({'w1_rigid': w1, 'w2_rigid': w2, 'w3_rigid': w3})

    df['I1'] = I1
    df['I2'] = I2
    df['I3'] = I3

    df['f'] = (I1*df['w1_rigid']**2 + I2*df['w2_rigid']**2 + I3*df['w3_rigid']**2) * 0.5
    df['g2'] = I1**2*df['w1_rigid']**2 + I2**2*df['w2_rigid']**2 + I3**2*df['w3_rigid']**2

    df['l1_rigid'] = I1*df['w1_rigid']
    df['l2_rigid'] = I2*df['w2_rigid']
    df['l3_rigid'] = I3*df['w3_rigid']


    # 実験室系に対する剛体系の立体角(rad)の時間発展を計算

    df['theta'] = np.arccos(df['l3_rigid']/g)
    df['phi'] = np.arctan2(df['l1_rigid'], df['l2_rigid'])
    df['dt_psi'] = g * (df['l1_rigid']*df['w1_rigid'] + df['l2_rigid']*df['w2_rigid']) / (df['l1_rigid']*df['l1_rigid'] + df['l2_rigid']*df['l2_rigid'])
    df['psi'] = ((df['dt_psi'] * dt).cumsum()) % (2*np.pi)     # psiの初期値は0


    # 回転行列で実験室系に変換

    df['R11'] = np.cos(df['psi'])*np.cos(df['phi']) - np.sin(df['psi'])*np.sin(df['phi'])*np.cos(df['theta'])
    df['R12'] = -np.cos(df['psi'])*np.sin(df['phi']) - np.sin(df['psi'])*np.cos(df['phi'])*np.cos(df['theta'])
    df['R13'] = np.sin(df['psi'])*np.sin(df['theta'])
    df['R21'] = np.sin(df['psi'])*np.cos(df['phi']) + np.cos(df['psi'])*np.sin(df['phi'])*np.cos(df['theta'])
    df['R22'] = -np.sin(df['psi'])*np.sin(df['phi']) + np.cos(df['psi'])*np.cos(df['phi'])*np.cos(df['theta'])
    df['R23'] = -np.cos(df['psi'])*np.sin(df['theta'])
    df['R31'] = np.sin(df['phi'])*np.sin(df['theta'])
    df['R32'] = np.cos(df['phi'])*np.sin(df['theta'])
    df['R33'] = np.cos(df['theta'])

    df['w1_lab'] = df['R11']*df['w1_rigid'] + df['R12']*df['w2_rigid'] + df['R13']*df['w3_rigid']
    df['w2_lab'] = df['R21']*df['w1_rigid'] + df['R22']*df['w2_rigid'] + df['R23']*df['w3_rigid']
    df['w3_lab'] = df['R31']*df['w1_rigid'] + df['R32']*df['w2_rigid'] + df['R33']*df['w3_rigid']

    df['l1_lab'] = df['R11']*df['l1_rigid'] + df['R12']*df['l2_rigid'] + df['R13']*df['l3_rigid']
    df['l2_lab'] = df['R21']*df['l1_rigid'] + df['R22']*df['l2_rigid'] + df['R23']*df['l3_rigid']
    df['l3_lab'] = df['R31']*df['l1_rigid'] + df['R32']*df['l2_rigid'] + df['R33']*df['l3_rigid']


    # エラーチェック
    
    if (abs(df['f'].iloc[0] - df['f'].iloc[-1])>1e-3):
        print('Error: f is not conserved')
    if (abs(df['g2'].iloc[0] - df['g2'].iloc[-1])>1e-3):
        print('Error: g2 is not conserved')
    if (abs(df['l1_lab'].iloc[0] - df['l1_lab'].iloc[-1])>1e-3) or (abs(df['l2_lab'].iloc[0] - df['l2_lab'].iloc[-1])>1e-3) or (abs(df['l3_lab'].iloc[0] - df['l3_lab'].iloc[-1])>1e-3):
        print('Error: l_lab is not conserved')
    if (abs(df['w3_lab'].iloc[0] - df['w3_lab'].iloc[-1])>1e-3):
        print('Error: w3_lab is not conserved')

    return df

# 2dグラフ作成
def draw_w_lab2(df, save=''):
    fig, ax = plt.subplots()
    ax.plot(df['w1_lab'], df['w2_lab'], c='blue', lw=0.5)

    str_I = f'({df['I1'].iloc[0]:.2f}, {df['I2'].iloc[0]:.2f}, {df['I3'].iloc[0]:.2f})'
    str_w0_rigid = f'({df['w1_rigid'].iloc[0]:.2f}, {df['w2_rigid'].iloc[0]:.2f}, {df['w3_rigid'].iloc[0]:.2f})'
    ax.set_title('I='+str_I+', w0_rigid='+str_w0_rigid, fontsize=10)

    ax.set_xlabel('x_lab')
    ax.set_ylabel('y_lab')
    ax.set_aspect('equal', adjustable='box')  # 軸のアスペクト比を揃える
    ax.grid(False)

    if save!='':
        fig.savefig(save+'.png', dpi=200, bbox_inches='tight')
    plt.show()


# main
df = calc_w(3.0, 2.0, 1.0, 2.0, 3.0, 4.0)
df.to_csv('output_df.csv')
draw_w_lab2(df, save='output_fig')
Expand

ハーポルホード曲線の例

実行時のパラメータは、画像上部に記載の通りです。

時間発展の最初のほう

g2case=1 の例　不規則な動きが見て取れます。

g2case=2 の例

長期間の時間発展

小さなずれが重なって、軌跡の濃淡が見えてきます。
保存量は高い精度で一致し続けているため、数値積分による誤差の影響は小さいと考えられます。

さまざまな模様

パラメータを色々変えてみると、美しい図形が得られます。

閉曲線になる場合

パラメータを上手に選ぶと、閉曲線になります。
論文には、楕円関数の周期を使った閉曲線が実現するための条件式も掲載されています。

下図は本サイトのロゴアイコンの元ネタです。

参考文献

タイトルとURLをコピーしました