diffeq2nd_planet_generator.py

プログラムの動作

このプログラムは、複数の惑星間の万有引力による軌道運動をシミュレーションするツールです。 NumPyを利用した高速な配列演算・ベクトル化処理を用いており、Pythonのジェネレータ機能（yield）を活用することで、逐次計算とデータ出力を並行して行う設計となっています。 Matplotlibを使用したリアルタイムでの軌道描画機能を備えており、指定された数値積分アルゴリズム（オイラー法、ホイン法、ベルレ法、速度ベルレ法）を切り替えることで、計算精度の違いやエネルギー・運動量の保存則の検証を行うことができます。

原理

本シミュレーションでは、ニュートンの運動方程式と万有引力の法則を用いて各惑星の加速度を計算します。惑星 \(i\) と惑星 \(j\) の間に働く引力による加速度 \(a_i\) は、以下の式で計算されます。

\[ a_i = \sum_{j \neq i} G_1 \frac{m_j}{|r_j - r_i|^3} (r_j - r_i) \]

ここで、\(G_1\) は計算用（単位系）にスケーリングされた万有引力定数、\(m_j\) は相手の惑星の質量、\(r_i, r_j\) はそれぞれの位置ベクトルです。

系の運動エネルギー \(K\) およびポテンシャルエネルギー \(U\) は次のように計算されます。

\[ K = \frac{1}{2} \sum_{i} m_i |v_i|^2 \]

\[ U = -\sum_{i < j} G_1 \frac{m_i m_j}{|r_j - r_i|} \]

力学的エネルギー \(E = K + U\) と全運動量 \(P = \sum m_i v_i\) の時間変化を追跡することで、選択した数値計算手法の安定性を評価します。時間積分には、オイラー法 (Euler)、ホイン法 (Heun)、ベルレ法 (Verlet)、速度ベルレ法 (vverlet) のいずれかを使用します。

必要な非標準ライブラリとインストール方法

本プログラムを実行するには、以下の非標準ライブラリが必要です。

numpy: 配列計算およびベクトル化された数値演算を行うため
matplotlib: リアルタイムで惑星の軌道を描画するため

これらのライブラリは、以下のコマンドを用いてインストールしてください。

pip install numpy matplotlib

必要な入力ファイル

シミュレーションの初期条件として、プログラムの実行ディレクトリに以下のファイルが必要です。

ファイル名: planet_db.csv

このCSVファイルはヘッダ行を持ち、各行に惑星のデータを記述します。少なくとも以下の列が含まれている必要があります。

Name: 惑星の名前（文字列）
Mass: 惑星の質量
Revolution Radius: 公転半径
Revolution Velocity: 公転速度

※ Name 列の文字列データ以外は、読み込み時に自動で浮動小数点数（float）に変換されてシミュレーションに利用されます。

生成される出力ファイル

プログラムの実行が完了すると、カレントディレクトリに以下の2つのCSVファイルが生成・保存されます（ファイル名の {solver} 部分には、使用したアルゴリズム名が小文字で入ります）。

diffeq2nd_Planet_{solver}.csv
- 内容: 各ステップの時間 \(t\) と、各惑星の二次元座標 \((x, y)\) の時間履歴。
diffeq2nd_Planet_{solver}_conservation.csv
- 内容: 時間 \(t\) に対する系のポテンシャルエネルギー \(U\)、運動エネルギー \(K\)、全力学的エネルギー \(E\)、全運動量の各成分 \(Px, Py, Pz\)、および運動量の二乗平均平方根 \(Pmsm\) の履歴。エネルギーや運動量の保存状態を確認するために用います。

コマンドラインでの使用例 (Usage)

コマンドライン引数を用いて、シミュレーションの各種パラメータを指定できます。引数は以下の順序で処理されます（引数を省略した場合はソースコード内のデフォルト値が適用されます）。

python D:/git/sphinx/tkProg/source/tiny_simulations/diffeq2nd_planet_generator.py [solver] [dt] [nt] [fplot] [yield_every] [enable_ctrlc]

solver : 数値積分アルゴリズム（Euler, Heun, Verlet, vverlet）。デフォルトは vverlet。
dt : 時間刻み幅。デフォルトは 0.1。
nt : 総シミュレーションステップ数。デフォルトは 20000。
fplot : 軌道のリアルタイム描画を有効にするか（1で有効、0で無効）。デフォルトは 1。
yield_every : 指定ステップごとにジェネレータから結果を出力し、描画を更新する。デフォルトは 10。
enable_ctrlc : Ctrl+C によるプログラムの中断を許可するか（1で許可、0で無視）。デフォルトは 1。

コマンドラインでの具体的な使用例

以下は、ホイン法（Heun）を用いて、時間刻みを 0.05、総ステップを 50000 とし、描画を有効にしたうえで、20ステップごとの出力更新を行い、Ctrl+C による中断も許可する実行例です。

python D:/git/sphinx/tkProg/source/tiny_simulations/diffeq2nd_planet_generator.py Heun 0.05 50000 1 20 1

実行結果の説明:

コマンドを実行すると、まず初期化情報やパラメータ、読み込まれた惑星のリストがコンソールの標準出力に表示されます。その後、fplot=1 が指定されているため Matplotlib のウィンドウが開き、惑星の軌道がリアルタイムで描画されます。コンソール上には一定のステップ間隔で現在の時間 \(t\) と各惑星の座標値が出力されていきます。指定した 50000 ステップの計算が完了するか、あるいはユーザーがターミナルで Ctrl+C を押して中断すると、プログラムの経過時間（Elapsed time）が表示され、「Press ENTER to exit>>」というプロンプトを出して終了待機状態になります。 Enterキーを押して終了すると、実行ディレクトリに diffeq2nd_Planet_heun.csv および diffeq2nd_Planet_heun_conservation.csv の2つのファイルが保存されています。