先日、「パソコンで見る天体の動き」という本を買いました。名前の通り、天体の軌道をパソコンで表示するための数式やサンプルコードなどが載っています。
内容はとても素晴らしいのですが、サンプルコードがBASICやFortranで書かれていたので、理解を深めるためにいくつかの解法をPythonで書き直してみました。
数値解法の詳しい説明を知りたい方は、是非「パソコンで見る天体の動き」をご覧ください。とても面白いですよ。
以下の微分方程式をオイラー法とルンゲ・クッタ法、フェールベルグ法を使って解いてみます。
オイラー法(Euler method)で求める
# f(t,x) = x**2 - t**2 - 2*t + 2 # t = 0 のとき、x = 0 であるものとする。 # 初期値 t = 0 x = 0 # 刻み幅 h = 0.125 f = lambda t, x : x**2 - t**2 2*t + 2 for i in range(10): # 「刻み幅 × 関数の傾き」で増加量が求まる。 x += h * f(t, x) t += h print("t:{0}, x:{1}".format(t, x))
ルンゲ・クッタ法(Runge-Kutta method)で求める
# f(t,x) = x**2 - t**2 - 2*t + 2 # t = 0 のとき、x = 0 であるものとする。 # 初期値 t = 0 x = 0 # 刻み幅 h = 0.125 f = lambda t, x: x**2 - t**2 - 2*t + 2 for i in range(10): f1 = f(t, x) f2 = f(t + h/2, x + h * f1/2) f3 = f(t + h/2, x + h * f2/2) f4 = f(t + h, x + h * f3) x += (h/6)*(f1 + 2*f2 + 2*f3 + f4) t += h print("t:{0}, x:{1}".format(t, x))
フェールベルグ法(Fehlberg method)で求める
# f(t,x) = x**2 - t**2 - 2*t + 2 # t = 0 のとき、x = 0 であるものとする。 # 初期値 t = 0 x = 0 # 刻み幅 h = 0.125 # 誤差条件 alfa_e = 0.0625 # 誤差がこの値に収まるように刻み幅を自動で調整する。 f = lambda t, x: x**2 - t**2 - 2*t + 2 for i in range(10): f1 = f(t, x) t1 = t + h/4 x1 = x + (h/4)*f1 f2 = f(t1, x1) t2 = t + 3*h/8 x2 = x + (h/32) * (3*f1 + 9*f2) f3 = f(t2, x2) t3 = t + 12*h/13 x3 = x + (h/2197) * (1932*f1 - 7200*f2 + 7296*f3) f4 = f(t3, x3) t4 = t + h x4 = x + h*(439*f1/216 - 8*f2 + 3680*f3/513 - 845*f4/4104) f5 = f(t4, x4) t5 = t + h/2 x5 = x + h*(-8*f1/27 + 2*f2 - 3544*f3/2565 + 1859*f4/4104 - 11*f5/40) f6 = f(t5, x5) # 以上より,1ステップ進んだ関数値が求まる.(5次の解) x += h*(16*f1/135 + 6656*f3/12825 + 28561*f4/56430 - 9*f5/50 + 2*f6/55) t += h # 誤差計算用に4次の精度を持つ解の近似値を求める. x_4 = x + h*(25*f1/216 + 1408*f3/2565 + 2197*f4/4104 - f5/5) # 誤差条件を満たさなくなったとき,新しい刻み幅を決める. k = 0.9 # 安全係数 beta_e = abs(x_4 - x) if beta_e > alfa_e: h = k*h * (alfa_e / beta_e)**(1/4) print("調整") print("t:{0}, h:{1}, x:{2}".format(t, h, x))
どれもまあ良くできたもんだなと思いますが、刻み幅を勝手に調整してくれるフェールベルグ法なんて特に凄いですよね。初めて知ったときは感動しました。
精度比較
因みに、今回使った微分方程式では以下の式で真値が求まります。
以上3つの解法によって求めた数値と真値を比較すると、以下のようになります。

色の割り当ては、
- オイラー法:青
- ルンゲ・クッタ法:赤
- フェールベルグ法:黄
- 真値:緑
としました。
オイラー法は t=0.75 付近から誤差が増え始め、t=1.50以降では爆発的に増大していることがわかります。
ルンゲ・クッタ法とフェールベルグ法に関しては、精度が高いので真値と重なって見えなくなっています。試しに少しずらして表示してみると以下のようになり、確かにプロットされていることが確認できます。
おわりに
フェールベルグ法は確かに精度が高くなりますが、刻み幅を調節すれば実用上はルンゲ・クッタ法で十分なように思います。
参考文献
- 長沢工,桧山澄子:「パソコンで見る天体の動き」,株式会社 地人書館(1992年9月20日)