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

プログラムの動作

diffeq2nd_heun.py は、2階の常微分方程式をヘウン法（Heun's method）を用いて数値的に解くPythonプログラムです。

このプログラムの主な機能は以下の通りです。

与えられた初期条件 ($x_0, v_0$) と時間ステップ幅 ($\Delta t$) を用いて、2階常微分方程式 $d^2x/dt^2 = \text{force}(t, x)$ の数値解を計算します。
特定の解析解（例: $x(t) = \sin(t)$）と比較できるよう、解析解も同時に計算します。
計算された時間 $t$、位置 $x$ (計算値)、位置 $x$ (解析解)、速度 $v$ (計算値) を、コンソールに出力するとともに、CSVファイルに保存します。

このプログラムは、解析的な解法が困難または存在しない2階常微分方程式の近似解を求める際に利用できます。

原理

このプログラムは、以下の2階常微分方程式をヘウン法で数値的に解きます。

\[\frac{d^2x}{dt^2} = F(t, x)\]

ここで、$F(t, x)$ はプログラム中の force(t, x) 関数によって定義されており、デフォルトでは $F(t, x) = -x$ です。 2階常微分方程式を解くために、まずこれを1階の連立常微分方程式に変換します。

\[ \frac{dx}{dt} = v \]

\[ \frac{dv}{dt} = F(t, x) \]

ヘウン法は、予測子・修正子法の一種で、予測ステップ（オイラー法）で次の点の傾きを予測し、その予測値を用いて現在の点と予測点の傾きの平均をとって修正するというものです。一般的なヘウン法（予測子・修正子法）のステップは以下の通りです。

ある時刻 $t_i$ における $x_i, v_i$ から、次の時刻 $t_{i+1} = t_i + \Delta t$ における $x_{i+1}, v_{i+1}$ を計算する場合：

予測子 (Predictor): オイラー法で次のステップの値を予測します。 $$ x^*_{i+1} = x_i + \Delta t \cdot v_i $$ v^*_{i+1} = v_i + \Delta t \cdot F(t_i, x_i) $$
修正子 (Corrector): 現在の点と予測点での傾きの平均を用いて、最終的な値を修正します。 $$ x_{i+1} = x_i + \frac{\Delta t}{2} (v_i + v^*_{i+1}) $$ v_{i+1} = v_i + \frac{\Delta t}{2} (F(t_i, x_i) + F(t_i+\Delta t, x^*_{i+1})) $$

プログラム diffeq2nd_heun.py におけるヘウン法の実装は、上記の一般的な形式とは異なる形で記述されています。コード内の変数 x0, v0, t0 はそれぞれ $x_i, v_i, t_i$ に対応し、x1, v1 が $x_{i+1}, v_{i+1}$ に対応します。

プログラムの実装では、以下のステップで数値積分を行っています。

速度 $v$ の更新: まず、force 関数を用いた2つの「傾き」関連の項を計算します。 $$ k_{v,0} = \Delta t \cdot F(t_0, x_0) $$ k_{v,1} = \Delta t \cdot F(t_0 + \Delta t, x_0 + k_{v,0}) $$ ここで、$x_0 + k_{v,0}$ の部分は、$x_0 + \Delta t \cdot F(t_0, x_0)$ を意味しており、一般的なヘウン法における $x$ の予測値 $x_0 + \Delta t \cdot v_0$ とは異なります。この式では、$x_0$ に速度の変化量（の推定値）を加えたものを $x$ の値として force 関数に渡しています。

これらの項を用いて、新しい速度 $v_1$ を計算します。 $$ v_1 = v_0 + \frac{1}{2}(k_{v,0} + k_{v,1}) $$
位置 $x$ の更新: 次に、速度 $v$ を用いた2つの「傾き」関連の項を計算します。 $$ k_{x,0} = \Delta t \cdot v_0 $$ k_{x,1} = \Delta t \cdot v_1 $$ ここで、$k_{x,1}$ の計算には、上記のステップで既に修正された速度 $v_1$ を使用しています。

これらの項を用いて、新しい位置 $x_1$ を計算します。 $$ x_1 = x_0 + \frac{1}{2}(k_{x,0} + k_{x,1}) $$

上記の計算を時間ステップ数 nt だけ繰り返すことで、指定された時間範囲における解を求めます。

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

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

numpy: 高度な数値計算を効率的に行うためのライブラリです。

numpy は pip コマンドでインストールできます。

pip install numpy

必要な入力ファイル

このプログラムは、実行時に特定の入力ファイルを必要としません。初期条件 ($x_0, v_0, \Delta t, nt, \text{iprint\_interval}$) はプログラムのソースコード内で直接定義されており、dt, nt, iprint_interval はコマンドライン引数で上書き可能です。

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

プログラムは実行されると、以下のCSVファイルを生成します。

diffeq2nd_heun.csv: 計算された各時刻 $t$ における数値解と解析解、および速度が出力されます。ファイルはヘッダー行を含み、カンマ区切りで以下のデータが格納されます。
- t: 時刻
- x(cal): 計算された位置 $x$ の値
- x(exact): 解析的に求められた位置 $x$ の値
- v(cal): 計算された速度 $v$ の値
例:
```
t,x(cal),x(exact),v(cal)
0.0,-0.000000,-0.000000,1.000000
0.01,0.010000,0.010000,1.000000
0.02,0.020000,0.020000,0.999900
...
```

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

基本的な実行コマンドと引数の説明は以下の通りです。

python diffeq2nd_heun.py [dt] [nt] [iprint_interval]

dt (オプション): 時間ステップ幅。浮動小数点数で指定します。デフォルト値は 0.01 です。
nt (オプション): シミュレーションの総時間ステップ数。整数で指定します。デフォルト値は 501 です。
iprint_interval (オプション): 計算結果をコンソールに出力する間隔。整数で指定します。デフォルト値は 20 です（1ステップ目と、20ステップごとに表示）。

引数を指定しない場合、プログラムはデフォルト値を使用します。引数を複数指定する場合、左から順に dt, nt, iprint_interval の順で解釈されます。例えば、dt のみを変更したい場合は python diffeq2nd_heun.py 0.005 のように指定し、nt も変更したい場合は python diffeq2nd_heun.py 0.005 1000 のように指定します。

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

例1: デフォルト引数での実行

時間ステップ幅 $\Delta t=0.01$、総ステップ数 $N_t=501$、コンソール出力間隔 $20$ ステップのデフォルト設定でプログラムを実行します。

python diffeq2nd_heun.py

実行結果（コンソール出力の一部）：

Solve second order diffrential equation by Heun method
Write to [diffeq2nd_heun.csv]
  t        x(cal)      x(exact)      v(cal)   
t= 0.00    -0.000000    -0.000000     1.000000
t= 0.20     0.198668     0.198669     0.980067
t= 0.40     0.389417     0.389418     0.921061
t= 0.60     0.564639     0.564642     0.825336
t= 0.80     0.717355     0.717356     0.696706
t= 1.00     0.841470     0.841471     0.540302
...

同時に diffeq2nd_heun.csv が生成され、全ての計算結果が保存されます。

例2: 引数を指定して実行

時間ステップ幅 $\Delta t=0.005$、総ステップ数 $N_t=1001$、コンソール出力間隔 $50$ ステップでプログラムを実行します。

python diffeq2nd_heun.py 0.005 1001 50

実行結果（コンソール出力の一部）：

Solve second order diffrential equation by Heun method
Write to [diffeq2nd_heun.csv]
  t        x(cal)      x(exact)      v(cal)   
t= 0.00    -0.000000    -0.000000     1.000000
t= 0.25     0.247404     0.247404     0.968911
t= 0.50     0.479426     0.479426     0.877583
t= 0.75     0.681639     0.681639     0.731688
t= 1.00     0.841471     0.841471     0.540302
t= 1.25     0.951600     0.951601     0.315322
...

この実行により、より細かいステップで計算が行われ、diffeq2nd_heun.csv も更新されます。