Euler method로 미분방정식 풀어보기

미분방정식을 컴퓨터를 이용해 수치적으로 계산해야 할 때가 있다. 이 때 미분방정식을 푸는 방법이 여러개 있는데 가장 쉬운 방법은 Euler method이다. 하지만 에러가 꽤 커서 실제 연구할 때에는 사용하지 않는 방법이기도 하다. 미분의 정의로부터 euler method를 간단히 유도할 수 있다. 미분의 정의에서 극한을 빼면, 

처럼 쓸 수 있다. 물리적인 상황을 떠올리기 위해 위치를 s(t), 속도를 v(t)로 나타내었다. 시간에 대한 위치의 변화량이 속도가 된다. 이 식을 다시 고쳐보면 다음과 같다.

다음 시간(t+dt)에서의 위치는 지금 시간의 위치에서 지금 시간의 속도만큼의 변화를 주어 예측하는 방법이 된다. 이게 바로 Euler method이다. 사실 dt에 대해 테일러 전개를 해보면, 

에서 최대 (dt)의 제곱항만큼의 오차가 생기게 된다. 

이제 Euler method를 프로그램으로 만들어보자. 그 전에 생각해볼 문제가 있다. 바로 euler method를 어떻게 반복시킬 것인가 하는 문제이다. 지금의 위치와 속도를 통해 다음 시간 스텝에서의 위치를 예측하는 것이 이 method의 목적임을 상기하자. 그럼 다음과 같이 순서를 짜야 한다.

1. 우선 현재 시간 t 에서의 위치 x를 보자. 

2. 그 뒤에 현재 시간 t에서의 속도 v(t)를 계산해야 한다. 

3. 그럼 현재 속도 v(t)와 x를  Euler method에 넣어 t+dt 때의 위치를 어림계산할 수 있다. (물론 오차는....ㄷㄷ) 

4. 이후에 이 때의 위치 x+dx 가 t+dt 때의 시간이라는걸 컴퓨터에게 알려준다.

이걸 반복하면 된다. 검은 화살표는 약간의 연산을 해야 하는 과정이다. (속도구하기, Euler method)

녹색은 그냥 컴퓨터에서 대입을 해주는 부분이다.

프로그래밍을 편하게 하기 위해서는 미리 자신이 원하는 dt(시간 간격)와 최종 시간을 정한 뒤 배열의 크기를 잡아준다. 배열의 크기는 (전체 시간 / 시간 간격) 으로 정할 수 있다. 예를들어 2초 동안을 0.1초로 나누어 보고 싶다면 2/(0.1) = 20개의 조각으로 시간이 나뉘게 된다.

이후 결과를 츨력할 때 편하게 하기 위해 초기값을 배열의 맨 앞에 넣어준다. 그럼 배열의 크기 역시 1개가 더 필요한데, 이를 위해 처음부터 배열의 크기를 (전체시간 / 시간간격) + 1 로 잡아주자. 위의 예시에서 배열의 크기는 21이 된다.

아래는 python으로 Euler method를 이용해 지구에서 자유낙하 운동을 만든 예제이다. 질량은 1이고 위로 30m/s 만큼 던졌을 때의 시간에 대한 모습을 보여준다. (NumPy, matplotlib 설치 필요)

import pylab as pl

g = 9.8

v0 = 30

def velocity(t):

return v0 - g * t

timeEnd = 8 ## second

dt = 0.01

g = 9.8

Nt = (int)(timeEnd/dt)+1

yList = pl.empty(Nt, float)

tList = pl.empty(Nt, float)

yList[0] = 60

tList[0] = 0

for i in range(Nt-1):

yList[i+1] = yList[i] + dt * velocity(tList[i])

tList[i+1] = tList[i] + dt

pl.plot(tList, yList)

pl.show()

계산을 해보면 7.71 ~ 7.72초 사이일 때 땅에 떨어진다고 하는데 등가속도 운동을 계산해보면 약 7.7105 초에 정확히 지면에 떨어진다. 

다음으로 Euler method의 오차를 생각해볼 수 있는데, 한 step에서 계산 결과는 (dt)^2 만큼 쌓이게 된다. 하지만 우리가 최종결과를 원래 이론값과 비교할 때의 오차는 조금 다르다. step마다 계산할 때 생기는 오차를 local error라고 부르며, 우리가 원하는 시간 T에서의 계산값과 이론값의 차이는 global error라고 부른다. 그래서 global error는 local error보다 더 크기 마련인데, 정확히는 (dt)^1 만큼이다. 왜일까? 한 step에서 (dt)^2 만큼 오차가 발생하는데 이 오차가 T/dt = N번 누적된다. 여기서 N은 결국 1/dt에 비례하고, global error는 local error가 N번 누적된 것이므로, 

만큼의 오차가 생긴다. 이를 확인해보자. dt를 바꾸어가면서 어떤 미분방정식을 풀어본 뒤에 log-scale로 global error 와 dt에 대한 그래프를 그려보면 된다. 분명 기울기가 1이 나온다. 위의 프로그램을 조금 고쳐서 확인해보자. T=8일 때의 위치는 정확히 -13.6m 이다.

import pylab as pl

from math import pow, fabs

g = 9.8

v0 = 30

def velocity(t):

return v0 - g * t

timeEnd = 8 ## second

GlobalError = []

dtList = []

for k in range(4):

dt = pow(0.1, k)

dtList.append(dt)

Nt = (int)(timeEnd/dt)+1

yList = pl.empty(Nt, float)

tList = pl.empty(Nt, float)

yList[0] = 60

tList[0] = 0

for i in range(Nt-1):

yList[i+1] = yList[i] + dt * velocity(tList[i])

tList[i+1] = tList[i] + dt

GlobalError.append( fabs(yList[Nt-1] + 13.6) )

pl.plot(dtList, GlobalError)

pl.yscale('log')

pl.xscale('log')

pl.show()

결과는 다음과 같다. 잘 보면 기울기가 log-log에서 -1이다. 

Global Error를 찍을 때 가끔 결과가 엉뚱하게 나오는 때가 있다.(기울기가 아무리해도 -1인 경우...심지어는 Runge-kutta 4th를 써도!!) 이번에 알게 되었는데 Nt의 크기에서 +1을 하지 않았을 경우 정확히 T=8 일 때가 아닌 T = 8-dt 때의 결과로 에러를 계산하기 때문이다.

이렇게 Euler method를 살펴보았다. 이해하기도 쉽고 만들기도 쉽지만 에러는 우리가 준 dt만큼밖에 줄어들지 않기 때문에 긴 시간뒤의 결과를 실제와 다르게 내주게 된다. 그렇다고 dt를 아주 작게 쪼갤 경우 loop를 많이 돌게 되어 컴퓨터 파워를 많이 써야 하는 부담이 있다. 이런 이유 때문에 Euler method를 상당히 초기조건에 민감한 비선형 미분방정식을 수치적으로 푸는데 적합하지 않다. 이를 극복하기 위해 Improved Euler, Leap-frog, Runge-Kutta 4th order 등의 여러 알고리듬이 있다.