Pythonで数値解法を試す - 宇宙電波実験室

先日、「パソコンで見る天体の動き」という本を買いました。名前の通り、天体の軌道をパソコンで表示するための数式やサンプルコードなどが載っています。

内容はとても素晴らしいのですが、サンプルコードが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))

どれもまあ良くできたもんだなと思いますが、刻み幅を勝手に調整してくれるフェールベルグ法なんて特に凄いですよね。初めて知ったときは感動しました。