はじめに

常微分方程式とは独立変数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の中の各種パラメータを変えると系の特性が見えてきます。