はじめに
常微分方程式とは独立変数t,とtに従属する関数y(t), 及び、そのn階導関数(n=0,1,2,...,N)を含む方程式のことです。
つまり、
g(t,y,y',y'',y^{(3)},...,y^{(N)}) = 0 \hspace{50pt}(1)
の形で記述可能な方程式のことです。
このような方程式はEular法(オイラー法)という数値解法により、解を近似できます。
実際の例を用いて説明していきます。
Eular法を使うための準備
上の微分方程式$(1)$は、実はオイラー法を使うには適していない形になっています。
これをどうにかして下の$(2)$のような形に変形しなければなりません。
\frac{dx(t)}{dt} = f(t,x(t)) \hspace{50pt}(2)
\\ただし, x(t)=\begin{pmatrix}y\\y'\\y''\\...\\y^{(N-1)}\end{pmatrix}
実際の例
ばね定数$k$のばねの先端に質量$m$の質点のついた系を考えます。系に振動的外力がかかり、また、空気抵抗を考慮する時の微分方程式は以下のようになります。
※解析的に解くことが可能ですが、数値計算法で解くことを考えます。
m\frac{d^2y(t)}{dt^2}+2\gamma \frac{dy(t)}{dt}+ky(t)=a\sin(\omega t)\hspace{50pt}(3)
$(3)$を変形すると、
\frac{d^2y(t)}{dt^2}=-\frac{2\gamma}{m} \frac{dy(t)}{dt}-\frac{k}{m}y(t)+\frac{a}{m} \sin(\omega t) \hspace{50pt}(3-1)
また、
\frac{dy(t)}{dt}=\frac{dy(t)}{dt}\hspace{50pt}(3-2)
ここで、ベクトル$x(t)$を
x(t)=\begin{pmatrix} x_1(t) \\ x_2(t) \end{pmatrix} = \begin{pmatrix} y(t) \\ \frac{dy(t)}{dt} \end{pmatrix}
と置くと、$(3-1)$,$(3-2)$から、
\frac{dx(t)}{dt}
=\frac{d}{dt} \begin{pmatrix} y(t) \\ \frac{dy(t)}{dt} \end{pmatrix}
=\begin{pmatrix} \frac{dy(t)}{dt} \\ \frac{d^2y(t)}{dt^2} \end{pmatrix}
=\begin{pmatrix} \frac{dy(t)}{dt} \\ -\frac{2\gamma}{m} \frac{dy(t)}{dt}-\frac{k}{m}y(t)+\frac{a}{m} \sin(\omega t) \end{pmatrix}
=\begin{pmatrix} x_2(t) \\ -\frac{2\gamma}{m} x_2(t)-\frac{k}{m}x_1(t)+\frac{a}{m} \sin(\omega t) \end{pmatrix}
よって、以下のように表せる。
\frac{dx(t)}{dt} = f(t,x(t))\hspace{50pt}(4)
\\where : f(t,x)=\begin{pmatrix} x_2 \\ -\frac{2\gamma}{m} x_2-\frac{k}{m}x_1+\frac{a}{m} \sin(\omega t) \end{pmatrix}
Eular法
ここまでの作業で、見栄えのいい形の微分方程式を作ることが出来ました。
では、ここからはEular法そのものについて説明します。
まず、微分の定義から
\frac{dx(t)}{dt} = f(t,x(t))
⇒\lim_{\Delta t \to 0} \frac{x(t+\Delta t)-x(t)}{\Delta t} = f(t,x(t))
ここで、十分小さい$\Delta t$を用いれば、
\frac{x(t+\Delta t)-x(t)}{\Delta t} \fallingdotseq f(t,x(t))
⇒x(t+\Delta t) \fallingdotseq x(t)+f(t,x(t))\Delta t\hspace{50pt}(5)
この式を使えば、最初の$x(t)$が分かれば$\Delta t$秒後の$x(t+\Delta t)$が計算できます。
プログラム(python)
Eular法をpythonで実装すると以下のようになります。
f.py
import numpy as np
a = 10 #a(外力の振幅)
(m,gm2,k) = (30,5,1)#m,2Γ,k
omg = np.pi #ω
def f(t,x):
f0 = x[1]
f1 = -(gm2/m)*x[1] -(k/m)*x[0] + (a/m)*np.sin(omg*t)
return np.array([f0,f1])
f_D.py
import numpy as np
from f import f as f
dt = 0.01
def f_D_eular(t,x):
return (t + dt, x + f(t,x)*dt)
main.py
import numpy as np
from f_D import f_D_eular as f_D
def main():
(t0,x0)=(0,np.array([1,0]))
t1 = 100
(t,x)=(t0,x0)
while(t < t1):
print(t,x[0])
(t,x)=f_D(t,x)
main()
という感じになります。
ちなみに、このプログラムの実行結果をグラフにすると、以下のようになりました。
初期条件である$t0,x0$やf.pyの中の各種パラメータを変えると系の特性が見えてきます。