Python / 古典的 Runge-Kutta 法 / 1 ステップだけ計算してみる

def rk(func, t0, init, h):
    f1 = func(t0        , init             )
    f2 = func(t0 + h / 2, init + h * f1 / 2)
    f3 = func(t0 + h / 2, init + h * f2 / 2)
    f4 = func(t0 + h    , init + h * f3    )
    return init + (h / 6) * (f1 + 2 * (f2 + f3) + f4)

# 次の微分方程式を 1/8 刻みで 1 ステップだけ解いてみる。真値は 0.236111111......
# x'(t) = x(t)^2 - t^2 - 2t + 2
# x(0) = 0
def xDot(t, x):
    return x**2 - t**2 - 2 * t + 2

print(rk(xDot, 0, 0, 1/8))

f:id:ti-nspire:20180105064105p:plain:w300
f:id:ti-nspire:20180105064748p:plain:w500
f:id:ti-nspire:20180105094712p:plain:w500
f:id:ti-nspire:20180105115822p:plain:w500