diffeq2nd_planet_euler.py の技術ドキュメント

プログラムの動作

diffeq2nd_planet_euler.py は、複数の天体間の重力相互作用によって引き起こされる運動をシミュレートするPythonプログラムです。このプログラムは、2階常微分方程式をEuler法を用いて数値的に解き、天体の位置の時間発展を追跡します。

主な機能は以下の通りです。

  • CSVファイルから惑星の初期データ（質量、初期位置、初期速度）を読み込みます。

  • 天体間の万有引力に基づいて各天体にかかる合力を計算します。

  • Euler法を用いて、各時間ステップにおける天体の位置と速度を更新します。

  • シミュレーション結果（各天体のx, y座標）をCSVファイルとして出力します。

  • シミュレーションの進行状況を標準出力に表示します。

このプログラムは、天体力学におけるN体問題の数値解を求める基本的な手法を提供し、惑星系のダイナミクスを理解するためのシンプルなモデルとして機能します。

原理

このプログラムは、N個の天体間の重力相互作用による運動を、Euler法という数値積分アルゴリズムを用いて計算します。

万有引力の法則

各天体は、他のすべての天体から万有引力の影響を受けます。質量 \(M_i\) と \(M_j\) を持つ2つの天体間に働く力 \(F_{ij}\) は、天体間の距離 \(r_{ij}\) を用いて以下のように表されます。

\[F_{ij} = G \frac{M_i M_j}{r_{ij}^2}\]

ここで \(G\) は万有引力定数です。プログラムでは、この力を各座標軸の成分に分解して計算します。天体 \(i\) が天体 \(j\) から受ける力 \(F_{ij}\) の \(x\) 成分は、天体間の相対位置 \(\Delta x_{ij}\) を用いて次のように計算されます。

\[F_{ij,x} = F_{ij} \frac{\Delta x_{ij}}{r_{ij}}\]

同様に \(F_{ij,y}\) と \(F_{ij,z}\) も計算されます。各天体にかかる合力は、他のすべての天体からの力のベクトル和として求められます。天体 \(i\) にかかる合力 \(F_i\) が分かれば、その加速度 \(a_i\) はニュートンの運動方程式 \(F_i = M_i a_i\) から得られます。

プログラムでは、単位変換された万有引力定数 G1 を使用しています。 G1 = G * DayToSecond * DayToSecond / AU / AU / AU これは、距離の単位を天文単位 (AU)、時間の単位を日 (Day) に変換するために導入されたものです。

Euler法

運動方程式は2階の常微分方程式であり、これを数値的に解くために最も単純な陽的Euler法が用いられます。Euler法では、現在の状態（位置 \(r\) と速度 \(v\)）と加速度 \(a\) を用いて、次の時間ステップ \(\Delta t\) 後の状態を近似的に予測します。

現在の時刻 \(t\) での位置、速度、加速度をそれぞれ \(r(t)\), \(v(t)\), \(a(t)\) とすると、次の時刻 \(t + \Delta t\) での位置 \(r(t+\Delta t)\) と速度 \(v(t+\Delta t)\) は次のように更新されます。

\[v(t + \Delta t) = v(t) + a(t) \cdot \Delta t\]
\[r(t + \Delta t) = r(t) + v(t) \cdot \Delta t\]

このプログラムでは、各天体について3次元（x, y, z）の位置と速度が計算されます。プログラム内の Fi 関数で現在の加速度 (力 / 質量) を計算し、それを用いて上記の式に従って次のステップの速度と位置を更新しています。

運動量保存の正規化

シミュレーションの開始時に、システム全体の総運動量がゼロになるように各天体の速度が調整されます (normalize_momentum 関数)。これは、システムの重心が初期位置で静止している（または一定速度で移動している）重心系でシミュレーションを行うことを意味し、シミュレーションの安定性を向上させる効果があります。

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

このプログラムは、標準ライブラリ以外に以下のライブラリを使用しています。

  • numpy: 数値計算を効率的に行うためのライブラリです。プログラム内で直接的な配列操作にはあまり使われていませんが、科学技術計算の標準的なツールとしてインポートされています。

インストールは pip コマンドを使用して行えます。

pip install numpy

必要な入力ファイル

プログラムの実行には、惑星の初期設定データが記述されたCSVファイルが必要です。デフォルトのファイル名は planet_db.csv です。

planet_db.csv

  • 形式: カンマ区切り（CSV）

  • 内容: シミュレーション対象となる各惑星の初期状態が記述されています。

  • データ構造: ヘッダー行を含み、各行が1つの惑星のデータに対応します。

    • Name: 惑星の名前（文字列）

    • Mass: 惑星の質量（浮動小数点数、単位はkg）

    • Revolution Radius: 惑星の初期公転半径（浮動小数点数、単位はm）

    • Revolution Velocity: 惑星の初期公転速度（浮動小数点数、単位はm/s）

例:

Name,Mass,Revolution Radius,Revolution Velocity
Sun,1.989e+30,0.0,0.0
Earth,5.972e+24,1.496e+11,29780.0
Mars,6.39e+23,2.279e+11,24130.0
Jupiter,1.898e+27,7.785e+11,13070.0

生成される出力ファイル

シミュレーションの結果は、指定されたCSVファイルに保存されます。デフォルトのファイル名は diffeq2nd_Planet_Euler.csv です。

diffeq2nd_Planet_Euler.csv

  • 形式: カンマ区切り（CSV）

  • 内容: 各時間ステップにおける各惑星のx座標とy座標が記録されます。

  • データ構造:

    • 最初の列はシミュレーション時刻 t (日単位) です。

    • それに続く列は、x(惑星名) と y(惑星名) のペアで、各惑星のx座標とy座標（天文単位 AU）を示します。

    • プログラム内の CSVOutputInterval 定数により、指定された時間ステップごとにデータが出力されます。

例（一部）:

t,x(Sun),y(Sun),x(Earth),y(Earth),x(Mars),y(Mars),x(Jupiter),y(Jupiter)
0.0,0.0,0.0,1.0,0.0,1.5225,0.0,5.2045,0.0
10.0,-0.0000,-0.0000,0.9815,0.1912,1.4792,0.3015,5.1952,0.2223
20.0,-0.0000,-0.0000,0.9262,0.3705,1.3533,0.5843,5.1670,0.4439
...

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

このプログラムは、引数を指定せずに直接Pythonインタープリタで実行します。

python diffeq2nd_planet_euler.py

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

planet_db.csv ファイルがプログラムと同じディレクトリに存在する場合、以下のコマンドで実行できます。

python diffeq2nd_planet_euler.py

実行結果 (標準出力):

実行すると、まずシミュレーションの基本情報（G, AU, G1 の値）が表示され、次に planet_db.csv から読み込まれた各惑星の詳細情報が出力されます。

Planet simulator: Solve simulataneous second order diffrential equations by Euler method
G = 6.67259e-11 Nm2/kg2
AU = 1.495978700000e+11 m
G1 = 3.963212871146665e-14

Planets:
  Sun
     Mass: 1.989e+30
     Revolution Radius: 0.0
     Revolution Velocity: 0.0
  Earth
     Mass: 5.972e+24
     Revolution Radius: 149600000000.0
     Revolution Velocity: 29780.0
  Mars
     Mass: 6.39e+23
     Revolution Radius: 227900000000.0
     Revolution Velocity: 24130.0
  Jupiter
     Mass: 1.898e+27
     Revolution Radius: 778500000000.0
     Revolution Velocity: 13070.0

Write to [diffeq2nd_Planet_Euler.csv]
  t      x(Sun)        y(Sun)        x(Earth)      y(Earth)      x(Mars)        y(Mars)      x(Jupiter)     y(Jupiter)
0.0     0.0000        0.0000        1.0000        0.0000        1.5245        0.0000        5.2066        0.0000
100.0  -0.0000e+00  -0.0000e+00   9.939e-01  -1.109e-01   1.517e+00  -2.193e-01   5.201e+00  -1.378e-01
200.0  -0.0000e+00  -0.0000e+00   9.754e-01  -2.152e-01   1.496e+00  -4.228e-01   5.187e+00  -2.753e-01
...

（上記はダミー出力例です。実際の数値はシミュレーションによって異なります。）

シミュレーションの進行中、iprint_interval で設定された時間間隔（デフォルトでは100ステップごと）で、現在の時刻 t と、nprint_planets で指定された数の惑星（デフォルトでは4つ）のx, y座標が標準出力に表示されます。

生成されるファイル:

実行後、diffeq2nd_Planet_Euler.csv という名前のCSVファイルがプログラムと同じディレクトリに生成されます。このファイルには、シミュレーション開始から終了までの全時間ステップにおける、各惑星のx座標とy座標が記録されています。このデータは、後でグラフ化したり解析したりするのに使用できます。