p1_ex07a

常微分方程式

モジュールのインポート

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

import scienceplots
#Warning : As of version 2.0.0, you need to add import scienceplots before setting the style (plt.style.use('science')).
plt.style.use(['science', 'notebook'])

　次のような微分方程式，初期条件に対して， $0 < x \le 2$ の範囲で，数値的に $y$ の値を求めよう． $$ y' = \frac{dy}{dx} = -4(x-1)y $$ ここで，初期条件は， $$ x_0 = 0, \ \ \ \ \ y_0 = e^{-2} $$ 微分方程式を解析的に解くと， $$ y = e^{-2(x-1)^2} $$ まず，この関数は，

def func7(x):
    y = np.exp(-2.0*(x-1.0)*(x-1.0))
    return y

微分方程式は，

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

関数 odeint

Scipy の関数 odeint
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html
より，

x00, y00 = 0.0, np.exp(-2.0) # 初期条件
m = 11
xe = 2.0
xx = np.linspace(x00, xe, m)
ode1 = odeint(dfunc7, y00, xx, tfirst=True)
ode1

array([[0.13533528],
       [0.27803726],
       [0.48675221],
       [0.72614895],
       [0.92311613],
       [0.99999972],
       [0.92311612],
       [0.72614885],
       [0.48675213],
       [0.27803723],
       [0.13533525]])

転置をとって，

ode1.T

array([[0.13533528, 0.27803726, 0.48675221, 0.72614895, 0.92311613,
        0.99999972, 0.92311612, 0.72614885, 0.48675213, 0.27803723,
        0.13533525]])

あるいは，

yy1 = ode1[:,0]
yy1

array([0.13533528, 0.27803726, 0.48675221, 0.72614895, 0.92311613,
       0.99999972, 0.92311612, 0.72614885, 0.48675213, 0.27803723,
       0.13533525])

出力フォーマットを整えて，

print(f"{'x':>8s} {'odeint':>10s} {'error':>10s}")
yy = func7(xx)
for i,x in enumerate(xx):
    print(f"{x:8.3f} {yy1[i]:10.6f} {yy1[i]-yy[i]:10.6f}")

x     odeint      error
0.000   0.135335   0.000000
0.200   0.278037  -0.000000
0.400   0.486752  -0.000000
0.600   0.726149  -0.000000
0.800   0.923116  -0.000000
1.000   1.000000  -0.000000
1.200   0.923116  -0.000000
1.400   0.726149  -0.000000
1.600   0.486752  -0.000000
1.800   0.278037  -0.000000
2.000   0.135335  -0.000000

プロットすると，

xxx = np.linspace(x00, xe, m*10)

fig = plt.figure() # グラフ領域の作成
plt.plot(xxx, func7(xxx))
plt.plot(xx, yy1, "o")
fig.savefig('p1_ex07a_1.pdf')

関数 solve_ivp

Scipy の関数 solve_ivp
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html
より，

x0 = x00
y_ini = np.array([y00])
yy2 = np.empty(0)
print(f"{'x':>8s} {'solve_ivp':>10s} {'error':>10s}")
for i,x in enumerate(xx):
    ode2 = solve_ivp(dfunc7, [x0, x], y_ini,  method='DOP853', dense_output=True)
    ss = ode2.sol(x)
    print(f"{x:8.3f} {ss[0]:10.6f} {ss[0]-yy[i]:10.6f}")
    x0 = x
    y_ini = ss.copy()
    yy2 = np.append(yy2, ss[0])

x  solve_ivp      error
0.000   0.135335   0.000000
0.200   0.278037   0.000000
0.400   0.486752   0.000000
0.600   0.726149   0.000000
0.800   0.923116   0.000000
1.000   1.000000   0.000000
1.200   0.923116   0.000000
1.400   0.726149   0.000000
1.600   0.486752   0.000000
1.800   0.278037   0.000000
2.000   0.135335   0.000000

プロットすると，

xxx = np.linspace(x00, xe, m*10)

fig = plt.figure() # グラフ領域の作成
plt.plot(xxx, func7(xxx))
plt.plot(xx, yy2, "o")
fig.savefig('p1_ex07a_2.pdf')

計算誤差をプロットして比較すると，

fig = plt.figure() # グラフ領域の作成
plt.plot(xx, yy1-yy, "o", label='odeint')
plt.plot(xx, yy2-yy, "o", label='solve_int')
plt.legend(ncol=2, loc='lower left', fancybox=False, frameon = True)
fig.savefig('p1_ex07a_3.pdf')