룽게-쿠타 방법

섬네일 미분방정식 사이드바 수치 해석에서 **룽게-쿠타 방법**(Runge-Kutta方法, Runge–Kutta method)은 적분 방정식 중 초기값 문제를 푸는 방법 중 하나이다.
# 역사
1900년 경 독일의 수학자 카를 룽게와 마르틴 빌헬름 쿠타가 개발하였다.
# 일반적인 4차 룽게-쿠타 방법
일반적으로 사용하는 룽게-쿠타 방법은 4차 룽게-쿠타 방법으로 보통 "RK4" 라고 쓴다. 별다른 수식어 없이 룽게-쿠타 방법을 쓴다고 말한다면 대체로 4차 룽게-쿠타 방법을 쓴다는 뜻이다. 다음과 같이 정의된 초기값 문제를 두자: : y' = f(t, y), \quad y(t_0) = y_0. y는 시간 t에 대한 미지의 함수이며, 우리가 근사하려는 것이다. y의 변화인 \dot{y}는 t와 y자신으로 이루어진 함수이고, 초기시간 t_0에 대응하는 y의 초기값은 y_0이며, 함수 f와 t_0, y_0의 값은 주어져 있다. 이제 h > 0 인 단계 크기로 다음을 정의한다. :\begin{align} y_{n+1} &= y_n + \tfrac{1}{6}h\left(k_1 + 2k_2 + 2k_3 + k_4 \right) \\\ t_{n+1} &= t_n + h \\\ \end{align} n = 0, 1, 2, 3, ... 에서, 다음을 사용한다. : \begin{align} k_1 &= f(t_n, y_n) \\\ k_2 &= f(t_n + \tfrac{1}{2}h, y_n + \tfrac{1}{2}h k_1) \\\ k_3 &= f(t_n + \tfrac{1}{2}h, y_n + \tfrac{1}{2}h k_2) \\\ k_4 &= f(t_n + h, y_n + h k_3) \end{align} 여기서 {y}_{n+1}이란 y(t_{n+1}) 의 RK4 근사 항이며, 다음 단계의 값({y}_{n+1})은 현재 단계의 값(y_n)에 네 구간의 크기 h의 증분치들의 가중평균을 더한 값이다. * k_1 은 구간의 시작부분의 기울기에 의한 증분값이며, y 를 사용한다.(오일러 방법) * k_2 는 구간의 중간의 기울기에 의한 증분값이며, y+\frac{h}{2}k_1 을 사용한다. * k_3 는 마찬가지로 구간의 중간의 기울기에 의한 증분값이나, y+\frac{h}{2}k_2 를 사용한다. * k_4 은 구간의 끝의 기울기에 의한 증분값이며, y+h k_3 를 사용한다. 네 증분을 평균을 낼 때, 중간 부분에서 더 많은 무게의 증분이 더해진다. 만약 함수 f 가 y 에 대하여 독립적이라면, 미분방정식은 단순한 적분이 되어 RK4는 심프슨 공식이 된다. 일반적인 미분방정식에서 RK4 방법과 다른 낮은 차원의 방법을 비교한 것이다. RK4 방법은 4차 해법이다. 이는 각 단계에서의 오류는 점근 표기법으로 O(h^5), 총 누적 오류는 O(h^4)라는 것을 의미한다.
# 양함수적(명시적) 룽게 - 쿠타 방법
양함수적(명시적, explicit) 룽게-쿠타 방법은 위에서 설명된 RK4 방법의 일반화이다. 이것은 다음과 같이 주어진다. :y_{n+1}=y_n+h\sum_{i=1}^s b_i k_i 이때 각 k_i 는 다음과 같다. :\begin{align} k_1&=f(t_n,y_n),\\\ k_2&=f(t_n+c_2h,y_n+h(a_{21}k_1)),\\\ k_3&=f(t_n+c_3h,y_n+h(a_{31}k_1+a_{32}k_2)),\\\ &\cdots\\\ k_s&=f(t_n+c_sh,y_n+h(a_{s1}k_1+a_{s2}k_2+\cdots+a_{s,s-1}k_{s-1})). \end{align} 특정한 방법을 지정하기 위해서는, 정수 *s* (각 단계들의 수) 가 주어져야 하고, 각각의 계수들 *aij* (1 ≤ *i* p 이다. s단계 explicit 룽게-쿠타 방법이 열린문제에서 p의 차수를 가지기 위해서는 최소한의 s가 필요하다. 몇 개의 알려진 값들은 다음과 같다. :\begin{array}{c|cc} p&1&2&3&4&5&6&7&8 \\\ \hline \min s&1&2&3&4&6&7&9&11\end{array}
## 예시
RK4는 이 방법에 맞아들어가며, 그 테이블은 다음과 같다. :\begin{array}{c|cc} 0 \\\ 1/2&1/2 \\\ 1/2&0&1/2\\\ 1&0&0&1\\\ \hline &1/6&1/3&1/3&1/6\end{array} 이 룽게-쿠타 방법에서의 약간의 변형을 준 방법은 쿠타에 의해서 연구되었고 3/8법칙 이라고 불린다. 이 방법의 가장 큰 장점은 다른 일반적인 방법들보다 오차계수가 작은 것이나, 이 방법은 매 단계마다 조금 더 많은 FLOP(부동소수점 연산)들이 필요하다. 그 Butcher 테이블은 다음과 같다. :\begin{array}{c|cc}0\\\1/3&1/3\\\2/3&-1/3&1\\\1&1&-1&1\\\ \hline &1/8&3/8&3/8&1/8\end{array} 그러나, 가장 간단한 룽게-쿠타 방법은 (전방)오일러 방법이다. 공식은 y_{n+1}=y_n+hf(t_n,y_n)이다. 이것은 오직 한 단계 뿐인 explicit 룽게-쿠타 방법이다. 대응되는 테이블은 다음과 같다. :\begin{array}{c|cc} 0 \\\ \hline &1\end{array}
## 두 단계의 이차 방법들
간단한 예로 두 단계의 2차 방법인 중간점 방법을 보면 다음과 같다. :y_{n+1}=y_n+hf\biggl(t_n+\frac12h,y_n+\frac12hf(t_n,y_n)\biggr). 대응되는 테이블은 다음과 같다. :\begin{array}{c|cc} 0\\\ 1/2&1/2\\\\\hline&0&1\\\ \end{array} 두 단계의 2차 룽게-쿠타 방법은 중간점 방법 밖에 없는 것은 아니다. 이런 종류의 방법은 α에 의해 다음과 같이 매개변수화되어 있다.y_{n+1}=y_n+h\Biggl(\biggl(1-\frac1{2\alpha}\biggr)f(t_n,y_n)+\frac1{2\alpha}f\bigl(t_n+\alpha h,y_n+\alpha hf(t_n,y_n)\bigr)\Biggr) 그 Butcher 테이블은 다음과 같다. :\begin{array}{c|cc}0\\\\\alpha&\alpha\\\ \hline&(1-\frac1{2\alpha})&\frac1{2\alpha}\\\ \end{array} 이 때 \alpha=\frac12일 때는 중간점 방법이 되며, \alpha=1일 때는 호인의 방법이 된다.
# 활용
예를 들어 α=2/3인 두 단계의 2차 룽게 쿠타 방법을 보자, 이것은 라슬톤 방법 이라고 알려져있다. 테이블은 다음과 같이 주어져 있다. :\begin{array}{c|cc}0 \\\ 2/3&2/3\\\\\hline&1/4&3/4\\\ \end{array} 대응되는 등식은 다음과 같다. :\begin{align}k_1&=f(t_n,y_n), \\\ k_2&=f(t_n+\frac{2}{3}h,y_n+\frac23hk_1),\\\y_{n+1}&=y_n+h\Bigl(\frac14k_1+\frac34k_2\Bigr). \end{align} 다음의 초기치 문제를 풀기 위하여 이 방법을 사용한다. :y^\prime=\tan(y)+1,\quad y_0=1,t\in[1,1.1] *h* = 0.025의 크기로 구간을 잡으면 이 방법은 네 단계가 필요하다. 이 방법의 결과는 다음과 같다: :\begin{array}{l} t_0=1:\\\&y;_0=1\\\ t_1=1.025:\\\&y;_0=1\qquad\qquad\qquad k_1=2.557407725\quad k_2=f(t_0+\frac23h,y_0+\frac23hk_1)=2.713898184\\\ &y;_1=y_0+h(\frac14k_1+\frac34k_2)=\underline{1.066869388}\\\ t_2=1.05:\\\&y;_1=1.066869388\quad\ k_1=2.813524695\quad k_2=f(t_1+\frac23h,y_1+\frac23hk_1)\\\ &y;_2=y_1+h(\frac14k_1+\frac34k_2)=\underline{1.141332181}\\\ t_3=1.075:\\\&y;_2=1.141332181\quad\ k_1=3.183536647\quad k_2=f(t_2+\frac23h,y_2+\frac23hk_1)\\\ &y;_3=y_2+h(\frac14k_1+\frac34k_2)=\underline{1.227417567}\\\ t_4=1.1:\\\&y;_3=1.227417567\quad\ k_1=3.796866512\quad k_2=f(t_3+\frac23h,y_3+\frac23hk_1)\\\ &y;_4=y_3+h(\frac14k_1+\frac34k_2)=\underline{1.335079087}.\\\ \end{array} 수치해석적 결과는 밑줄친 값들이다.
# 적응적(Adaptive) 룽게 - 쿠타 방법
Adaptive한 방법은 단일 룽게-쿠타 방법의 지역적인 절단 오차의 추정치를 찾기 위해 고안되었다. 이것은 테이블에서 차수가 p와 p-1인 두 방법들을 가져서 이루어졌다. 낮은 차수의 단계는 다음과 같이 주어진다. : y^*_{n+1}=y_n+h\sum_{i=1}^sb^*_ik_i, k_i가 높은 차원의 방법에서도 같다면 오차는 다음과 같다. : e_{n+1}=y_{n+1}-y^*_{n+1}=h\sum_{i=1}^s(b_1-b^*_i)k_i, 오차는 O(h^p)이다. 이런 종류의 Butcher 테이블은 b^8_i의 값을 얻기 위하여 확장 되었다. : \begin{array}{c|cc}0\\\ c_2&a;_{21}\\\ c_3&a;_{31}&a;_{32}\\\ \vdots&\vdots&&\ddots\\\ c_s&a;_{s1}&a;_{s2}&\cdots&a;_{s,s-1}\\\\\hline &b;_1&b;_2&\cdots&b;_{s-1}&b;_s\\\&b;^*_1&b;^*_2&\cdots&b;^*_{s-1}&b;^*_s\\\ \end{array} 룽게-쿠타-펠베르크 방법은 차수가 5와 4인 두 방법을 가지고 있다. 그 확장된 부처 태블로는 다음과 같다: : \begin{array}{l|ll}1\\\1/4&1/4\\\ 3/8&3/32&9/32\\\ 12/13&1932/2197&-7200/2197&7296/2197\\\ 1&439/219&-8&3680/513&-845/4104\\\1/2&-8/27&2&-3544/2565&1859/4104&-11/40\\\\\hline&16/135&0&6656/12825&28561/56430&-9/50&2/55\\\&25/216&0&1408/2565&2197/4104&-1/5&0\\\ \end{array} 하지만 가장 간단한 adaptive 룽게-쿠타 방법은 2차인 호인의 방법과 1차인 오일러 방법과 결합되어 있다. 그 확장된 Butcher 테이블은 다음과 같다: : \begin{array}{c|cc}0\\\1&1\\\\\hline&1/2&1/2\\\&1&0\\\ \end{array} 오차 추정은 단계 크기를 제어하기 위해 사용된다. 다른 adaptive 룽게-쿠타 방법은 Bogacki–Shampine 방법(3과 2의 차수), Cash–Karp 방법과 Dormand–Prince 방법(둘 다 5,4의 차수)이 있다.
# Implicit 룽게 - 쿠타 방법
위에서 언급된 모든 룽게-쿠타 방법들은 모두 explict 방법들이다. explict 룽게-쿠타 방법은 일반적으로 stiff equation에는 맞지 않는다. 왜냐하면 그들의 절대적인 안정성의 범위가 작기 때문이다; 특히, 그 식들은 한정되어 있다. 이 문제는 편미분방정식의 솔루션에서 특히 중요하다. explict 룽게-쿠타 방법의 불안정성은 음함수적(암시적, implicit)한 방법을 만들게 되었다. implicit 룽게-쿠타 방법은 다음의 형태를 가진다. : y_{n+1}=y_n+h\sum_{i=1}^sb_ik_i 이 때, : k_i=f\Biggl(t_n+c_ih,y_n+h\sum_{j=1}^sa_{ij}k_j\Biggr),\quad i=1,\dots,s. explict 방법간의 차이는 explict 방법은 시그마의 위에 *i* - 1이 들어간다는 것이다. 이것은 Butcher 테이블에서도 나타난다: explict 방법의 계수들의 행렬은 하 삼각 행렬이다. implicit 방법은 시그마 위에 *s*가 올라가고 계수들의 행렬은 다음과 같이 삼각행렬이 아니다. : \begin{array}{c|cc}c_1&a;_{11}&a;_{12}&\cdots&a;_{1s}\\\ c_2&a;_{21}&a;_{22}&\cdots&a;_{2s}\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\ c_s&a;_{s1}&a;_{s2}&\cdots&a;_{ss}\\\ \hline &b;_1&b;_2&\cdots&b;_s\\\ &b;^*_1&b;^*_2&\cdots&b;^*_s \\\\\end{array} =\begin{array}{c|c}c&A;\\\\\hline&b;^T\\\\\end{array} 이 차이는 매 단계마다, 대수적인 방정식의 시스템을 풀어야 한다는 점이다. 이것은 계산비용을 상당히 증가시킨다. *m*개의 성분이 있는 미분방정식을 *s*단계의 방법을 사용하면 대수적인 방정식은 *ms*개의 성분이 있다.이것은 implict한 선형 다단계 방법(ODE들의 다른 방법)과 대조된다: implicit한 *s*단계의 선형 다단계 방법은 *m*개의 성분을 가진 방정식을 풀기 위해서 단계수가 늘어나도 시스템의 크기가 늘어나지 않는다
## 예시
가장 간단한 implicit 룽게-쿠타 방법은 역 오일러 방법이다: : y_{n+1}=y_n+hf(t_n+h,y_{n+1}). Butcher 테이블은 다음과 같다. : \begin{array}{c|c}1&1\\\\\hline&1\\\\\end{array} 이 Butcher 테이블에 대응하는 식은 다음과 같다. : k_1=f(t_n+h,y_n+hk_1)\quad and\quad y_{n+1}=y_n+hk_1 이것은 위에서 나온 역 오일러 방법의 형태로 바꿀 수 있다. 다른 implicit 룽게-쿠타 방법의 예는 사다리꼴 법칙이다. 그 Butcher 테이블은 다음과 같다. : \begin{array}{c|cc}0&0&0\\\1&\frac12&\frac12\\\\\hline&\frac12&\frac12\\\&1&0\\\\\end{array} 사다리꼴 법칙은 배열 방법이다. 모든 배열 방법들은 모두 implicit 룽게-쿠타 방법이지만 모든 implicit룽게-쿠타 방법은 배열 방법이 아니다. Gauss-Legendre 방법은 가우스 구적법에 기반한 배열 방법이다. *s*단계의 Gauss-Legendre 방법은 차수가 2*s*이다(따라서 임의의 높은 차수의 방법들이 구성될 수 있다).두 단계의 방법(차수는 4이다)의 Buther 테이블은 다음과 같다: : \begin{array}{c|cc}\frac12-\frac16\sqrt3&\frac14&\frac14-\frac16\sqrt3\\\\\frac12+\frac16\sqrt3&\frac14+\frac16\sqrt3&\frac14\\\\\hline&\frac12&\frac12\\\&\frac12+\frac12\sqrt3&\frac12-\frac12\sqrt3\\\\\end{array}
## 안정도
explicit한 룽게-쿠타 방법에 비교하여 implicit 방법의 장점은 안정도가 특히 딱딱한 방정식에서 높다는 것이다. 다음의 선형 방정식 *y*'=λ*y*를 고려해 보자. 룽게-쿠타 방법을 가용하면 이 식은 다음y_{n+1}=r(h\lambda)y_n의 반복으로 줄며, *r*은 다음과 같다. :r(z)=1+zb^T(I-zA)^{-1}e=\frac{\det(I-zA+zeb^T)}{\det(I-zA)}, 여기서 *e*는 1들의 벡터를 나타낸다. 함수 *r*은 *안정도 함수*라고 불린다. 만약 이 방법이 *s*단계라면 *r*은 차수 *s*의 두 다항식의 몫의 형태를 띈다. explicit 방법은 순상삼각행렬인 행렬 *A*를 가진다. 이것은 det(*I*-*zA*) = 1이라는 것을 의미하며 그 안정도 함수는 다항함수이다. 위의 선형방정식은 |r(z)|<1이고 z=h\lambda일 때, 수치적 해법은 0이 된다. 이런*z*들의 집합은 *절대 안정성의 영역*이라고 부른다. 특히 이 방법은 모든 실수부가 음인 z가 절대 안정성의 영역 안에 있을 때 A-안정하다고 불린다. explicit 룽게-쿠타 방법의 안정도 함수는 다항함수이기 때문에 explicit 룽게 쿠타 방법은 A-안정할 수 없다.
# 알고리즘

# 코드
* C/C++ 코드 //C/C++ code of 4th order Runge-Kutta method # include # include # include float f(float t,float y) { return(t*y+1); } void main() { float t0,y0,t1,y1,h,i,at,k1,k2,k3,k4; clrscr(); printf("\n Enter Value Of t0 and y0"); scanf("%f %f",&t0;,&y0;); printf("\n Enter Value of h"); scanf("%f",&h;); printf("\n Enter Final Value of t"); scanf("%f",&at;); printf("\nt\ty"); printf("\n_______________________________________\n"); do { k1=h*f(t0,y0); k2=h*f(t0+h/2,y0+k1/2); k3=h*f(t0+h/2,y0+k2/2); k4=h*f(t0+h,y0+k3); y1=y0+h/6*(k1+(2*k2)+(2*k3)+k4); t1=t0+h; printf("\n%.4f\t%.4f",t1,y1); t0=t1; y0=y1; }while(t0 * Python # Program Of Runge-Kutta 4th Order from math import * def f(t, y): return sin(y**t) def runge_kutta(): t0, y0 = float(input("initial time:")), float(input("initial value:")) h = float(input("interval size:")) t = float(input("terminal time:")) print("t y") print("______________") while t0 < t: k1 = f(t0, y0) k2 = f(t0 + h / 2, y0 + (k1 * h) / 2) k3 = f(t0 + h / 2, y0 + (k2 * h) / 2) k4 = f(t0 + h, y0 + (k3 * h)) y1, t1 = y0 + h / 6 * (k1 + (2 * k2) + (2 * k3) + k4), t0 + h print(f"{t1:5.4f}, {y1:5.4f}") t0, y0 = t1, y1 if __name__ == "__main__": runge_kutta()