미분방정식에 대한 수치해석학적 해(Runge-Kutta) - 원리

   이 방법은 원리는 조금 어려운데 알고리즘의 구현은 어렵지 않다. 따라서 원리보다는 알고리즘이 문제가 된다면 우선 알고리즘만 보면 된다. 또한, C++의 가상함수를 사용하여 조금 일반화시켰다. 만약 수천수만번 미방을 풀어야 한다면 이렇게 일반화된 방법은 좀 많이 느리므로^[각주:1] 직접 코드를 구현해 사용할 것을 추천한다.

글이 너무 길어서 두 부분으로 나눈다. 구현 부분은 따로 뺀다.

   풀어야 할 문제는 주어진 함수의 미분된 식을 알고 있을 때, 미분되지 않은 식을 찾아 내는 것이다. 즉,

dy/dt = f(t, y(t))

로 y의 t 에 대한 미분식이 주어져 있을 때 y(t) 를 찾아 내는 것이 목표이다. 예를 들면 이렇다.

y' = t - 1

일 때, y는 무엇일까? 이것은 매우 간단하다. 양변을 t 로 적분하고, 적분상수는 초기 조건, 즉, 흔히 y(0) = k 로 주어지는 것으로부터 알 수 있다. 그러나 만약

y' = exp(-t^2) + y*ln(y)

로 되어 있다면 어떻게 할까? 우리는 몇 가지 경우에 대하여 알려진 몇 가지 해결책을 미분방정식 수업 시간에 배운다. 그러나 실제 상황에서 만나게 되는 문제들은 그러한 analytic solution 이 존재하는 경우는 드물다. 그렇게 딱 떨어지는 답을 구할 수 없는 경우, 수치해석학적으로 다음과 같이 구할 수 있다.

유도 과정에 나오는 식 중 최종식은 나도 직접 계산해 보지 않았다. 따라서 이 글은 기본 원리 정도를 설명하는 것을 목표로 한다.

문제 dy/dt = f(t, y), a<= t <= b, y(a) = k 가 주어져 있을 때, a <= t <= b 에서의 y(t) 를 구한다.

즉, t 에 대한 함수 y 에 대하여, y 를 t 로 미분한 식이 구간 [a,b] 에 대하여 주어져 있을 경우, 그 구간 내에서 y(t)를 찾아 내는 것.

대부분의 알고리즘이 테일러 전개에서 출발을 한다. 원리는 다음과 같다.

위에서 굵고 검은 선을 생각해 보자. 아직 우리는 이 함수 y 를 모른다고 가정하자. t_i 에서의 값 y(t_i)를 알고 있다고 가정하자. 또한, 우리에게는 y'는 주어져 있기 때문에 t_i 에서이 미분값 즉 dy/dt|t=t_i 값을 알고 있다. 그렇다면 t를 t_i 에서 살짝 옆으로 간 t_i+1 에서의 y 값인 y(t_i+1) 은 무엇일까? 즉, 다음과 같다.

는 주어진 조건이고(즉, 우리는 함수의 미분식을 알고 있으므로 각 시각에서의 미분값이 주어진 것이다),

는 미분(기울기)의 정으로부터 유도되는 식이다. 따라서 위의 두 식을 합하면

이고, 우리가 알고 싶은 것은 y(t_i+1) 이므로 위 식을 다음과 같이 정리한다.

그리하여, 우리는 한 점 t_i 에서의 y(t_i) 를 알면 다음 점 t_i+1 에서의 y(t_i+1) 값을 위와 같이 찾아낼 수 있는 것이다. 바로 그것 때문에 조건에 a <= t <= b 라는 t 에 대한 구간과, y(a) = k 라는 초기 조건이 필요한 것이다. 따라서 시작 구간에서의 y 값을 알고 있으므로, 이 값을 이용하여 위의 방법을 이용해서 시작 구간 a 의 다음 점 a + h 에서의 y의 값 y(a+h) 를 구할 수 있는 것이다. h 가 작으면 작을수록 오차가 줄어들 수 있겠지만^[각주:2] 계산 시간이 오래 걸릴 것이다. 위 식에서 h 는 t_i+1 - t_i 로, 즉 step size 인 것이다.

위의 방법을 아주 약간 formal 하게 표현하면 다음과 같다.

주어진 조건

에 대하여 테일러 전개를 하여 다음과 같이 표현할 수 있고,

이 식에서 2차 항을 없애면

위와 같은 식을 얻게 된다. 바로 위 식이 Euler 방법이다. 역시 h 는 step size가 된다. 이 방법은 간단한만큼 오차가 좀 크다. 특히 y 값이 갑자기 크게 변하는 곳에서 그 오차가 심하게 크게 나온다. 이 오차를 줄이기 위해서는 테일러 전개시 2차 이상까지 전개한 이후 그것들까지 쓰면 된다. (왜 이 방법은 많이 안 쓰는 것일까?) 하지만 다음과 같이 2변수에 대한 테일러 전개를 이용하면 더 쉽게 해결할 수 있게 된다. 이렇게 다변수 테일러 전개를 이용하는 방법이 Runge-Kutta 방법인 것이다.

2변수에 대한 테일러 전개는 다음과 같다. 함수 f(t,y) 를 t 와 y 로 테일러 전개를 해서 다음과 같은 형태로 만들 것이다.

위의 식은 바로

을 간략히 표현한 것이다. 위가 바로 함수 f(t,y) 를 t, y 두 변수에 대하여 n 차까지 테일러 전개한 식이 되는 것이다. 이제 우리의 목표는

위처럼 2차 항까지 전개한 식에 대하여, 이 값과 매우 유사한 값을 만들어 주는 a_i, α₁,β₁ 을 찾는 것이다. 우선 위 식의 첫 번째 줄을 좀 더 풀어 보면 다음과 같다.

첫 번째 줄은 테일러 전개의 첫 번째 항까지 표현한 것이고,

두 번째 줄은 f 의 미분을 다시 표현해 본 것이고,

세 번째 줄은 f 의 미분을 chain rule 을 써서 풀어 쓴 것이고(현재 f 는 t 와 y 에 대한 함수이니까),

네 번째 줄은 dy/dt = f(t,y) 로 주어진 조건을 이용해서 세 번째 줄을 다시 쓴 것이다.

이제

이 식에서 두 번째 줄에 있는 a₁ f(t + α₁, y + β₁) 을 좀 더 풀어 써 보자.

이다. 이것은 편미분의 정의로부터 나오는 것이다. 이제 이 두 식을 합하면,

이다. 이것을 다시 정리해서, 위 식 중 첫 번째 줄과 마지막 줄만 남기면,

이다. 이 식에서 f(t,y)의 계수와 f 를 t 와 y 각각으로 편미분한 항의 계수를 맞추면,

가 나온다. 약간 정신이 없는데, 왜 이런 일을 했는지 잠깐 생각해 보자.

(음... tex 에 한글을 쓸 수 있도록 설정하지 않아서...) 오일러 방법과 Runge-Kutta 방법을 비교해 보자. 오일러 방법에서 y' 를 테일러 전개한 식으로부터 y(t_i+1) 을 찾아 내었다. 마찬가지로 Runge-Kutta 에서도 y' 를 찾아 내서 그 식으로부터 y(t_i+1) 을 찾아낼 것이다. y' 를 다변수로 2차 항까지 테일러 전개한 항과 나머지 항으로 나누어서 표현할 수 있다는 것이 제일 마지막 식의 아래서 세 번째 줄(P_n(t,y) + R_n(t,y)) 과 네 번째 줄이고, 우리는 이것을 제일 아래 식으로 approximation 시킨 후 a₁, α₁, β₁ 을 찾고자 하였다. 그 값들을 바로 위에서 찾았었다.  따라서 y'를 approximation 하기 위한 a₁ f(t+α₁, y+β₁) 을 이용해서, 오일러 때와 같은 식으로 y(t_i+1) 을 표현하면,

이 된다.

지금까지는 2변수 테일러 전개항의 2차항까지만 이용한 것을 살펴 보았다. 똑같은 방법으로 4차 항까지 전개하여 계수를 맞춰서 변수를 찾아내어 사용하는 것이 일반적으로 많이 사용하는 Runge-Kutta of order four 방법이다. 즉,

이 되는 것이다. 이 때, w 는 approximation 된 y 값이고 alpha는 시작 구간에서의 y 값이다.

1. C++의 virtual 함수는 생각보다 많이 느리더군. [본문으로]
2. h를 줄인다고 오차가 항상 주는 것이 아니다. 오차 항, 즉 테일러 전개 시 뒤쪽에 날려 버리는 항이 h 에 대하여 선형이 아니라면 h를 반으로 줄인다고 오차가 반이 되지는 않는다 [본문으로]