TL,DR
リー群の自由度をもつ系のラグランジュ力学を考えるとき、リー群の座標を張らなくてもいいらしい。
本編
剛体にSO(3)の回転自由度を与えた系を考える。
剛体の慣性主軸座標を考える。これをローカル座標とする。これと別にワールド座標を考える。ローカル座標系基底ベクトルをワールド座標系の基底ベクトルで展開したときの成分を縦に並べて配置した行列をRとする。したがってこれは能動変換であり、逆にあるベクトルの成分をワールド座標系のそれからローカル座標へ座標変換するときはR^{-1}が必要になる。
今剛体には、R(t)\in SO(3)という配位自由度があることになる。この時間微分を考えると、
\frac{d}{dt} R(t) = R(t)\Omega(t)
とできるようなリー代数\Omega \in \mathfrak{so}(3)を取ることができる。つまり、接バンドルの実体 = SO(3)の微分係数は、\mathfrak{so}(3)の(その場所への)左移動だと思うことができる。
\Omegaと\omegaは、前者からリー代数基底で成分を取り出した物を後者だとする。つまり\Omega_{ij}=-\epsilon_{ijk}\omega_kとする。
さて、回転群自由度の「速度」がR(t)\Omega(t)であるとき、運動エネルギーはどうなっているだろうか？これは、今すでにR(t)回っている状況から、更にR(t)\Omega(t)継ぎ足すように動いているということだ。言い換えれば、この瞬間のローカル座標系では\Omega(t)だけ回ろうとしている。
例えばローカル座標系で位置ベクトル\mathbf{r}の場所にあるものは、\Omega \mathbf{r}の速度を得ている。したがって系の運動エネルギーとしては
\frac{1}{2}\int dr^3 \rho(\mathbf{r})\mathbf{r}^\top \Omega^\top \Omega \mathbf{r}
になっている。二次モーメントテンソルを
M = \int dr^3 \mathbf{r} \rho(\mathbf{r})\mathbf{r}^\top
とすれば
\frac{1}{2}\mathrm{Tr}\Omega M \Omega^\top
である。これは、慣性モーメントを使った伝統的な
\frac{1}{2}\omega^\top I \omega
とは互換性がある。\Omega_{ij} = -\epsilon_{ijk} \omega_kなので、完全反対称テンソルの計算から
I = \mathrm{Tr}M \mathbb{I} - M
のように二次モーメントと慣性モーメントは対応している。
もし、回転の中心が重心ではなく、重力を無視できないなら、ポテンシャル関数が作れる。初期配置のときのワールド座標系での重心位置ベクトルを\mathbf{r}_0、総重量をm、重力加速度方向を\mathbf{g}としよう。ポテンシャル関数は以下のようにできるが、もっと変なものを考えてもよい。
-m\mathbf{g}\cdot R \mathbf{r}_0
さて、運動方程式を考えよう。ラグランジュ形式だ。ふつうのラグランジュ形式は、配位自由度のパラメータの微分しか一般化速度として扱わない。しかし、SO(3)のような複雑な群を自由度にとっているときは、いちいち座標系を貼りたくない。座標系としてはオイラー角などがあるが、対称性が失われるし、ジンバルロックが発生することもある。そこで、リー代数およびリー群をそのまま一般化速度と配位自由度として使い、その座標を明示しないまま議論すること ができないか？と思えてくる。
まず、系の自由度はリー群なのだから、変分も当然リー群値で行う。変分リー群は微小なので、それはリー代数で近似できるはずである。\Theta(t)をその変分を与えるリー代数の値とする。
V(t) = (1+ \Theta(t)\epsilon)
そうすると変分を受けた軌道の微分(速度の変分)は
\begin{align*}
&\frac{d}{dt}V(t)R(t)\\
= &(1+\Theta(t)\epsilon)R(t)\Omega(t) + \dot{\Theta}(t)R(t)\epsilon\\
= &R(t)\Omega(t)+\epsilon\left(\Theta(t)R(t)\Omega(t) + \dot{\Theta}(t)R(t)\right)
\end{align*}
これをリー代数値に巻き戻すと、\epsilonの一次で
\begin{align*}
&R^{-1}(t)V^{-1}(t)\frac{d}{dt}V(t)R(t)\\
=& \Omega(t) + \epsilon R^{-1}(t)\dot{\Theta}(t)R(t)
\end{align*}
これを使ってラグランジアンを変分する。ただし、ラグランジアンは、R,\Omegaの関数になっているとする。つまり、配位自由度がRの状態にあり、回転運動がその点の接空間要素R\OmegaであるときのラグランジアンがL(R,\Omega)ということだ。
\begin{align*}
&\delta S \\
=&\epsilon \int dt \frac{\partial L}{\partial R}\Theta(t)R(t) + \frac{\partial L}{\partial \Omega}R^{-1}(t)\dot{\Theta}(t)R(t)\\
=&\epsilon \int dt \frac{\partial L}{\partial R}\Theta(t)R(t)-\frac{d}{dt}\frac{\partial L}{\partial \Omega} R^{-1} \Theta R - \frac{\partial L}{\partial \Omega} \left[R^{-1}\Theta R, \Omega \right]\\
=& \epsilon \int dt \frac{\partial L}{\partial R}RR^{-1}\Theta R-\frac{d}{dt}\frac{\partial L}{\partial \Omega} R^{-1} \Theta R + \frac{\partial L}{\partial \Omega} \mathrm{ad}_{\Omega}R^{-1} \Theta R
\end{align*}
変分項に比例する積分となったので、その係数を取り出す。これより、coordinate freeなオイラーラグランジュ方程式
\frac{d}{dt}\frac{\partial L}{\partial \Omega} - \mathrm{ad}^{\ast}_{\Omega} \frac{\partial L}{\partial \Omega} - \frac{\partial L}{\partial R}R = 0
が得られる。これはすでに知られており、オイラー=ポアンカレの方程式 と呼ばれているようだ。
この式の型が一見謎だが、これは\mathfrak{so}(3)^\astつまりリー代数の双対空間に値をとる方程式となる。例えば、\Omegaによる偏微分は、リー代数引数のラグランジアンをリー代数で微分している。Rによる微分が悩ましいが、この微分係数についてるRは必須である。Rで直接微分すると、Rの微小差分、つまりR\Omegaのような形を期待することになるが、この左移動Rをデフォルトにすることでリー代数の双対になるようにしている。
自由回転する剛体ラグランジアンについて具体的に考えてみよう。
L(R,\Omega) = \frac{1}{2}\mathrm{Tr}\Omega M \Omega^\top + m\mathbf{g}\cdot R\mathbf{r}_0
まず、リー代数は反対称行列のなす3次元空間だが、内積\frac{1}{2}\mathrm{Tr}(-^\top,-)を使うことで\mathfrak{so}(3)\sim \mathfrak{so}(3)^\astの同一視ができる。このときリー代数の基底はこの対応で自己双対になる。この同型を\mapstoであらわすと
\begin{align*}
\frac{\partial L}{\partial \Omega} = \mathrm{Tr}(M\Omega^\top -) \mapsto -\frac{1}{2}(\Omega M +M\Omega)
\end{align*}
第二項は
\begin{align*}
&\mathrm{ad}^{\ast}_{\Omega} \frac{\partial L}{\partial \Omega} (A)\\
= &\mathrm{Tr}(M\Omega^\top (\Omega A - A\Omega))
\end{align*}
なので
\mathrm{ad}^{\ast}_{\Omega} \frac{\partial L}{\partial \Omega} \mapsto \frac{1}{2}(\Omega\Omega M - M\Omega \Omega)
ポテンシャル項は
\frac{\partial L}{\partial R} R = m \mathrm{Tr}(\mathbf{r}_ 0\mathbf{g}^\top R-)\mapsto \frac{m}{2}( \mathbf{r}_0\mathbf{g}^\top R - R^{-1}\mathbf{g}\mathbf{r}_0^\top)
よって
-(\dot{\Omega}M + M\dot{\Omega}) - (\Omega\Omega M - M\Omega \Omega) = \frac{m}{2}( \mathbf{r}_0\mathbf{g}^\top R - R^{-1}\mathbf{g}\mathbf{r}_0^\top)
というリー代数値の方程式が立つ。これは
J = \frac{1}{2}(\Omega M + M\Omega)
としたら
\dot{J} + [\Omega,J] = - \frac{m}{2}( \mathbf{r}_0\mathbf{g}^\top R - R^{-1}\mathbf{g}\mathbf{r}_0^\top)
と同値である。I,Mがローカル座標系で対角化されている(かならず同時対角化できる)とき、
\begin{align*}
I = & \begin{bmatrix}
I_1 & 0 & 0\\
0 & I_2 & 0\\
0 & 0 & I_3
\end{bmatrix} = \begin{bmatrix}
M_2 + M_3 & 0 & 0\\
0 & M_3+M_1 & 0\\
0 & 0 & M_1 + M_2
\end{bmatrix} \\
M = & \begin{bmatrix}
M_1&0&0\\
0&M_2&0\\
0&0&M_3
\end{bmatrix}
\end{align*}
実はJはローカル座標系で見た時の角運動量そのものだ。ローカル座標系での角運動量ベクトルをj =I\omegaとしたときJ_{ij} = -\epsilon_{ijk}j_kということだ。
\begin{align*}
J = & \begin{bmatrix}
0 & -I_3 \omega_3 & I_2 \omega_2\\
I_3 \omega_3 & 0 & -I_1\omega_1\\
-I_2 \omega_2 & I_1 \omega_1 & 0
\end{bmatrix}
\end{align*}
したがってやはり完全反対称テンソルの計算により、方程式は
\begin{align*}
I\dot{\omega} + \omega \times I\omega = f(R)
\end{align*}
とオイラーの運動方程式に帰着する。成分毎に書くなら
\begin{align*}
I_1 \dot{\omega}_1 - \omega_2 \omega_3(I_2-I_3) &= f(R)_{1}\\
I_2 \dot{\omega}_2 - \omega_3 \omega_1(I_3-I_1) &= f(R)_{2}\\
I_3 \dot{\omega}_3 - \omega_1 \omega_2(I_1-I_2) &= f(R)_{3}
\end{align*}
ただしここで
\begin{align*}
F(R) =& \frac{m}{2}( \mathbf{r}_0\mathbf{g}^\top R - R^{-1}\mathbf{g}\mathbf{r}_0^\top)\\
F(R)_{ij} =& -\epsilon_{ijk}f(R)_{k}
\end{align*}
とおいた。
あとは\Omega(t)によって駆動されるR(t)の方程式
\frac{d}{dt}R(t) = R(t)\Omega(t)
を連立すればよい。シミュレーションをするなら、\Omega(t)の更新と(1+\Omega(t) dt)の右作用を繰り返すことになる。
Discussion