常微分方程式の数値解法(5)

SciPyのodeintを用いて以下のOrdinary Differential Equations (ODE)を解き，文献[1,2]のChemical Akzo Nobel problemの数値解を再現する．

[1] Amat, S. et al., On a Variational Method for Stiff Differential Equations Arising from CHemistry Kinetics, Mathematics, 7, 459 (2019) [2] Kafash, B. et al., Chemical Akzo Nobel Problem; Mathematical Description and Its Solution Via Model Maker Software, MATCH Commum. Math. Comput. Chem., 71, 279-286 (2014)

文献[1]の再現

#
# Chemical Akzo Nobel Problem in
#
# [1] Amat, S. et al.,
#     On a Variational Method for Stiff Differential Equations
#     Arising from CHemistry Kinetics,
#     Mathematics, 7, 459 (2019)
#
# [2] Kafash, B. et al.,
#     Chemical Akzo Nobel Problem; Mathematical Description
#     and Its Solution Via Model Maker Software,
#     MATCH Commum. Math. Comput. Chem., 71, 279-286 (2014)
#
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# right-hand side
def rhs(y, x, k0, k1, k2, K, k3, klA, Ks, p, H):

    r0 = k0 * y[0]**4 * np.sqrt(y[1])
    r1 = k1 * y[2] * y[3]
    r2 = k1 * y[0] * y[4] / K
    r3 = k2 * y[0] * y[3]**2
    y5 = Ks * y[0] * y[3]
    r4 = k3 * y5**2 * np.sqrt(y[1])
    Fin = klA * (p / H - y[1])

    yp0 = - 2.0 * r0 + r1 - r2 - r3
    yp1 = - 0.5 * r0 - r3 - 0.5e0 * r4 + Fin
    yp2 = r0 - r1 + r2
    yp3 = - r1 + r2 - 2.0 * r3
    yp4 = r1 - r2 + r4
    return [yp0, yp1, yp2, yp3, yp4]

# parameters
k0 = 18.7
k1 = 0.58
k2 = 0.09
K = 34.4
k3 = 0.42
klA = 3.3
Ks = 115.83
p = 0.9
H = 737.0

# initial condition
# 0: FLB, 1: CO2, 2: FLBT, 3: ZHU, 4: ZLA, 5: FLB.ZHU
FLB     = 0.444
CO2     = 0.00123
FLBT    = 0.0
ZHU     = 0.007
ZLA     = 0.0
FLB_ZHU = 0.367
ZHU     = FLB_ZHU/(Ks*FLB)
y0      = [FLB, CO2, FLBT, ZHU, ZLA]

# output interval
x = np.arange(0.0, 200.0, 0.1)

# solve ode
y = odeint(rhs, y0, x, args = (k0, k1, k2, K, k3, klA, Ks, p, H))

# output
plt.xlabel('t')
plt.ylabel('y')
plt.plot(x, y[:,0],           label='FLB')
plt.plot(x, y[:,1]*1e2,       label='CO2 (\u00d7$10^2$)')
plt.plot(x, y[:,2],           label='FLBT')
plt.plot(x, y[:,3]*1e1,       label='ZHU (\u00d7$10^1$)')
plt.plot(x, y[:,4]*1e1,       label='ZLA (\u00d7$10^1$)')
plt.plot(x, Ks*y[:,0]*y[:,3], label='FLB.ZHU')
plt.legend()
#plt.show()
plt.savefig('ChemicalAkzoNobelProblem.png')