※Amazonのアソシエイトとして、当メディア(Nラボ備忘録)は適格販売により収入を得ています。
今回のターゲット：オイラー法（Euler Method）
前回は減衰振動の微分方程式を計算してみました。
今回から複数回に分けて、常微分方程式の数値計算法を勉強してみようかと思います。
今回は最も基本的な公式である オイラー法 を勉強してみました。
勉強の際に参考にさせていただいた書籍もご紹介させていただきます。具体的なPythonのコードも載っており、大変勉強になる１冊です。
オイラー法（Euler Method）
オイラー法は以下の式で定義されます。
$$
y_{n+1} = y_{n} + hf\left(t_{n}, y_{n}\right)
$$
ここで、$y_{n}$ は $t_{n}$ の時の解であり、$t_{n} = t_{0} + nh$ です。
このままの状態だとよくわかりませんが、以下のように書き換えるとわかりやすいかもしれません。
$$
y\left(t_{n} + h\right) = y\left(t_{n}\right) + h\cfrac{dy\left(t_{n}\right)}{dt}
$$
上の式は $y\left(t_{n} + h\right)$ を $t = t_{n}$ の周りでTaylor展開して、$h^{2}$ 以上の項を無視すると得られます。このように考えれば、少しわかりやすく感じました。
オイラー法で数値計算をしてみる
それではオイラー法を使って数値計算をしてみます。今回は、下の式をオイラー法で数値計算してみます。
$$
\cfrac{dy}{dt} = y
$$
これは普通に厳密解 $y$ を求めることもできるので、求めてみると以下のようになります。
$$
y = y_{0}\exp\left(t\right)
$$
ただし、 $y_{0}$ は定数です。
Pythonで計算する
コードは下のようにしました。詳細はいつも通りコメントアウトで書いています。
今回は厳密解を求められる常微分方程式だったので、厳密解も同じグラフに描画するコードにしています。
import numpy as np
import matplotlib.pyplot as plt
#数値計算したい微分方程式の定義
def func(y):
return y
#オイラー法
def euler_method(f, y, h):
return y + f(y) * h
#常微分方程式を計算
def ode_solve(f, y0, T, N):
t = np.linspace(0, T, N)
h = t[1] - t[0]
y = [y0, ]
for i in range(N - 1):
y.append(euler_method(f, y[-1], h))
return t, np.array(y)
#初期条件などを設定して、上で定義した関数たちを呼び出す
y0 = 0.01
T = 10
N =101
t, y = ode_solve(func, y0, T, N)
y_exact = y0 * np.exp(t) #math.exp()だと配列ではないので計算不可
#以下は可視化のためにグラフを書くコード
fig, ax = plt.subplots()
ax.plot(t, y, color='crimson', label='Euler Method')
ax.plot(t, y_exact, color='cornflowerblue', label='Exact')
ax.set_xlabel('t')
ax.set_ylabel('y')
ax.legend()
ax.tick_params(direction='in', axis='x')
ax.tick_params(direction='in', axis='y')
#以下の範囲設定は適当です
ax.set_xlim(0, 5)
ax.set_ylim(0, 1)
plt.show()
fig.savefig("EulerMethod.jpg")実行結果
実行結果は下のグラフのようになると思います。オイラー法で求めたものは厳密解と比べると、誤差がことがわかります。これは、先ほど述べたようにオイラー法は $h^{2}$ 以下の項を無視していることが要因です。
最後に
今回はオイラー法を勉強してみました。実装自体は簡単でしたが、実用性には欠けそうですね。
次回以降は、ルンゲ-クッタ法なども勉強していこうかなと思っています。
今回はこの辺で終わりにします。
他の記事もぜひご覧ください！
コメント