p1_ex07e

連立常微分方程式

モジュールのインポート

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 1/2$の範囲で，数値的に$y$，$z$の値を計算しよう． $$ y' = \frac{dy}{dx} = 3y+4z, \ \ \ \ \ z' = \frac{dz}{dx} = 2y+5z $$ ここで，初期条件は， $$ x_0 = 0, \ \ \ \ \ y_0 = 3, \ \ \ \ \ z_0 = 0 $$ 微分方程式を解析的に解くと[数学演習，p.272 [11]（4）]， $$ y = 2c_1 e^x + c_2 e^{7x}, \ \ \ \ \ z = -c_1 e^x + c_2 e^{7x} $$ ただし，$x=0$のとき， $$ c_1 = \frac{y(0)-z(0)}{3}, \ \ \ \ \ c_2 = \frac{y(0)+2z(0)}{3} $$ より，$y(0)=3$，$z(0)=0$とおくと，$c_1 =1$，$c_2=1$となるので， $$ y = 2 e^x + e^{7x}, \ \ \ \ \ z = -e^x + e^{7x} $$ まず，この関数は，

def func7(x):
    return [2*np.exp(x)+np.exp(7*x), -np.exp(x)+np.exp(7*x)]

微分方程式は，

def dfunc7(x, s): # Coupled Second order ODEs
    y, z = s
    return [3*y+4*z, 2*y+5*z]

関数 odeint

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

x00, y00, z00 = 0.0, 3.0, 0.0 # 初期条件
s00 = (y00, z00)
m = 11
xe = 0.5
xx = np.linspace(x00, xe, m)
ode1 = odeint(dfunc7, s00, xx,  tfirst=True)
yy1, zz1 = ode1.T[0], ode1.T[1]

解析的に得られた式から誤差を求めると，

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

x   odeint y      error   odeint z      error
0.000   3.000000   0.000000   0.000000   0.000000
0.050   3.521610  -0.000000   0.367796  -0.000000
0.100   4.224095   0.000000   0.908582   0.000000
0.150   5.181320   0.000000   1.695817   0.000000
0.200   6.498006   0.000000   2.833797   0.000000
0.250   8.322654   0.000000   4.470577   0.000000
0.300  10.865888   0.000000   6.816311   0.000000
0.350  14.426482   0.000000  10.169279   0.000000
0.400  19.428297   0.000001  14.952823   0.000001
0.450  26.472690   0.000001  21.767753   0.000001
0.500  36.412896   0.000001  31.466732   0.000001

プロットすると，

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

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

関数 solve_ivp

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

x0 = x00
y_ini = s00
yy2 = np.empty(0)
zz2 = np.empty(0)
print(f"{'x':>8s} {'solve_ivp y':>12s} {'error':>12s} {'solve_ivp z':>12s} {'error':>12s}")
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]:12.8f} {ss[0]-yy[i]:12.8f} {ss[1]:12.8f} {ss[1]-zz[i]:12.8f}")
    x0 = x
    y_ini = ss.copy()
    yy2 = np.append(yy2, ss[0])
    zz2 = np.append(zz2, ss[1])

x  solve_ivp y        error  solve_ivp z        error
0.000   3.00000000   0.00000000   0.00000000   0.00000000
0.050   3.52160974  -0.00000000   0.36779645  -0.00000000
0.100   4.22409454  -0.00000000   0.90858179  -0.00000000
0.150   5.18131960  -0.00000000   1.69581688  -0.00000000
0.200   6.49800548  -0.00000000   2.83379721  -0.00000000
0.250   8.32265351  -0.00000000   4.47057726  -0.00000000
0.300  10.86588753  -0.00000000   6.81631110  -0.00000000
0.350  14.42648182  -0.00000000  10.16927917  -0.00000000
0.400  19.42829617  -0.00000000  14.95282207  -0.00000000
0.450  26.47268895  -0.00000000  21.76775239  -0.00000000
0.500  36.41289450  -0.00000000  31.46673069  -0.00000000

プロットすると，

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

yy1,yy2,yy

(array([ 3.        ,  3.52160974,  4.22409459,  5.18131964,  6.49800556,
         8.32265363, 10.86588772, 14.42648211, 19.42829669, 26.47268983,
        36.41289594]),
 array([ 3.        ,  3.52160974,  4.22409454,  5.1813196 ,  6.49800548,
         8.32265351, 10.86588753, 14.42648182, 19.42829617, 26.47268895,
        36.4128945 ]),
 array([ 3.        ,  3.52160974,  4.22409454,  5.1813196 ,  6.49800548,
         8.32265351, 10.86588753, 14.42648182, 19.42829617, 26.47268895,
        36.4128945 ]))

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

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

fig = plt.figure() # グラフ領域の作成
plt.plot(xx, zz1-zz, "o", label='odeint z')
plt.plot(xx, zz2-zz, "o", label='solve_int z')
plt.legend(ncol=2, loc='upper right', fancybox=False, frameon = True)
fig.savefig('p1_ex07e_4.pdf')