소개

상미분 방정식은 독립 변수 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 방법을 파이썬으로 구현하면 다음과 같습니다.

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안의 각종 파라미터를 바꾸면 계의 특성이 보입니다.