方程式?
一般的に方程式といえば以下のような式を思い浮かべます。
方程式は xという「数」が未知数であり、xを求めます。
では、微分方程式はどうでしょうか
微分方程式？
微分方程式では未知数が方程式であり、方程式を求めます。
微分方程式の解は一意に定まらなかったり、そもそも解が存在しないことがある。
一意に定まるときは以下の解曲線が無数にあり、初期値などを定めた時に定まる。
微分方程式数値解析手法
オイラー法とルンゲ・クッタ法という手法を紹介します。
オイラー法
数式的意味
微分 $$ \frac{dy}{dt} = \lim_{Δt \to 0} \frac{y(t+Δt)-y(t)}{Δt} $$
微分商 $$ \frac{dy}{dt} = \lim_{Δt \to 0} \frac{y(t+Δt)-y(t)}{Δt} \simeq f(t, y)$$
$$ y(t+Δt) = y(t) + Δt\frac{y’}{1!}+ Δt^{2}\frac{y^{”}}{2!} +Δt^{3}\frac{y^{(3)}}{3!}+Δt^{4}\frac{y^{(4)}}{4!} + ・・・・$$
$$(Δt^{2}\frac{y^{”}}{2} +Δt^{3}\frac{y^{(3)}}{3!}+Δt^{4}\frac{y^{(4)}}{4!} + ・・・・)$$
は、微小なので無視する。
$$ Δt=h$$とすると、
$$ y(t+h) = y(t)+hf(t, y(t))$$
図形的意味
ただし、
$$ \frac{dy}{dt} = \frac{y(t+Δt)-y(t)}{Δt} \simeq f(t, y)$$
として扱ったため、傾きに誤差が生じて、少しずつずれていく（以下pythonで実装時に出現する）
ルンゲ・クッタ法
$$ y(t+Δt) = y(t) + Δt\frac{y’}{1!}+ Δt^{2}\frac{y^{”}}{2!} +Δt^{3}\frac{y^{(3)}}{3!}+Δt^{4}\frac{y^{(4)}}{4!} + ・・・・$$
Δt=hとして、hでまとめると、
$$ y(t+h) = y(t) + h( \frac{y’}{1!}+ h\frac{y^{”}}{2!} +h^{2}\frac{y^{(3)}}{3!}+h^{3}\frac{y^{(4)}}{4!} + ・・・・ )$$
となる。ここで、
$$ k1 = hf(t, y(t))$$
$$ k2 = hf(t+\frac{h}{2}, y(t)+\frac{k_1}{2})$$
$$ k3 = hf(t+\frac{h}{2}, y(t)+\frac{k_2}{2})$$
$$ k4 = hf(t+h, y(t)+k_3)$$
$$ y(t+h) = y(t)+\frac{1}{6}(k_1+2k_2+2k_3+k_4)$$
と表すことができる
計算の中身の話を書くと長くなるので、割愛します。
オイラー法の傾きをより細かく刻んで、部分的に求めて足し合わせていく方法です。
Pythonで実装
オイラー法実装
元の関数を$$y=e^x$$とします。
# ライブラリのインポート
import numpy as np
import matplotlib.pyplot as plt
# f(t,y)
def f(t, y):
return np.exp(t)
# 初期条件
y0 = 1
yt = y0
# 区間
dt = 0.1
t = np.arange(0, 1.1, 0.1)
# オイラー法の変数y
y = np.zeros(len(t))
y[0] = y0
# 実際の関数
y_exact = np.exp(t)
for i in range(len(t) - 1):
y[i+1] = yt + dt *f(t[i], yt)
yt = y[i + 1]
plt.scatter(t, y, color='orange', label="Euler")
plt.plot(t, y_exact, label="function")
plt.legend()
plt.show()
オレンジの点でプロットした部分がオイラー法による解析結果です。
本来の関数とは少しずれがあります。
ずれを無くすにはどうするべきでしょうか？？
答えは、Δtの値を小さくする！です
# 区間
dt = 0.01 # dt = 0.1を0.01にしてみます。
t = np.arange(0, 1.1, 0.1)Δtの値を小さくすることで、区間が細かくなるのでより正確になります。
少しみにくいですが、細かくすることで、より元の関数に近い値をとっていることがわかると思います。
しかし、これは値を大きくしていくと、傾きの誤差がちょっとずつ出てくるので正直意味ない？って感じもします。ww
ルンゲ・クッタ法実装
元の関数を$$y=e^x$$とします。
import numpy as np
import matplotlib.pyplot as plt
# f(t,y)
def f(t, y):
return np.exp(t)
# 初期条件
y0 = 1
yt = y0
dt = 0.1
# 区間
t = np.arange(0, 1.1, dt)
# ルンゲ・クッタ法の変数
y = np.zeros(len(t))
y[0] = y0
# 実際の関数
y_exact = np.exp(t)
for i in range(len(t) - 1):
k1 = dt * f(t[i], y[i])
k2 = dt * f(t[i] + dt/2, y[i] + k1/2)
k3 = dt * f(t[i] + dt/2, y[i] + k2/2)
k4 = dt * f(t[i] + dt, y[i] + k3)
y[i+1] = y[i] + (1/6) * (k1 + 2*k2 + 2*k3 + k4)
plt.scatter(t, y, color='orange', label="Euler")
plt.plot(t, y_exact, label="function")
plt.legend()
ルンゲ・クッタ法では、dt=0.1で実行しましたが、グラフの通り、オイラー法と比較しても、関数に近い値をとっていることがわかります。
以上
コメント