![見出し画像](https://assets.st-note.com/production/uploads/images/95907745/rectangle_large_type_2_3ea89afc5a0a5a03a6e04d0c9003539e.png?width=1200)
Python : 常微分方程式
odeint, BDFの利用について
Stiff な方程式なのか、うまく解けなかったときのメモです。
連続槽型反応器の濃度変化をグラフ化したく、常微分方程式を数値的に解くことを考えました。
問題
反応器サイズ 1L のCSTRに、反応速度式 $${-r_A = k C_A^{1.558}}$$となる原料を 100mmol/Lの濃度で、流量 10, 3, 1.2, 0.5 L/hrの各流量変化させて供給した。各流量切り替え前には、タンク内の反応成分は押し出して、ゼロの状態からスタートする。各流量におけるタンクの濃度変化を確認する。
$$
{ V\frac{dC}{dt} = QC_{A0} - QC_{A} - kC_{A}^n }
$$
以下、計算条件
$${V = 1 [L] \\ Q = (10, 3.0, 1.2, 0.5) [L/hr] \\ k = 0.1397 [(L/mol)^{0.558} / hr] \\ n = 1.558 \\ C_{A0} = 100.0 [mmol/L] }$$
計算
odeint での計算
python scipy にあるodeint でまずは解いてみました。
pythonのodeintに渡す関数の引数の中の並びが標準では、$${(y, t)}$$ で、他のソフトでは、$${ (t, y) }$$ の方が多いかなと思い、tfirst オプションをtrueにして後者になるようにしています。
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy.integrate import odeint
plt.rcParams['figure.figsize'] = (8,8)
plt.rcParams['font.size'] = 22
# モデル式の定義
def rA(t, CA, k,n,Q):
if t < 0.5:
CA0 = 0.0
else:
CA0 = 100
dCA_dt = (Q*CA0 - Q*CA - k * CA**n) / V #
return dCA_dt
V = 1.0 # L
k = 0.1397 # [(L/mol)^0.558/ hr]
n = 1.558 #
CA0 = 100.0 # mmol/L
time = np.linspace(0, 5, 101)
i = 0 # run No.
for Q in np.array([10.0, 3.0, 1.2, 0.5]):
CA = odeint(rA,0, time, args=(k,n,Q), tfirst=True)
i = i + 1
plt.plot(time, CA, label='Run '+ str(i) + ' Q=' + str(Q) + '[L/hr]')
plt.plot([0,0.5,0.5,5],[0,0,100,100], label='step injection')
plt.grid()
plt.xlabel('time[hr]')
plt.ylabel('Concentraion[mmol/L]')
plt.legend()
計算結果のグラフです。流量3.0L/hrの時に計算がうまくいかないようです。
![](https://assets.st-note.com/img/1674011781057-eLpCHP55TM.png)
scipy.integrate.odeint — SciPy v1.10.0 Manual
Scipyのマニュアルとよむと、lsoda法が使われている。
lsoda法は、Adams法とBDF法を切り替えて使うとのこで、Stiff, nonstiffの両方に対応していそうなのですが、うまくいきませんでした。
BDFでの計算
BDF(陰解法)でのトライに変えてみます。
scipy.integrate.BDF — SciPy v1.10.0 Manual
を使ってみました。
以下に変更しています。
from scipy.integrate import BDF
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy.integrate import BDF
plt.rcParams['figure.figsize'] = (8,8)
plt.rcParams['font.size'] = 22
V = 1.0 # L
k = 0.1397 # [(L/mol)^0.558/ hr]
n = 1.558 #
CA0 = 100.0 # mmol/L
t_end = 5.0
i = 0
for Q in np.array([10.0, 3.0, 1.2, 0.5]):
def rA(t, CA):
if t < 0.5:
CA0 = 0.0
else:
CA0 = 100
dCA_dt = (Q*CA0 - Q*CA - k * CA**n) / V #
return dCA_dt
# CA = odeint(rA,0, time)
ODE = BDF(rA, t0=0.0, y0=[0.0], t_bound=t_end)
log_t = []
log_CA = []
while ODE.status == 'running':
ODE.step()
log_t.append(ODE.t)
log_CA.append(ODE.y)
i = i + 1
plt.plot(log_t, log_CA, label='Run '+ str(i) + ' v0=' + str(v0)+'L/hr') #plt .plot(time, CA, linewidth=2, label='Run '+ str(i) + ' v0=' + str(v0))
plt.plot([0,0.5,0.5,5],[0,0,100,100], label='step injection')
plt.grid()
plt.xlabel('time[hr]')
plt.ylabel('Concentraion[mmol/L]')
plt.legend()
![](https://assets.st-note.com/img/1674013050534-4gprp1giWs.png)
うまくいったようです。定常に落ち着く濃度も、問題ないようです。(別途確認しています。)
BDFに渡す常微分方程式内のパラメータを外から渡せないようになるので、for文の中で、反応速度式を毎回定義しなおす形になっています。
まとめ
lsoda法で自動で切り替えてもらうのでは、対応できない事象があることを体験しました。
より一般的な実装になっているode(今回は使っていません)、BDF関数は、使い方が少しわかりにくいですが、慣れれば問題なさそうです。
以上、メモです。
よろしければサポートをお願いします。講習会、有料情報の取得などにあてたいと考えています。やっぱり、単純にうれしいです。