記号の世界ゟ

このブログでは, 数学書などの書評を書きます。また、受験などの勉強法をまとめます。

二体問題を解く話（四元数の応用）

可積分 数学 物理学

二体問題は解けて三体問題は解けないと言いますが，二体問題の解法の一つを紹介します．これは三体問題に関する次の記事の準備というのが最大のモチベーションです。せっかくなので以下の論文で紹介されている四元数を使った解法を説明します．
Jan Vrbik, "A novel solution to Kepler’s problem", 2003.

二体問題とケプラー問題

二体問題とケプラー問題の関係が曖昧になっている印象なので，まずはそれらを整理しておきます．

三次元空間の二つの質点の位置を v = (v_1, v_2, v_3), w = (w_1, w_2, w_3) とし，それぞれの質量を m_1, m_2 とします．万有引力定数を G とすると，二体問題の運動方程式は
\displaystyle \qquad m_1 \ddot{v} = - \frac{G m_1 m_2}{|v - w|^3} (v-w) \\ \displaystyle \qquad m_2 \ddot{w} = - \frac{G m_1 m_2}{|v - w|^3} (w-v)
です．

ここで，質点２から見た質点１の相対座標 r = v - w は
\displaystyle \qquad \ddot{r} = \ddot{v} - \ddot{w} \\ \displaystyle \qquad \, = \left( - \frac{G m_1 m_2}{m_1 |v - w|^3} (v-w) + \frac{G m_1 m_2}{m_2 |v - w|^3} (w-v) \right) \\ \displaystyle \qquad \, = - \left( \frac{G m_1 m_2}{m_1 |r|^3} r + \frac{G m_1 m_2}{m_2 |r|^3} r \right) \\ \displaystyle \qquad \, = - \frac{G m_1 m_2}{M |r|^3} r
ここで，
\displaystyle \qquad M = \left(\frac{1}{m_1} + \frac{1}{m_2} \right)^{-1} = \frac{m_1 m_2}{m_1 + m_2}
である．ここで， \mu = G m_1 m_2 / M と置けば定数は一つになる．このように導かれた以下の方程式をケプラー問題という．
\displaystyle \qquad \ddot{r} = -\frac{\mu}{ |r|^3} r

一方，質量中心
\displaystyle \qquad u = \frac{m_1 v + m_2 w}{m_1 + m_2}
をとると，
\displaystyle \qquad \ddot{u} = \frac{1}{m_1 + m_2} \left( m_1 \ddot{v} + m_2 \ddot{w} \right) \\ \displaystyle \qquad \, = \frac{1}{m_1 + m_2} \left( - \frac{G m_1 m_2}{|v - w|^3} (v-w) - \frac{G m_1 m_2}{|v - w|^3} (w-v) \right) \\ \qquad \, = 0
となるので，質量中心は等速直線運動する．

このことから，二体問題はケプラー問題の解を使って書けることが分かる．

（命題）
ケプラー問題の解を用いて二体問題は書ける．つまり，ケプラー問題が解ければ二体問題も解ける．

証明．質量中心 u = (m_1 v + m_2 w) / (m_1 + m_2) は等速直線運動するので，
\displaystyle \qquad u(t) = \frac{m_1 v(t) + m_2 w(t)}{m_1 + m_2} = u_0 + u_1 t
のように書ける．ここで u_0, u_1 は定数ベクトル．ここでケプラー問題の解を r(t) とすると， r(t) = v(t) - w(t) なので，計算すると
\displaystyle \qquad v(t) = \frac{m_2}{m_1 + m_2} r(t) + u(t) \\ \displaystyle \qquad w(t) = \frac{-m_1}{m_1 + m_2} r(t) + u(t)
のように二体問題の解を求めることができる．\square

ケプラー問題は，二体問題で一つの質点を固定した問題と見ることもできる．質量の差が大きい場合には二体問題をそのように捉えることは妥当である．しかし，二体問題の近似としてケプラー問題を考えて，ケプラー問題の性質から二体問題の解が楕円や双曲線になると述べてるものがあるが，そういう説明では論理に飛躍がある．上で述べたような，ケプラー問題の解を使って二体問題が書けるというのが正しい説明の仕方だと思う．

（もっと言えば，ケプラー問題は制限二体問題とも言える．二体問題と制限二体問題を同一視していい理由は自明ではない．例えば，三体問題と制限三体問題では結構違いが出てくる．）

四元数を使ったケプラー問題の解法

ケプラー問題を解こう．様々な解法や説明の仕方があると思うが，今回は四元数による説明をする．

四元数について

ここからは三次元ベクトルを表すのに太字 \mathbf{a} を用いて，四元数を表すのにブラックボード体 \mathbb{a} を用いる．また，実部を持たない四元数 \mathbb{a} = xi + yj +zk をベクトル \mathbf{a} = (x,y,z) と同一視する．

四元数の知識は既知とするが今回使うのはそれほど難しいものではない．複素数の場合， i^2 = -1 となる i を導入するが，それに加えて j^2 = k^2 = -1 となる j,k を加えたものが四元数である．ただし，積に関して可換ではなく， ij = k, jk = i, ki = j となる一方で ji = -k, kj = -i, ik = -j となる．慣れてなければ， i,j,k はこういう性質を満たす行列と思えば計算ミスが少なくなる．

行列と同様，指数も
\displaystyle \qquad e^{\mathbb{a} } = \sum_{l=0}^\infty \frac{1}{n!} \mathbb{a}^n
と定義する．特に，
\displaystyle \qquad e^{x i} = \cos (x) + i \sin (x) \\ \displaystyle \qquad e^{x j} = \cos (x) + j \sin (x) \\ \displaystyle \qquad e^{x k} = \cos (x) + k \sin (x) \\
となる．今回の内容では，指数は三角関数に直して計算すればうまくいくと思う． \mathbb{a}, \mathbb{b} が可換でなければ e^{\mathbb{a} } e^{\mathbb{b} } = e^{\mathbb{b} } e^{\mathbb{a} } とはならないことに注意せよ ．

共役とノルムを \bar{\mathbb{a}}, |\mathbb{a} | で表す．論文に合わせて， \mathbb{w} \bar{\mathbb{w}} = 1 を満たす四元数を回転四元数と呼ぶことにする．実部を持たない四元数，つまり，ベクトル \mathbf{u} を用いて，回転四元数を \mathbb{w} = e^{\mathbf{u}} と書くことができる．三次元での回転は回転四元数による変換
\quad \displaystyle \mathbf{r} \mapsto \mathbb{w} \mathbf{r} \bar{\mathbb{w}} = \exp (\mathbf{u} ) \mathbf{r} \exp (-\mathbf{u})
と書くことができる．

ケプラー問題を解く準備

さて，ケプラー問題
\displaystyle \qquad \ddot{\mathbf{r} } = -\frac{\mu}{ |\mathbf{r}|^3} \mathbf{r}
を四元数を用いて解いていこう．ポイントはベクトル \mathbf{r} を変数変換
\quad \displaystyle \mathbf{r} = \mathbb{u} k \bar{\mathbb{u} }
により四元数 \mathbb{u} と変換することである．右辺の共役をとると，
\quad \displaystyle \overline{\mathbb{u} k \bar{\mathbb{u} } } = \mathbb{u} \bar{k} \bar{\mathbb{u} } = - \mathbb{u} k \bar{\mathbb{u} }
となり，実部を持たないから，任意の \mathbb{u} で \mathbf{r} はちゃんとベクトルになっている．注意として，ベクトルは３次元で四元数は４次元だからこの変換は一対一ではない．別の言い方をすれば \mathbb{u} には１次元ぶん自由に決めることができる．この性質をうまく使う．

少し記号がややこしくなるが \displaystyle r = |\mathbf{r} | = \sqrt{\mathbf{r} \bar{\mathbf{r}} } と表す．計算すれば， r = \mathbb{u} \bar{\mathbb{u}} = \bar{\mathbb{u} } \mathbb{u} であることも分かる．（計算の途中で， \mathbb{u} \bar{\mathbb{u}} が実数だから k と可換であるという性質を使う．）

非可換なので r を微分してみると，
\quad \displaystyle \dot{\mathbf{r}} = \dot{\mathbb{u}} k \bar{\mathbb{u}} + \mathbb{u} k \dot{ \bar{\mathbb{u}} }
となり，微分方程式に代入してみても非常に複雑な式が出る．そこで，
\quad \displaystyle \Gamma = \mathbb{u} k \dot{ \bar{\mathbb{u}} }-\dot{\mathbb{u}} k \bar{\mathbb{u}}
とおくと，
\quad \displaystyle \dot{\mathbf{r}} = 2\dot{\mathbb{u}} k \bar{\mathbb{u}} + \Gamma
となる． \mathbb{u} の定め方に不定性があったので， \Gamma = 0 となるようにすれば，
\quad \displaystyle \dot{\mathbf{r}} = 2\dot{\mathbb{u}} k \bar{\mathbb{u}} \, (= 2 \mathbb{u} k \dot{ \bar{\mathbb{u}} } )
となり，計算が簡単になる． \Gamma をゲージという．さらに， t も
\quad \displaystyle \frac{d t}{d s} = 2r \sqrt{\frac{a}{\mu} }
により s に変換する．ここで a は初期値に応じて定める定数である．以下では s による微分を \mathbb{u}' のようにプライムで表示する．するとケプラー問題は以下のように書ける．

（命題）
ゲージ \Gamma = 0 の下でケプラー問題
\displaystyle \qquad \ddot{\mathbf{r} } = -\frac{\mu}{ |\mathbf{r}|^3} \mathbf{r}
は
\displaystyle \qquad \mathbb{u}'' - \mathbb{u} \left( \frac{ \bar{\mathbb{u} }' \mathbb{u}' -2a}{r} \right) = 0
と書ける．

証明． \dot{\mathbf{r}} = 2\dot{\mathbb{u}} k \bar{\mathbb{u}} の両辺に右から 2\sqrt{a/\mu} u k をかけると，
\displaystyle \qquad 2 \sqrt{\frac{a}{\mu}} \dot{r} u k = - 2 \sqrt{\frac{a}{\mu}} 2 \dot{u} r = - 2 u'
となる．さらにこれを s で微分すると，
\displaystyle \qquad -2 \mathbb{u}'' = 2 \sqrt{\frac{a}{\mu}} r \frac{d}{dt} \left(2 \sqrt{\frac{a}{\mu}} \dot{r} \mathbb{u} k \right) \\ \displaystyle \qquad \qquad =4r \frac{a}{\mu} ( \ddot{\mathbf{r} } \mathbb{u} k + \dot{\mathbf{r}} \dot{\mathbb{u} } k) \\ \displaystyle \qquad \qquad =4a \frac{\mathbb{ u } }{r^2} + \frac{2}{r} \mathbb{u}' k \bar{\mathbb{u} } \mathbb{u}' k
ここで，ケプラー問題の方程式を用いた．さらに \Gamma = 0 の式を用いると最後の項は
\displaystyle \qquad \frac{2}{r} \mathbb{u}' k \bar{\mathbb{u} } \mathbb{u}' k = \frac{2}{r} \mathbb{u} k \bar{\mathbb{u} }' \mathbb{u}' k = - \frac{2}{r} \mathbb{u} \bar{\mathbb{u} }' \mathbb{u}'
となる．これを用いて整理すると，
\displaystyle \qquad \mathbb{u}'' - \mathbb{u} \left( \frac{ \bar{\mathbb{u} }' \mathbb{u}' -2a}{r} \right) = 0
が分かる． \square

ここで登場した
\displaystyle \qquad \frac{ \bar{\mathbb{u} }' \mathbb{u}' -2a}{r}
を微分すると０になるので一定であることが分かるが，実は計算するとハミルトニアンの定数倍になっていることも分かる．

（命題）
ケプラー問題のハミルトニアン
\displaystyle \qquad H = \frac{ |\dot{\mathbf{r}} |^2}{2} - \frac{\mu}{r}
を用いれば，
\displaystyle \qquad \frac{ \bar{\mathbb{u} }' \mathbb{u}' -2a}{r} = 2a \mu H
と書ける．
証明．
\displaystyle \qquad \frac{ \bar{\mathbb{u} }' \mathbb{u}' }{r} = 4\frac{a}{\mu} r\bar{\dot{ \mathbb{u} } } \dot{ \mathbb{u}}
であり，
\displaystyle \qquad |\dot{\mathbf{r} }|^2 = \dot{\bar{\mathbf{r} } } \dot{\mathbf{r}} =- \dot{\mathbf{r}} \dot{\mathbf{r}} = -4 \dot{\mathbb{u}} k \bar{\mathbb{u}} \dot{\mathbb{u}} k \bar{\mathbb{u}} = 4 r \dot{\mathbb{u}} \dot{ \bar{\mathbb{u}} } = 4 r \dot{ \bar{\mathbb{u}} } \dot{\mathbb{u}}
であることから証明できる． \square

まとめると以下が示された．

（定理）
ゲージ \Gamma = 0 の下でケプラー問題
\displaystyle \qquad \ddot{\mathbf{r} } = -\frac{\mu}{ |\mathbf{r}|^3} \mathbf{r}
は
\displaystyle \qquad \mathbb{u}'' - 4 a \mu H \mathbb{u} = 0
と書ける．

楕円軌道の場合の解法

線形微分方程式に帰着されたので，解自体は簡単に求まる．問題は初期値に応じてうまく a を決めることである． a は s の定義の平方根の中に入っているため， a > 0 としなければならない．

・ H < 0 となる初期値の場合， -4 a \mu H = 1 となるように a を定めれば，微分方程式は
\displaystyle \qquad \mathbb{u}'' + \mathbb{u} = 0
だから，解は
\displaystyle \qquad \mathbb{u} = \mathbb{p} \cos s + \mathbb{q} \sin s
と書ける．
・ H = 0 となる初期値の場合，微分方程式は
\displaystyle \qquad \mathbb{u}'' = 0
だから，解は
\displaystyle \qquad \mathbb{u} = \mathbb{p} + \mathbb{q} s
と書ける．
・ H > 0 となる初期値の場合， 4 a \mu H = 1 となるように a を定めれば，微分方程式は
\displaystyle \qquad \mathbb{u}'' - \mathbb{u} = 0
だから，解は
\displaystyle \qquad \mathbb{u} = \mathbb{p} \cosh s + \mathbb{q} \sinh s
と書ける．

これから \mathbf{r} を求めれば解は完全に分かる．ただし，ゲージが \Gamma = 0 となるように四元数の初期値を定める必要がある．さらに，この四元数による表示では実際の \mathbf{r} は複雑な形になってしまうため，それほど嬉しさを感じられない．（解けることはこの表示で一目瞭然ではあるけども）

そこで， H < 0 の場合に上でお見せした解の表示
\displaystyle \qquad \mathbb{u} = \mathbb{p} \cos s + \mathbb{q} \sin s
よりも意味がはっきりとする表示を与える．

（定理）
四元数 \mathbb{p}, \mathbb{q} を任意に取って定まる s の関数
\displaystyle \qquad \mathbb{u} = \mathbb{p} \cos s + \mathbb{q} \sin s
は
\displaystyle \qquad \mathbb{u} = \mathbb{w} \left[ A \exp ( i ( s - s_p) ) (1 + j \gamma) + B \exp ( - i (s - s_p)) \right] \exp (k \delta)
と書ける．ここで， A, B, \delta, \gamma, s_p は実数のパラメータで， \mathbb{w} は回転四元数である．

単純な計算で証明ができるがヒントなしではかなり難しいのでおまけの節で証明する．

さて，この表示
\displaystyle \qquad \mathbb{u} = \mathbb{w} \left[ A \exp ( i ( s - s_p) ) (1 + j \gamma) + B \exp ( - i (s - s_p)) \right] \exp (k \delta)
をとりあえず楕円表示とでも呼ぶことにして，この表示の意味を考えよう．

まず，ゲージ \Gamma = 0 が成り立たないとケプラー問題の解ではないため，この表示でゲージが０となるためのパラメータの条件を与えよう．

\displaystyle \qquad \mathbb{u}' k \mathbb{u} \\ \displaystyle \qquad = \mathbb{w} \left[ A \exp ( i (s- s_p) ) i ( 1 + j \gamma) + B \exp ( -i (s-s_p) ) (-i) \right] \\ \displaystyle \qquad \qquad \qquad \cdot k \left[ A ( 1 - j \gamma) \exp ( -i (s- s_p) ) + B \exp ( i (s-s_p) ) \right] \bar{\mathbf{w} } \\ \displaystyle \qquad = \mathbb{w} \left[ -A^2 (1+\gamma^2)\exp(2 i(s-s_p) ) + B^2 \exp(-2 i (s-s_p) )+2AB \cos(s - s_p) \right] j \bar{ \mathbb{w} } - 2 A \gamma
ここで，実数 \xi に対して ke^{i \xi} = e^{-i \xi} k となることを用いている．一方
\displaystyle \qquad \mathbb{u} k \mathbb{u}' \\ \displaystyle \qquad = - \overline{\mathbb{u}' k \mathbb{u}　} \\ \displaystyle \qquad = - \mathbb{w} (-j) \left[ -A^2 (1+\gamma^2)\exp( - 2 i(s-s_p) ) + B^2 \exp(2 i (s-s_p) )+2AB \cos(s - s_p) \right] \bar{ \mathbb{w} } + 2 A \gamma \\ \displaystyle \qquad = \mathbb{w} \left[ -A^2 (1+\gamma^2)\exp( 2 i(s-s_p) ) + B^2 \exp( -2 i (s-s_p) )+2AB \cos(s - s_p) \right] j \bar{ \mathbb{w} } + 2 A \gamma
となるので，
\displaystyle \qquad \mathbb{u}' k \mathbb{u} - \mathbb{u} k \mathbb{u}' = -4A\gamma
となる．よって， \Gamma = 0 となるには A\gamma = 0 とすればよい．

（命題）
楕円表示において \gamma = 0 と置けば， \Gamma = 0 となる．

これにより，ケプラー問題の解が得られる． \mathbf{r} = \mathbb{u} k\bar{\mathbb{u} } の計算を見てみれば，真ん中にある計算で e^{k \delta} k e^{-k \delta} = k となるため，実は \delta は \mathbf{r} では現れないパラメータである．よって， \delta = 0 としてよい．

最後に，
\displaystyle \qquad \frac{ \bar{\mathbb{u} }' \mathbb{u}' -2a}{r}
の値が s = s_p で
\displaystyle \qquad \frac{ (A-B)^2 -2a}{(A+B)^2}
となることが簡単な計算で分かる．この値が -1 となるように a を定めたのだから， A^2 + B^2 = a となっている．

以上をまとめると，楕円表示は
\displaystyle \qquad \mathbb{u} = \mathbb{w} \left[ A \exp ( i ( s - s_p) ) + B \exp ( - i (s - s_p)) \right]
となるが， \beta = B/A とおき，A をくくり出すと，
\displaystyle \qquad \mathbb{u} = \mathbb{w} \left[ \exp ( i ( s - s_p) ) + \beta \exp ( - i (s - s_p)) \right] \sqrt{\frac{a}{1 + \beta^2} }
これを使って \mathbf{r} を計算すると，
\displaystyle \qquad \mathbf{r} = \mathbb{u} k \bar{\mathbb{u} } \\ \displaystyle \qquad \, = \mathbb{w} \left[ \exp ( i ( s - s_p) ) + \beta \exp ( - i (s - s_p) ) \right] k \frac{a}{1 + \beta^2} \left[ \exp ( -i ( s - s_p) ) + B \exp ( i (s - s_p) ) \right] \bar{\mathbb{w} } \\ \displaystyle \qquad \, = \mathbb{w} \left[ \exp ( i ( s - s_p) ) + \beta \exp ( - i (s - s_p) ) \right]^2 \frac{a k}{1 + \beta^2} \bar{\mathbb{w} } \\ \displaystyle \qquad \, = \mathbb{w} \left[ \exp (2 i ( s - s_p) ) + \beta^2 \exp ( - 2 i (s - s_p) ) +2 \beta \right] \frac{a k}{1 + \beta^2} \bar{\mathbb{w} } \\ \displaystyle \qquad \, = \mathbb{w} \left[- (1-\beta^2) \sin 2 (s -s_p) j + \{(1+\beta^2) \cos 2(s-s_p) +2\beta \} k\right] \frac{a}{1 + \beta^2} \bar{\mathbb{w} } \\
となる．つまり， y-z 平面で
\displaystyle \qquad (y,z) = \left(-\frac{a}{1 + \beta^2} (1-\beta^2) \sin 2 (s -s_p) , \frac{a}{1 + \beta^2}\{(1+\beta^2) \cos 2(s-s_p) +2\beta \} \right)
と動くものを \mathbb{w} で回転させたものが解になっている． x-y 平面ではなく y -z 平面になっているところや， (\cos, \sin) ではなく (\sin, \cos) のような解の表示になっているところは不満が残るが，論文では少し変わった四元数の記法を使っているためそれが解消されている．（だから本当は楕円表示もこの定式化にあったもっといい表現方法があるはずである．）

少なくともこの表示により，解が楕円軌道上で動くことやその形や離心率も計算することができる．

まとめ

この記事の目的は，二体問題がケプラー問題に帰結することとケプラー問題が解けることを確認することでした．ケプラー問題は初期値に応じて軌道が大きく変わるわけですが，その解法が変わるのが少し問題ではあります．

本当は論文にはない双曲線の場合をやりたかったのですが，なかなか難しいです．なにが問題かというと，上で導入した楕円表示のようなものを双曲線の場合にも作る必要があるところです．安直に \cos, \sin を \cosh, \sinh に変更してもうまくいきそうにありませんでした．線形微分方程式の全ての解が表せるだけのパラメータを持っていて，かつパラメータに意味があるように表示する必要があります．例えば，楕円表示の場合には \gamma = 0 がゲージが０になることと対応していたのでした．答えから逆算すればいい表示が見つけられるかもしれません。

おまけ

疲れたので時間がある時に書きます．