数値計算③ SciPyライブラリ【Pythonと物理#7】

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint 

#数値計算したい微分方程式の定義（オイラー法、ルンゲクッタ法用）
def func(y):
    return y
#数値計算したい微分方程式の定義（scipy.integrate.odeint用）
def func_scipy(y, t):
    return y

#ルンゲ-クッタ,法
def RungeKutta_method(f, y, h):
    k1 = f(y)
    k2 = f(y + k1*h/2)
    k3 = f(y + k2*h/2)
    k4 = f(y + k3*h/2)
    return y + (k1 + 2*k2 + 2*k3 + k4)*h/6

#オイラー法
def euler_method(f, y, h):
    return y + f(y) * h

#常微分方程式を計算
def ode_solve_RKM(f, y0, T, N):
    t = np.linspace(0, T, N)
    h = t[1] - t[0]
    y = [y0, ]
    for i in range(N - 1):
        y.append(RungeKutta_method(f, y[-1], h))
    return t, np.array(y)
    
def ode_solve_EM(f, y0, T, N):
    t = np.linspace(0, T, N)
    h = t[1] - t[0]
    y = [y0, ]
    for i in range(N - 1):
        y.append(euler_method(f, y[-1], h))
    return t, np.array(y)
    
#初期条件などを設定して、上で定義した関数たちを呼び出す
y0 = 0.01
T = 10
N =101
t_scipy = np.linspace(0, T, N) #本来ならtとt_scipyで分けて定義するべきではないと思います

t, y_RKM = ode_solve_RKM(func, y0, T, N) #ルンゲ-クッタ法
t, y_EM = ode_solve_EM(func, y0, T, N) #オイラー法
y_exact = y0 * np.exp(t) #math.exp()だと配列ではないので計算不可
y_scipy = odeint(func_scipy, y0, t_scipy)


#以下は可視化のためにグラフを書くコード
fig, ax = plt.subplots(dpi=240) #解像度上げないと分かりづらかったので、240dpiまで上げています
ax.plot(t_scipy, y_scipy, color='green', label='Scipy Library')
ax.plot(t, y_RKM, color='orange', label='Runge-Kutta Method')
ax.plot(t, y_EM, color='crimson', label='Euler Method')
ax.plot(t, y_exact, "--", color='cornflowerblue', label='Exact')
    
ax.set_xlabel('t')
ax.set_ylabel('y')
ax.legend()
    
ax.tick_params(direction='in', axis='x')
ax.tick_params(direction='in', axis='y')

#以下の範囲設定は適当です
ax.set_xlim(0, 5)
ax.set_ylim(0, 1)

plt.show()