力学の例題

惑星の運動

太陽の周りをまわる惑星の運動は、太陽を焦点の一つとする楕円軌道を描くことが知られている（ケプラーの第１法則）。 太陽を原点とする平面内で惑星の軌道を極座標 \( (r, \phi) \) で記述すると、動径 \( r \) は \( \phi \) の関数として \[ r(\phi) = \frac{l}{1 + \varepsilon \cos \phi}, \] と表される。 \( \varepsilon \) は離心率、\( l = a (1 - \varepsilon^2 ) \)、また、\( a \) は楕円の長軸半径である。
時々刻々変化する惑星の位置は、楕円軌道上を、 面積速度一定の法則（ケプラーの第２法則）「太陽と惑星を結ぶ線が単位時間に掃過する面積は一定である」、 に従って運動する。

ここまではよくある説明。 では、実際に時々刻々変化する運動の様子をどう与えるか。

様々な方法があるでしょうが、ここでは、以下の、ケプラーの方程式と呼ばれる、時刻 \( t \) と、離心近点離角 \( \theta \) の関係式 \[ \frac{2 \pi}{T} t = \theta - \varepsilon \sin \theta \] を用いて、運動の様子を求めます（ケプラーの方程式に関する説明は他をご参照下さい）。 なお、 動径 \( r \) と \( \theta \) の間には \[ r(\theta) = a ( 1 - \varepsilon \cos \theta ) \] の関係がある。 また、時刻 \( t = 0 \) において惑星は 近日点（ \( r = a ( 1 - \varepsilon ) \). このとき \( \phi = 0, \theta = 0. \) ） にある。

以下では、上記方程式を解くことで得られる楕円運動をアニメーションで示す。
\( 0 \le \varepsilon < 1 \) :

数値計算を行う上で必要な詳細情報を追記

どのように時々刻々と変化する惑星の位置を求めるかはほぼ上に書いた通り。 作成したコードでは、 各時刻 \( t \) における \( \theta \) を求めた後、 以下の式 \[ \phi = \arccos \left( \frac{ \cos \theta - \varepsilon }{1 - \varepsilon \cos \theta} \right) \] により \( \phi \) を求め、これを用いて惑星の位置を定めている。 この式は、上記の \( r(\phi) \) と \( r(\theta) \) の式を等しいとおくことより得られる。

数値計算で \( \arccos \) の値域が \( [0, \pi ] \) であることより、 \( \pi < \phi < 2 \pi \) の領域では、算出された結果を \( \phi \rightarrow 2 \pi - \phi \) と補正している。
また、\( 0 \le t < T \) で一周するので、元の位置に戻ってきたら、 \( t \rightarrow t - T \) と処理している。

後記

\( \theta \) は \( t \) の簡単な関数で解析的には表されていないので、 上記ケプラーの方程式を解くことで、各 \( t \) における \( \theta \) を数値的に得ています。 方程式の解法に関してよさげなライブラリでもないかと検索してみるが、使いやすそうなものが見つからず、 結局自分で作成することに。ケプラーの方程式の右辺は \( \theta \) の単調増加関数なので、 単純に２分法を用いました。
慣れていない言語でのプログラミングは、文法を調べたりバグとりやら、あれこれなんだかんだでいろいろと、 やはり結構時間かかります。（2020/04/22 記）
離心率の選択範囲のバグを修正。（2021/04/19）
閲覧して下さった方からのコメントをきっかけにページを読み直したところ、私自身でもわかりにくい箇所が散見していたので、加筆修正。（2022/10/12）
（2022/10/12 最終更新）

Copyright © Wataru Izumida, All Rights Reserved.