diffeq2nd_planet.py 技術ドキュメント

プログラムの動作

diffeq2nd_planet.py は、複数の天体（惑星）間の重力相互作用に基づく運動をシミュレーションするためのPythonプログラムです。主に太陽系の惑星の軌道をシミュレートすることを想定しています。

主な機能:

多体問題の解決: 複数の惑星間にかかる万有引力を計算し、それによって生じる運動を数値的に解きます。
数値積分アルゴリズムの選択: 運動方程式を解くための数値積分アルゴリズムとして、オイラー法 (Euler method) またはベレ法 (Verlet method) を選択できます。
惑星データの読み込み: 惑星の質量、初期位置、初期速度などのデータは外部のCSVファイルから読み込みます。
軌道データの出力: 各惑星の時系列での位置 (x, y) データをCSVファイルに出力します。
保存量の出力: システム全体の運動エネルギー、ポテンシャルエネルギー、総エネルギー、および総運動量をCSVファイルに出力し、保存則の検証を可能にします。
リアルタイムグラフ表示: シミュレーションの進行と同時に、各惑星の軌跡をmatplotlibを用いてリアルタイムでグラフ表示します（オプション）。
コマンドライン引数による制御: シミュレーションのソルバー、時間刻み、総ステップ数、グラフ表示の有無をコマンドライン引数で指定できます。

解決する課題:

このプログラムは、古典的な物理学における多体問題（特に天体力学）を数値的に解析する手段を提供します。異なる数値積分アルゴリズムがシミュレーションの安定性や精度にどのように影響するかを比較・検討することができます。

原理

このプログラムは、ニュートンの万有引力の法則と運動の法則に基づいて、各惑星の運動をシミュレートします。

1. ニュートンの万有引力の法則

質量 \(M_i\) の惑星 \(i\) と質量 \(M_j\) の惑星 \(j\) の間に働く万有引力 \(\vec{F}_{ij}\) は、二つの惑星間の距離 \(r_{ij}\) を用いて以下のように表されます。

\[\vec{F}_{ij} = G \frac{M_i M_j}{r_{ij}^2} \frac{\vec{r}_j - \vec{r}_i}{|\vec{r}_j - \vec{r}_i|}\]

ここで \(G\) は万有引力定数です。惑星 \(i\) にかかる全重力は、他のすべての惑星からの引力の総和となります。

\[\vec{F}_i = \sum_{j \neq i} \vec{F}_{ij}\]

2. 運動方程式

惑星 \(i\) の位置ベクトルを \(\vec{r}_i\)、速度ベクトルを \(\vec{v}_i\)、加速度ベクトルを \(\vec{a}_i\) とすると、運動方程式は以下のようになります。

\[M_i \vec{a}_i = M_i \frac{d^2\vec{r}_i}{dt^2} = \vec{F}_i\]

これを加速度について解くと、

\[\vec{a}_i = \frac{d^2\vec{r}_i}{dt^2} = \frac{1}{M_i} \sum_{j \neq i} G \frac{M_i M_j}{r_{ij}^2} \frac{\vec{r}_j - \vec{r}_i}{|\vec{r}_j - \vec{r}_i|} = \sum_{j \neq i} G \frac{M_j}{r_{ij}^2} \frac{\vec{r}_j - \vec{r}_i}{|\vec{r}_j - \vec{r}_i|}\]

プログラム内では、加速度の計算には Fi 関数が用いられ、引力定数 \(G\) は単位変換された \(G_1\) を使用します。

\[G_1 = G \cdot \left(\frac{\text{DayToSecond}}{\text{AstronomicalUnit}}\right)^2 \cdot \frac{1}{\text{AstronomicalUnit}}\]

この \(G_1\) を用いることで、距離の単位をAU (天文単位)、時間の単位を日 (Day) として計算を行うことができます。

3. 数値積分アルゴリズム

2階の常微分方程式を数値的に解くために、以下のいずれかのアルゴリズムを使用します。

a) オイラー法 (Euler method)

最も基本的な数値積分法です。現在の位置と速度、加速度を用いて次のステップの位置と速度を計算します。

位置: \(\vec{r}(t+\Delta t) = \vec{r}(t) + \vec{v}(t)\Delta t\)
速度: \(\vec{v}(t+\Delta t) = \vec{v}(t) + \vec{a}(t)\Delta t\)

b) ベレ法 (Verlet method)

エネルギー保存性に優れるとされる数値積分法です。現在の位置と過去の位置、現在の加速度を用いて次のステップの位置を計算し、その結果から速度を推定します。

位置: \(\vec{r}(t+\Delta t) = 2\vec{r}(t) - \vec{r}(t-\Delta t) + \vec{a}(t)(\Delta t)^2\)
速度: \(\vec{v}(t+\Delta t) = \frac{\vec{r}(t+\Delta t) - \vec{r}(t-\Delta t)}{2\Delta t}\)

4. 保存量

シミュレーションの精度や安定性を評価するために、システム全体の総エネルギーと総運動量の保存を確認します。

a) 総運動量 (Total Momentum)

システム全体の総運動量 \(\vec{P}\) は、各惑星の運動量の総和です。

\[\vec{P} = \sum_i M_i \vec{v}_i\]

プログラムでは、初期段階でシステム全体の総運動量がゼロになるように各惑星の速度を調整 (normalize_momentum 関数) しています。

b) 総エネルギー (Total Energy)

システム全体の総エネルギー \(E\) は、運動エネルギー \(K\) とポテンシャルエネルギー \(U\) の和です。

運動エネルギー: \(K = \sum_i \frac{1}{2} M_i |\vec{v}_i|^2\)
ポテンシャルエネルギー: プログラムの Utot 関数内では、引力ポテンシャルエネルギーの符号が反転した値、すなわち \(U_{code} = \sum_{i<j} G_1 \frac{M_i M_j}{r_{ij}}\) として計算されます。
総エネルギー: プログラムでは、総エネルギー \(E\) は \(K + U_{code}\) として計算されます。これは、物理的なポテンシャルエネルギーが負の値であることを考慮すると、見かけ上は \(E = K - |U_{pot}|\) となりますが、ここではプログラムの定義に従います。

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

diffeq2nd_planet.py の実行には、以下の非標準ライブラリが必要です。

numpy: 数値計算を効率的に行うためのライブラリです。
matplotlib: データの可視化、特にグラフ描画に用いられます。

これらのライブラリは、Pythonのパッケージマネージャーである pip を使用してインストールできます。

pip install numpy matplotlib

必要な入力ファイル

プログラムは、惑星の初期条件を記述したCSVファイルを読み込みます。デフォルトのファイル名は planet_db.csv です。

ファイル形式:

ファイルはヘッダ行を含み、各行が1つの惑星のデータに対応します。以下の列を持つ必要があります。

Name: 惑星の名前 (文字列)
Mass: 惑星の質量 (kg)
Revolution Radius: 初期公転半径 (m)
Revolution Velocity: 初期公転速度 (m/s)

planet_db.csv の例:

Name,Mass,Revolution Radius,Revolution Velocity
Sun,1.989e+30,0.0,0.0
Mercury,3.3011e+23,5.790905e+10,4.7362e+04
Venus,4.8675e+24,1.08208e+11,3.502e+04
Earth,5.9722e+24,1.49598023e+11,2.978e+04
Mars,6.4171e+23,2.279392e+11,2.413e+04
Jupiter,1.8982e+27,7.7857e+11,1.307e+04
Saturn,5.6834e+26,1.43353e+12,9.69e+03
Uranus,8.6813e+25,2.8772e+12,6.81e+03
Neptune,1.02413e+26,4.4983964e+12,5.43e+03

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

プログラムは、シミュレーション結果を2つのCSVファイルに出力します。ファイル名は選択されたソルバーによって動的に決定されます。

軌道データファイル: diffeq2nd_Planet_{solver}.csv
- {solver} は Euler または Verlet のいずれか。
- 内容: 各時間ステップにおける各惑星のX座標とY座標。
  - t: 時間 (日)
  - x(PlanetName): 惑星のX座標 (AU)
  - y(PlanetName): 惑星のY座標 (AU)
  - 例: t,x(Sun),y(Sun),x(Mercury),y(Mercury),...
保存量データファイル: diffeq2nd_Planet_{solver}_conservation.csv
- {solver} は Euler または Verlet のいずれか。
- 内容: 各時間ステップにおけるシステム全体の運動エネルギー、ポテンシャルエネルギー、総エネルギー、および総運動量の成分。
  - t: 時間 (日)
  - U: ポテンシャルエネルギー（プログラム内で定義された値）
  - K: 運動エネルギー
  - E: 総エネルギー (U + K)
  - Px: 総運動量のX成分
  - Py: 総運動量のY成分
  - Pz: 総運動量のZ成分
  - Pmsm: 各惑星の運動量の二乗和の平均の平方根に相当する値

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

diffeq2nd_planet.py は、以下の書式でコマンドラインから実行できます。

python diffeq2nd_planet.py [solver] [dt] [nt] [fplot]

引数:

solver (文字列, オプション): 使用する数値積分アルゴリズム。
- 'Euler' (デフォルト) または 'Verlet'。
dt (浮動小数点数, オプション): 時間刻み幅 (日)。
- デフォルトは 0.1 (日)。
nt (整数, オプション): 計算する総ステップ数。
- デフォルトは 20000 ステップ。
fplot (整数, オプション): グラフ表示のフラグ。
- 0: グラフを表示しない。
- 1: グラフを表示する (デフォルト)。

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

1. デフォルト設定での実行

ソルバーにオイラー法、時間刻み 0.1 日、総ステップ数 20000、グラフ表示ありで実行します。

python diffeq2nd_planet.py

実行結果:

コンソールにシミュレーションの進行状況（各ステップでの時間と主要惑星の座標）が表示されます。
diffeq2nd_Planet_Euler.csv ファイルが生成され、各惑星の軌道データが保存されます。
diffeq2nd_Planet_Euler_conservation.csv ファイルが生成され、エネルギーと運動量の保存データが保存されます。
matplotlibによる軌道グラフがリアルタイムで表示され、シミュレーションが終了するとグラフウィンドウが表示されたまま停止し、エンターキーを押すことで終了します。

2. ベレ法を使用し、時間刻みを細かく、ステップ数を増やしてグラフ表示ありで実行

ソルバーにベレ法、時間刻み 0.05 日、総ステップ数 40000、グラフ表示ありで実行します。

python diffeq2nd_planet.py Verlet 0.05 40000 1

実行結果:

コンソールにシミュレーションの進行状況が表示されます。
diffeq2nd_Planet_Verlet.csv ファイルが生成され、軌道データが保存されます。
diffeq2nd_Planet_Verlet_conservation.csv ファイルが生成され、エネルギーと運動量の保存データが保存されます。
matplotlibによる軌道グラフがリアルタイムで表示され、シミュレーション終了後にエンターキーで終了します。この設定では、デフォルトよりもシミュレーション時間が長くなりますが、より高精度な結果が期待できます。

3. オイラー法を使用し、グラフ表示なしで実行

ソルバーにオイラー法、時間刻み 0.1 日、総ステップ数 20000、グラフ表示なし (fplot=0) で実行します。

python diffeq2nd_planet.py Euler 0.1 20000 0

実行結果:

コンソールにシミュレーションの進行状況のみが表示されます。
diffeq2nd_Planet_Euler.csv と diffeq2nd_Planet_Euler_conservation.csv のファイルは生成されます。
グラフウィンドウは表示されず、バックグラウンドで計算が進みます。シミュレーション終了後、エンターキーでプログラムが終了します。この方法は、大規模な計算やサーバー環境での実行に適しています。