化学反応でオイラー法

#数値計算やりたい
数学弱いけど数値計算やってみたい
数値計算は暗黙知的なノウハウが研究室に脈々と伝承され続けてるって聞くけど
極めるには独学じゃ無理かな

#化学反応のシミュレーション
物質Aが一次反応で消失していく場合を考える

\frac{d[A]}{dt} = -k[A]

##オイラー法
オイラー法は次の漸化式を解く
$h$は刻み幅、$ \ f(t, y)=\frac{dy}{dt} $である

y_{n+1}=y_{n}+h*f(t,y)

上記の化学反応に当てはめると、次のように表せる

\begin{align}
[A]_{t+1}&=[A]_{t}+h(-k[A]_{t}) \\
&=[A]_t(1-kh)
\end{align}

pythonコードはこうなる

euler.py

log=[A]
for i in range(step):
    A=A*(1-k*h)
    log.append(A)

##修正オイラー法
修正オイラー法は次の漸化式を解く
$h$は刻み幅、$ \ f(t, y)=\frac{dy}{dt} $である

y_{n+1}=y_{n}+\frac{h}{2}(f(t,y)+f(t+1,y+1))

上記の化学反応に当てはめると、次のように表せる

[A]_{t+1}=[A]_{t}+\frac{h}{2}(-k[A]_{t}-k[A]_{t+1})

(1+\frac{h}{2}k)[A]_{t+1}=(1-\frac{h}{2}k)[A]_{t}

[A]_{t+1}=\frac{(1-\frac{h}{2}k)}{(1+\frac{h}{2}k)}[A]_{t}

pythonコードはこうなる

modified_euler.py

log2=[AA]
for _ in range(step):
    AA=AA*(1-k*h/2)/(1+k*h/2)
    log2.append(AA)

#結果
以下のコードを実行

main.py

import matplotlib.pyplot as plt

#初期値・定数記述
A=1
AA=1
h=1
k=0.3
step=20

log=[A]
log2=[AA]

for i in range(step):
    A=A*(1-k*h)
    log.append(A)

    AA=AA*(1-k*h/2)/(1+k*h/2)
    log2.append(AA)

plt.plot([i for i in range(step+1)],[math.e**(-k*t*h) for t in range(step+1)],label="Exact_solution")
plt.plot([i for i in range(step+1)],log,label="Euler")
plt.plot([i for i in range(step+1)],log2,label="Modified-Euler")
plt.legend()

オイラー法は若干下方向にずれているのに対して、修正オイラー法は厳密解とほぼ重なっている

#まとめ
修正オイラー法の威力が実感できてよかった
式変形あってるか自信ない
あとでルンゲクッタとかもやるかも