常微分方程式を解く

はじめに

記事の対象者

• Pythonの記法がわかる人
• 常微分方程式を解きたい人

実行環境

• Python 3.10.9
• Matplotlib 3.6.2

原理

いろんな人が書いているので省略．

Code

def func(x, y):
    return x * x - y

def func2(x):
    return x ** 2 - 2 * x + 2

def Euler(h, f, x, y):
    return x + h, y + h * f(x, y)

def Heun(h, f, x, y):
    _x = x + h
    k1 = h * f(x, y)
    k2 = h * f(x + h, y + k1)
    _y = y + 0.5 * (k1 + k2)
    return _x, _y

def RungeKutta(h, f, x, y):
    _x = x + h
    k1 = h * f(x, y)
    k2 = h * f(x + h / 2.0, y + k1 / 2.0)
    k3 = h * f(x + h / 2.0, y + k2 / 2.0)
    k4 = h * f(x + h, y + k3)
    _y = y + (k1 + 2 * (k2 + k3) + k4) / 6.0
    return _x, _y

def calc(method, name):
    h = 0.125
    x0 = 0
    y0 = 2
    xi, yi = method(h, func, x0, y0)
    print(name, "method")
    for i in range(1, 8+1):
        print("{} : {:0.6f} {:0.6f}, err = {:0.10f}".format(i, xi, yi, func2(xi) - yi))
        xi, yi = method(h, func, xi, yi)
    print()

if __name__ == "__main__":
    calc(Euler, "Euler")
    calc(Heun, "Heun")
    calc(RungeKutta, "Runge-Kutta")