*Runge-Kutta method (설명, 예제)

Numerical Analysis (수치 해석학): 어떤 함수나 방정식의 해를 컴퓨터를 이용해 수치적으로 근사해서 근삿값을 구하는 알고리즘에 대한 연구를 하는 학문

수치해석을 사용하는 이유는 사람의 손으로는 풀 수가 없는 문제를 풀기 위해서 사용하며, 특히 미분 방정식이나 다양한 문제를 컴퓨터를 이용하여 풀 수 있다는 장점이 존재합니다. 

 

Runge-Kutta method 

독일의 수학자 카를 `다비트 톨메 룽게`와 `마르틴 빌헬름 쿠타`가 개발한 수치 해석학의 알고리즘 중 하나입니다. 상미분 방정식의 해를 구하기 위해 고차원의 미분 식이 필요 없이 정확하게 근사할 수 있는 방법 중 하나입니다. 

장점: Tayler series와 같이 고차원의 미분이 필요 없으며, 상당수의 temporal discretization 상미분 방정식을 풀 수 있음. 

단점: stiff equation에 대해서는 보통 효과적이지 못하며, stable 하지 못한 결과를 가져옴. 계산 과정이 다른 과정보다 많이 시간 소모가 필요함. 

일반적으로 `룽게-쿠타 방법` 또는 보통 4차 식을 활용하기 때문에 `RK4`라고 간단하게 쓰기도 합니다. RK1 이나 RK2와 같이 차수를 적게 하여서 사용이 가능하지만 정밀도 (근사)는 차수에 따라서 더 증가하기 때문에 RK4를 제일 많이 사용합니다.  

*Euler method = 1st Runge Kutta,    Improved Euler method = 2nd Runge Kutta. 

RK는 오일러 방법 (Euler method)보다 더 정확하게 근사를 할 수 있으며 쉽게 implementation이 가능하다는 특징이 있습니다.

오일러 방법의 전개 중에서 아래와 같이 dx^2에 따라서 비례하여 오차가 누적되게 됩니다. 

이를 개선시키기 위해서 고차의 항까지 고려를 하여서 식을 사용하게 되면, 즉 dx^n 에서 n이 커질수록 오차가 작아지게 되기에 해당 아이디어에서 시작되었습니다 (여기서 dx는 보통 작다고 가정을 하고 시작됩니다).

RK는 아래의 그림과 같이 4개 지점 (dt = 0, dt = h/2, dt = h) 의 기울기 평균을 통해서 다음 step의 지점을 찾아내는 계산 방법입니다. 4개의 slope를 토대로 평균을 내어서 예를 들어서 속도를 고려한 위치 변화를 계산할 때 아래의 오른쪽과 같이 slope를 근사 시킬 수가 있습니다.

[1] RK concept	[2] RK Overview

우리가 풀고자 하는 문제가 아래와 같은 초기값 문제일때, y는 시간 t에 대한 미지의 함수이며 우리가 근사를 하려는 식입니다. y의 변화인 y dot은 t와 y 자신으로 이뤄졌으며 대응하는 y의 초기값은 y_0입니다. 

 이때 해당 식에 근사하기 위해서 RK는 아래와 같이 사용됩니다. 

h는 step size이며 Sampling time과 비슷한 의미입니다. step size가 작아질수록 정밀한 근사가 가능해집니다. 아래의 표는 다른 수치해석 방법과 비교해서 step size에 따른 global error 결과를 표시한 표입니다. 결과와 같이 모든 방법이 step size가 작아질수록 error값은 작아지며,  RK4의 결과가 제일 낮은 에러 값을 가집니다. 

The maximum global error committed by Euler, Heun, and RK4 methods [3]

각각의 k 값들은 증분값 (Increment value)을 의미하며, k_1은 구간의 시작 부분의 기울기에 의한 증분 값을 의미하기에 Euler method와 동일합니다. 

k_2는 구간의 중간 기울기에 의한 증분 값이며, y+h/2 * k_1를 사용합니다. 

k_3는 구간의 중간 기울기에 의한 증분값이며 k_2와 동일한 의미이지만 y+h/2 * k_2를 사용합니다. 

k_4는 구간의 끝의 기울기에 대한 증분값이며 y+hk_3를 사용합니다. 

즉 네 증분 (k_1,2,3,4)의 평균을 통해서 아래와 같은 식으로 표현이 되어서 최종식에 더해지는 형식으로 진행됩니다. 

 

 

RK Example 

MATLAB을 이용하여서 다음과 같은 문제를 풀어봅시다. 

본 함수는 y(x)이며, 초기값은 x=0에서 y값은 -1을 가지는 방정식입니다. 해결을 위해서 먼저 미분 방정식을 저는 따로 MATLAB의 function 파일을 하나 만들어주었습니다. 

이름은 target_func이며 dY를 계산하여 넘겨주는 함수입니다. 

function dY = target_func(t, Y) 
dY = sin(t);
end

완료가 되었다면 아래와 같이 기본적인 변수들 정보와 RK 식을 그대로 for문에 implementation 해주면 됩니다. 

clc;
clear all;
close all force;

Y0 = -1;

t_start = 0;
t_end = 6.24;

dt = 0.01;
n = t_end/dt + 1; 
Y = Y0;

Y_data = zeros(1,n); 
Y_data(1,1) = Y;
t_data = zeros(1,n); 

% exact solution
fn_y = @(t) -cos(t);
exact_value = zeros(1,n); 
exact_value(1,1) = Y;

% iteration start
for i = 2:n 
   t = t_start + (i-1) * dt; 
   
   K1 = target_func(t, Y); 
   K2 = target_func(t + dt/2, Y + K1 * dt/2); 
   K3 = target_func(t + dt/2, Y + K2 * dt/2); 
   K4 = target_func(t + dt, Y + K3*dt); 
   
   Y = Y + dt*(K1 + 2*K2 + 2*K3 + K4)/6; 
   
   Y_data(:,i) = Y; 
   t_data(:,i) = t;
   
   exact_value(:,i) = fn_y(t); 
    
end

figure; hold on; grid on;
plot(t_data, exact_value, 'black-','LineWidth',2)
plot(t_data, Y_data, 'r--','LineWidth',2)
legend('Exact solution','RK 4th order')
xlabel('time (s)'); ylabel('y(t)');

결과는 아래와 같습니다. 

 

Reference

1. http://www.physics.drexel.edu/~steve/Courses/Comp_Phys/Integrators/rk4.html

2. https://www.haroldserrano.com/blog/visualizing-the-runge-kutta-method

3. Runge-Kutta Method - Handbook of Dynamical Systems, 2002