微分方程式(オイラー法とルンゲ・クッタ法)

大学・大学院生活

方程式?

一般的に方程式といえば以下のような式を思い浮かべます。

方程式は 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で実行しましたが、グラフの通り、オイラー法と比較しても、関数に近い値をとっていることがわかります。

以上

コメント

タイトルとURLをコピーしました