化学反応でルンゲクッタ

#ルンゲクッタ法
オイラー法に比べて良い近似解を与える(?)数値解法
やっていることはオイラー法の発展形なので恐れる必要はない

##ルンゲクッタ法の解法
次の漸化式を解く
$h$は刻み幅、$ \ f(t, y)=\frac{dy}{dt} $である

y_{n+1}=y_{n}+\frac{h}{6}(k_{1}+2k_{2}+2k_{3}+k_{4}) \\

ただし、

\begin{align}
k_{1}&=f(t_{n},y_{n})=-ky_{n} \\
k_{2}&=f(t_{n}+\frac{dt}{2},y_{n}+\frac{dt}{2}k_{1})=-k(y_{n}+\frac{dt}{2}k_{1})　\\
k_{3}&=f(t_{n}+\frac{dt}{2},y_{n}+\frac{dt}{2}k_{2})=-k(y_{n}+\frac{dt}{2}k_{2}) \\
k_{4}&=f(t_{n}+dt,y_{n}+dtk_{3})=-k(y_{n}+k_{3})
\end{align}

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

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

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

\begin{align}
k_{1}&=f(t_{n},y_{n})=-k[A]_{n} \\
k_{2}&=f(t_{n}+\frac{dt}{2},y_{n}+\frac{dt}{2}k_{1})=-k([A]_{n}+\frac{dt}{2}k_{1})　\\
k_{3}&=f(t_{n}+\frac{dt}{2},y_{n}+\frac{dt}{2}k_{2})=-k([A]_{n}+\frac{dt}{2}k_{2}) \\
k_{4}&=f(t_{n}+dt,y_{n}+dtk_{3})=-k([A]_{n}+k_{3})
\end{align}

pythonコードはこうなる

runge.py

A=1
log=[A]

for _ in range(step):
    k1=-k*A
    k2=-k*(A+k1/2*h)
    k3=-k*(A+k2/2*h)
    k4=-k*(A+k3*h)
    A=A+h/6*(k1+2*k2+2*k3+k4)
    log.append(A)

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

main.py

import math
import matplotlib.pyplot as plt

#初期値・定数記述
A=1
AA=1
AAA=1

h=3
k=0.5
step=10

log=[A]
log2=[AA]
log3=[AAA]

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

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

    k1=-k*AAA
    k2=-k*(AAA+k1/2*h)
    k3=-k*(AAA+k2/2*h)
    k4=-k*(AAA+k3*h)
    AAA=AAA+h/6*(k1+2*k2+2*k3+k4)
    log3.append(AAA)

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.plot([i for i in range(step+1)],log3,label="runge")
plt.legend()

オイラー法は誤差がでかくなりすぎて発散しかけてる
修正オイラー法とルンゲクッタ法の性能差を知りたいならもっと複雑な式の方がよさそう

#まとめ
前回のオイラー法に続いてルンゲクッタ法をやった
よく分からないけど陽的に解くにはこれが最強?