p2_ex07

常微分方程式

モジュールのインポート

import numpy as np
from scipy.integrate import solve_ivp
# from scipy.integrate import odeint

オイラー法

# x0, y0: 初期条件，y0 = y'(x0)
# xs: yを求めるxの値
# dfunc: 導関数
# n: 区間数
# y: 数値解，y = y(xs)
def euler(dfunc, x0, y0, xs, n): # オイラー法
    h = (xs-x0)/n
    y = y0
    for i in range(n):
        x = x0+h*i
        y = y+h*dfunc(x,y)
    return y

修正オイラー法（2次のルンゲ・クッタ法）

# x0, y0: 初期条件，y0 = y'(x0)
# xs: yを求めるxの値
# dfunc: 導関数
# n: 区間数
# y: 数値解，y = y(xs)
def rk2(dfunc, x0, y0, xs, n): # 修正オイラー法（2次のルンゲ・クッタ法）
    h = (xs-x0)/n
    y = y0
    for i in range(n):
       x = x0+h*i
       k1 = h*dfunc(x,y)
       k2 = h*dfunc(x+h,y+k1)
       y = y+(k1+k2)*0.5
    return y

4次のルンゲ・クッタ法

# x0, y0: 初期条件，y0 = y'(x0)
# xs: yを求めるxの値
# dfunc: 導関数
# n: 区間数
# y: 数値解，y = y(xs)
def rk4(dfunc, x0, y0, xs, n): # 4次のルンゲ・クッタ法
    h = (xs-x0)/n
    h2 = h*0.5
    y = y0
    for i in range(n):
       x = x0+h*i
       k1 = h*dfunc(x,y)
       k2 = h*dfunc(x+h2,y+k1*0.5)
       k3 = h*dfunc(x+h2,y+k2*0.5)
       k4 = h*dfunc(x+h,y+k3)
       y = y+(k1+2.0*(k2+k3)+k4)/6.0
    return y

Scipy の常微分方程式の関数の一つは，
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html
まず，次のような関数を考える．

def dfunc7(x,y): # 微分方程式 y' = dy/dx = f(x,y)
    dydx = -4.0*(x-1.0)*y
    return dydx

これは解析的に解け，次のような関数で表される．

def func7(x): # 微分方程式を解析的に解いて得られた関数
    y = np.exp(-2.0*(x-1.0)*(x-1.0))
    return y

これらより常微分方程式を数値的に解いてみよう．

eps = 1.0e-9
n, m = 5, 11
x00, y00 = 0.0, np.exp(-2.0)# 初期条件 for func7
xe = 2.0
dx = (xe-x00)/(m-1)
x0 = x00
y0a, y0b, y0c = y00, y00, y00
y_ini = np.array([y00])
print(f"{'x':>8s} {'solve_ivp':>10s} {'er1':>10s} {'ek2':>10s}"\
      f"{'er2':>10s} {'rk4':>10s} {'er4':>10s}")
for i in range(m):
    x = x00+dx*i
    y = func7(x)
    y1 = euler(dfunc7, x0, y0a, x, n)
    y2 = rk2(dfunc7, x0, y0b, x, n)
    y4 = rk4(dfunc7, x0, y0c, x, n)
    solver = solve_ivp(dfunc7, [x0, x], y_ini, dense_output=True)
    y5 = solver.sol(x)
    print(f"{x:8.3f} {y5[0]:10.6f} {y5[0]-y:10.6f} {y2:10.6f}"\
          f"{y2-y:10.6f} {y4:10.6f} {y4-y:10.6f}")
    x0 = x
    y0a, y0b, y0c = y1, y2, y4
    y_ini = y5.copy()

x  solve_ivp        er1        ek2       er2        rk4        er4
0.000   0.135335   0.000000   0.135335  0.000000   0.135335   0.000000
0.200   0.278038   0.000000   0.277404 -0.000633   0.278037  -0.000001
0.400   0.486753   0.000001   0.485100 -0.001652   0.486751  -0.000002
0.600   0.726150   0.000001   0.723366 -0.002783   0.726147  -0.000002
0.800   0.923118   0.000001   0.919466 -0.003651   0.923113  -0.000003
1.000   1.000002   0.000002   0.996014 -0.003986   0.999997  -0.000003
1.200   0.923118   0.000001   0.919419 -0.003697   0.923113  -0.000003
1.400   0.726150   0.000001   0.723298 -0.002851   0.726147  -0.000002
1.600   0.486753   0.000001   0.485056 -0.001696   0.486751  -0.000002
1.800   0.278038   0.000001   0.277421 -0.000616   0.278037  -0.000001
2.000   0.135337   0.000001   0.135411  0.000076   0.135335   0.000000