力学の例題

惑星の運動

太陽の周りをまわる惑星の運動は、太陽を焦点の一つとする楕円軌道を描くことが知られている（ケプラーの第１法則）。 太陽を原点とする平面内で惑星の軌道を極座標 \( (r, \phi) \) で記述すると、動径 \( r \) は \( \phi \) の関数として \[ r(\phi) = \frac{l}{1 + \varepsilon \cos \phi}, \] と表される。 \( \varepsilon \) は離心率、\( l = a (1 - \varepsilon^2 ) \)、また、\( a \) は楕円の長軸半径である。 時々刻々変化する惑星の位置は、 面積速度一定の法則（ケプラーの第２法則）「太陽と惑星を結ぶ線が単位時間に掃過する面積は一定である」、 および、ケプラーの第３法則（調和の法則）「惑星の公転周期（\( T \)）の2乗は、軌道長半径（\( a \)）の3乗に比例する」 に従う。 すなわち、公転周期（\( T \)）が与えられると、面積速度一定の法則より、 時刻 \( t \) における惑星の位置が定められる。 これは、 動径 \( r \) を \[ r(\theta) = a ( 1 - \varepsilon \cos \theta ) \] と記述する変数 \( \theta \) を用いて、 ケプラーの方程式 \[ \frac{2 \pi}{T} t = \theta - \varepsilon \sin \theta \] により定められる。 なお、時刻 \( t = 0 \) において惑星は 近日点（ \( r = a ( 1 - \varepsilon ) \). このとき \( \phi = 0, \theta = 0. \) ） にある。

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

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

どのように時々刻々と変化する惑星の位置を求めるかはほぼ上に書いてあるとおりだが、 各時刻 \( t \) における \( \theta \) を求めた後は 以下の式 \[ \phi = \arccos \left( \frac{ \cos \phi - \varepsilon }{1 - \varepsilon \cos \phi} \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）
（2021/04/19 最終更新）

Copyright © Wataru Izumida, All Rights Reserved.