常微分方程式の数値解法(4) - 流体解析のあれこれ

SciPyのodeintを用いて以下のOrdinary Differential Equations (ODE)を解き，文献[1]のProblem HIRESの数値解を再現する．

$\displaystyle y_1^{\prime} = - 1.7 y_1 + 0.43 y_2 + 8.32 y_3 + 0.0007$

$\displaystyle y_2^{\prime} = 1.7 y_1 - 8.75 y_2$

$\displaystyle y_3^{\prime} = - 10.03 y_3 + 0.43 y_4 + 0.035 y_5$

$\displaystyle y_4^{\prime} = 8.32 y_2 + 1.71 y_3 - 1.12 y_4$

$\displaystyle y_5^{\prime} = - 1.745 y_5 + 0.43 y_6 + 0.43 y_7$

$\displaystyle y_6^{\prime} = - 280 y_6 y_8 + 0.69 y_4 + 1.71 y_5 - 0.43 y_6 + 0.69 y_7$

$\displaystyle y_7^{\prime} = 280 y_6 y_8 - 1.81 y_7$

$\displaystyle y_8^{\prime} = - 280 y_6 y_8 + 1.81 y_7$

[1] Amat, S. et al., On a Variational Method for Stiff Differential Equations Arising from Chemistry Kinetics, Mathematics, 7, 459 (2019)

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

# right-hand side
def rhs(y, x):
        return[- 1.71*y[0] + 0.43*y[1] + 8.32*y[2] + 0.0007,
                 1.71*y[0] - 8.75*y[1],
               - 10.03*y[2] + 0.43*y[3] + 0.035*y[4],
                 8.32*y[1] + 1.71*y[2] - 1.12*y[3],
               - 1.745*y[4] + 0.43*y[5] + 0.43*y[6],
               - 280.0*y[5]*y[7] + 0.69*y[3] + 1.71*y[4] - 0.43*y[5] + 0.69*y[6],
                 280.0*y[5]*y[7] - 1.81*y[6],
               - 280.0*y[5]*y[7] + 1.81*y[6]]

# initial condition
y0 = [1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0057]

# output interval
x = np.arange(0.0,400.0,0.01)

# solve ode
y = odeint(rhs, y0, x)

# outout
figs = plt.figure(figsize=(6.4,4.8*2))
fig1 = figs.add_subplot(4, 2, 1)
fig2 = figs.add_subplot(4, 2, 2)
fig3 = figs.add_subplot(4, 2, 3)
fig4 = figs.add_subplot(4, 2, 4)
fig5 = figs.add_subplot(4, 2, 5)
fig6 = figs.add_subplot(4, 2, 6)
fig7 = figs.add_subplot(4, 2, 7)
fig8 = figs.add_subplot(4, 2, 8)
fig1.set_xlim(0.0,5.0)
fig2.set_xlim(0.0,5.0)
fig3.set_xlim(0.0,5.0)
fig4.set_xlim(0.0,5.0)
fig5.set_xlim(0.0,400.0)
fig6.set_xlim(0.0,400.0)
fig7.set_xlim(0.0,5.0)
fig8.set_xlim(0.0,5.0)
fig1.set_ylabel('y1')
fig2.set_ylabel('y2')
fig3.set_ylabel('y3')
fig4.set_ylabel('y5')
fig5.set_ylabel('y5')
fig6.set_ylabel('y6')
fig7.set_ylabel('y7')
fig8.set_ylabel('y8')
fig7.set_xlabel('t')
fig8.set_xlabel('t')
fig1.plot(x, y[:,0])
fig2.plot(x, y[:,1])
fig3.plot(x, y[:,2])
fig4.plot(x, y[:,3])
fig5.plot(x, y[:,4])
fig6.plot(x, y[:,5])
fig7.plot(x, y[:,6])
fig8.plot(x, y[:,7])
figs.tight_layout()
#plt.show()
figs.savefig('tmp.png')