この記事は現代化学(東京化学同人)の2022年4月号から2024年3月号まで連載された安藤耕司先生の「Pythonによる化学シミュレーション入門」を読んでやってみた記事です。詳しい内容は上記の「Pythonによる化学シミュレーション入門」をご参照ください。今回は、2022年6月号の反応速度論(1)、2022年7月号の反応速度論(2)の内容です。
目標
単分子一次反応の解析解と数値解の比較。
・単分子一次反応
A \longrightarrow B
のような反応で、反応速度$v$が次のように表されるもの
v = k A(t) = - \frac{dA(t)}{dt} = \frac{dB(t)}{dt}
ここで$A(t)$、$B(t)$は時刻$t$でのAとBの濃度、$k$は反応速度定数。
単分子一次反応の解析解と数値解
解析解
$k A(t) = - \frac{dA(t)}{dt}$より、この微分方程式を解いて、
A(t) = A(0) exp(-kt)
また、$\frac{dB(t)}{dt} = k A(t) = k A(0) exp(-kt)$より、この微分方程式を解いて
B(t) = B(0) + A(0) - A(0) exp(-kt)
import numpy as np
import matplotlib.pyplot as plt
k, t_max, h = 0.5, 5.0, 0.01
t, a, b = 0.0, 1.0, 0.0
ttn = np.arange(t, t_max, h)
a_exact = a*np.exp(-k*ttn)
b_exact = b + a - a*np.exp(-k*ttn)
apb_exact = a_exact + b_exact
plt.plot(ttn, a_exact, label='[A] exact')
plt.plot(ttn, b_exact, label='[B] exact')
plt.plot(ttn, apb_exact, label='[A]+[B] exact')
plt.legend(fontsize=14)
plt.ylabel("consentration", fontsize=16)
plt.xlabel("time", fontsize=16)
plt.show()
プロットは次のようになります。
数値解1 オイラー法（一次のルンゲ・クッタ法）
$h$が十分小さいとき、$v = - \frac{dA(t)}{dt} \simeq - (\frac{ A(t + h) - A(t)}{h} )$と近似して、
A(t + h) = A(t) - v h
同様に$ v \simeq \frac{ B(t + h) - B(t)}{h} $とするとき
B(t + h) = B(t) + v h
import matplotlib.pyplot as plt
k, t_max, h = 0.5, 5.0, 0.01
t, a, b = 0.0, 1.0, 0.0
tt, aa, bb, cc = [], [], [], []
while t <= t_max:
tt.append(t)
aa.append(a)
bb.append(b)
cc.append(a + b)
t += h
v = k*a
a -= v*h
b += v*h
plt.plot(tt, aa, label='[A]')
plt.plot(tt, bb, label='[B]')
plt.plot(tt, cc, label='[A]+[B]')
plt.legend(fontsize=14)
plt.ylabel("consentration", fontsize=16)
plt.xlabel("time", fontsize=16)
plt.show()
プロットを示します。オイラー法は簡単な近似ですが$h$が十分小さければ解析解との差は小さいです。
数値解2 二次のルンゲ・クッタ法
ここまでは濃度を$A(t)$、$B(t)$としてましたが、より一般化して、次のような変数$t$の関数$x$について考えます。
\displaylines{
x_{n+1} = x_n + f(x_n)*h \\
t_{n+1} = t_n + h
}
このとき、二次のルンゲ・クッタ法では、$s_1$と$s_2$を次のように定め、$x_{n+1}$を近似します。
\displaylines{
s_1 = f(x_n) \\
s_2 = f(x_n + s_1\frac{h}{2}) \\
x_{n+1} \simeq x_n + \frac{s_1 + s_2}{2}*h
}
import numpy as np
import matplotlib.pyplot as plt
k, t_max, h = 0.5, 5.0, 0.01
t, a, b = 0.0, 1.0, 0.0
tt, aa, bb, cc = [], [], [], []
while t <= t_max:
tt.append(t)
aa.append(a)
bb.append(b)
cc.append(a + b)
t += h
s1 = k*a
s2 = k*(a-s1*h/2)
v = (s1 + s2)/2
a -= v*h
b += v*h
plt.plot(tt, aa, label='[A]')
plt.plot(tt, bb, label='[B]')
plt.plot(tt, cc, label='[A]+[B]')
plt.legend(fontsize=14)
plt.ylabel("consentration", fontsize=16)
plt.xlabel("time", fontsize=16)
plt.show()
数値解3 四次のルンゲ・クッタ法
四次のルンゲ・クッタ法についてはこの解説が詳しいです。
\displaylines{
s_1 = f(x_n) \\
s_2 = f(x_n + s_1\frac{h}{2}) \\
s_3 = f(x_n + s_2\frac{h}{2}) \\
s_4 = f(x_n + s_3*h) \\
x_{n+1} \simeq x_n + \frac{s_1 + 2s_2 + 2s_3 + s_4}{6}*h
}
import numpy as np
import matplotlib.pyplot as plt
k, t_max, h = 0.5, 5.0, 0.01
t, a, b = 0.0, 1.0, 0.0
tt, aa, bb, cc = [], [], [], []
while t <= t_max:
tt.append(t)
aa.append(a)
bb.append(b)
cc.append(a + b)
t += h
s1 = k*a
s2 = k*(a-s1*h/2)
s3 = k*(a-s2*h/2)
s4 = k*(a-s3*h)
v = (s1 + 2*s2 + 2*s3 + s4)/6
a -= v*h
b += v*h
plt.plot(tt, aa, label='[A]')
plt.plot(tt, bb, label='[B]')
plt.plot(tt, cc, label='[A]+[B]')
plt.legend(fontsize=14)
plt.ylabel("consentration", fontsize=16)
plt.xlabel("time", fontsize=16)
plt.show()
解析解と数値解の比較
それでは、解析解と数値解の比較を行います。
import numpy as np
import matplotlib.pyplot as plt
def unimolecular_reaction(k, t_max, h, t, a, b):
tt, aa1, bb1, cc1 = [], [], [], []
aa2, bb2, cc2 = [], [], []
aa4, bb4, cc4 = [], [], []
a1, b1 = a, b
a2, b2 = a, b
a4, b4 = a, b
while t <= t_max:
tt.append(t)
aa1.append(a1)
bb1.append(b1)
cc1.append(a1 + b1)
aa2.append(a2)
bb2.append(b2)
cc2.append(a2 + b2)
aa4.append(a4)
bb4.append(b4)
cc4.append(a4 + b4)
t += h
v1 = k*a1
a1 -= v1*h
b1 += v1*h
s2_1 = k*a2
s2_2 = k*(a2-s2_1*h/2)
v2 = (s2_1 + s2_2)/2
a2 -= v2*h
b2 += v2*h
s4_1 = k*a4
s4_2 = k*(a4-s4_1*h/2)
s4_3 = k*(a4-s4_2*h/2)
s4_4 = k*(a4-s4_3*h)
v4 = (s4_1 + 2*s4_2 + 2*s4_3 + s4_4)/6
a4 -= v4*h
b4 += v4*h
ttn = np.array(tt)
a_exact = aa[0]*np.exp(-k*ttn)
b_exact = bb[0] + aa[0] - aa[0]*np.exp(-k*ttn)
apb_exact = a_exact + b_exact
numerical1 = [tt, aa1, bb1, cc1]
numerical2 = [tt, aa2, bb2, cc2]
numerical4 = [tt, aa4, bb4, cc4]
exact = [ttn, a_exact, b_exact, apb_exact]
return numerical1, numerical2, numerical4, exact
numerical1, numerical2, numerical4, exact = unimolecular_reaction(1, 5.0, 0.1, 0, 1, 0)
plt.plot(numerical1[0], numerical1[1], label='[A] Euler')
plt.plot(numerical2[0], numerical2[1], label='[A] RK-2')
plt.plot(numerical4[0], numerical4[1], label='[A] RK-4')
plt.plot(exact[0], exact[1], "--", label='[A] exact')
plt.legend(fontsize=14)
plt.ylabel("consentration", fontsize=16)
plt.xlabel("time", fontsize=16)
plt.yscale("log")
plt.show()
それぞれの$A(t)$を比較しましょう。差がわかりやすいように縦軸を対数スケールで図示します。プログラムでは$h=0.01$ですが、以下に$h=0.01, 0.1, 0.5$の三種類のプロットを示します。
Eulerがオイラー法、RK-2が二次のルンゲ・クッタ法、RK-4が四次のルンゲ・クッタ法、exactが解析解です。$h=0.01$では全ての数値解が解析解とピッタリ一致しているように見えます。$h$が小さいので近似が十分成り立っています。一方、$h=0.5$ではオイラー法と二次のルンゲ・クッタ法は解析解ずれています。$h$が大きいので近似が成り立ちません。しかしながら$h=0.5$においても四次のルンゲ・クッタ法は、解析解と一致しているように見えます。このように、四次のルンゲ・クッタ法はある程度$h$が大きくとも、解析解と近い値が得られることが知られています。