룽게-쿠타 방법 문서 원본 보기

{{위키데이터 속성 추적}}
[[파일:Runge-kutta.svg|섬네일]]
{{미분방정식 사이드바}}
[[수치 해석]]에서 '''룽게-쿠타 방법'''(Runge-Kutta方法, {{llang|en|Runge–Kutta method}})은 [[적분 방정식]] 중 [[초기값 문제]]를 푸는 방법 중 하나이다.

== 역사 ==
1900년 경 독일의 수학자 [[카를 룽게]]와 [[마르틴 빌헬름 쿠타]]가 개발하였다.

== 일반적인 4차 룽게-쿠타 방법 ==
일반적으로 사용하는 룽게-쿠타 방법은 4차 룽게-쿠타 방법으로 보통 "RK4" 라고 쓴다. 별다른 수식어 없이 룽게-쿠타 방법을 쓴다고 말한다면 대체로 4차 룽게-쿠타 방법을 쓴다는 뜻이다.

다음과 같이 정의된 [[초기값 문제]]를 두자:
:<math> y' = f(t, y), \quad y(t_0) = y_0. </math>

y는 시간 t에 대한 미지의 함수이며, 우리가 근사하려는 것이다. y의 변화인 <math>\dot{y}</math>는 t와 y자신으로 이루어진 함수이고, 초기시간 <math>t_0</math>에 대응하는 y의 초기값은 <math>y_0</math>이며, 함수 f와 <math>t_0</math>, <math>y_0</math>의 값은 주어져 있다.

이제 h > 0 인 단계 크기로 다음을 정의한다.

:<math>\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}</math>

n = 0, 1, 2, 3, ... 에서, 다음을 사용한다.

:<math>
\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}
</math>

여기서 <math>{y}_{n+1}</math>이란 <math>y(t_{n+1})</math> 의 RK4 근사 항이며, 다음 단계의 값(<math>{y}_{n+1}</math>)은 현재 단계의 값(<math>y_n</math>)에 네 구간의 크기 h의 증분치들의 가중평균을 더한 값이다.
* <math>k_1
</math>은 구간의 시작부분의 기울기에 의한 증분값이며, <math>y
</math>를 사용한다.([[오일러 방법]])
* <math>k_2

</math>는 구간의 중간의 기울기에 의한 증분값이며, <math>y+\frac{h}{2}k_1
</math>을 사용한다.
* <math>k_3
</math>는 마찬가지로 구간의 중간의 기울기에 의한 증분값이나, <math>y+\frac{h}{2}k_2
</math>를 사용한다.
* <math>k_4
</math>은 구간의 끝의 기울기에 의한 증분값이며, <math>y+h k_3
</math>를 사용한다.
네 증분을 평균을 낼 때, 중간 부분에서 더 많은 무게의 증분이 더해진다. 만약 함수 <math>f
</math>가 <math>y
</math>에 대하여 독립적이라면, 미분방정식은 단순한 적분이 되어  RK4는 [[심프슨 공식]]이 된다.
[[파일:Runge Kutta method.png|대체글=File:Runge Kutta method.png|섬네일|일반적인 미분방정식에서 RK4 방법과 다른 낮은 차원의 방법을 비교한 것이다.]]
RK4 방법은 4차 해법이다. 이는 각 단계에서의 [[절단 오차|오류]]는 [[점근 표기법]]으로 <math>O(h^5)</math>, 총 누적 오류는 <math>O(h^4)</math>라는 것을 의미한다.

== 양함수적(명시적) 룽게 - 쿠타 방법 ==
양함수적(명시적, [[:en:Explicit_and_implicit_methods|explicit]]) 룽게-쿠타 방법은 위에서 설명된 RK4 방법의 일반화이다. 이것은 다음과 같이 주어진다.

:<math>y_{n+1}=y_n+h\sum_{i=1}^s b_i k_i</math>

이때 각 <math>k_i
</math>는 다음과 같다.

:<math>\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}
</math>

특정한 방법을 지정하기 위해서는, 정수 ''s'' (각 단계들의 수) 가 주어져야 하고, 각각의 계수들 ''a<sub>ij</sub>'' (1 ≤ ''i'' < ''j'' ≤ ''s''), ''b<sub>i</sub>'' (''i'' = 1, 2, ..., ''s''), ''c<sub>i</sub>'' (''i''=2, 3, ..., ''s'') 가 필요하다. 행렬 [''a<sub>ij</sub>'']는 ''룽게-쿠타 행렬''이며, ''b<sub>i</sub>'' 와''c<sub>i</sub>''  는 알려져 있듯이 ''무게'' 와 ''마디''이다.

이 값들은 ''[[부처 태블로]]''(Butcher tableau)이라는 기억장치 내에서 배열되어있다. [[존 찰스 부처]]의 이름을 딴 것이다.

:<math>\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 \\ \end{array}
</math>

룽게-쿠타 방법은 다음을 만족할 때 일관성이 있다.

:<math>\sum_{j=1}^{i-1}a_{ij}=c_i\ for\ i=2,\dots,s.
</math>

만약 하나가 특정 차수 ''p'',즉 지역 절단 오차가 <math>O(h^{p+1})
</math>이 되도록 이 방법이 요구되면 거기에는 또한 부수적인 요구사항이 있다. 이것은 절단오차의 정의에 의해서 나타난다. 예를 들어 2차 방법은 <math>b_1+b_2=1
</math>, <math>b_2c_2=1/2
</math>,<math>a_{21}=c_2
</math> 일 때 차수가 2이다.

일반적으로, s단계의 explicit 룽게-쿠타 방법이 p의 차수를 가진다면, <math>s\geq p
</math>이고, 만약 <math>p\geq 5
</math>일 때는 <math>s>p
</math>이다. s단계 explicit 룽게-쿠타 방법이 열린문제에서 p의 차수를 가지기 위해서는 최소한의 s가 필요하다. 몇 개의 알려진 값들은 다음과 같다.

:<math>\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}
</math>
===예시===
RK4는 이 방법에 맞아들어가며, 그 테이블은 다음과 같다.

:<math>\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}
</math>

이 룽게-쿠타 방법에서의 약간의 변형을 준 방법은 쿠타에 의해서 연구되었고 3/8법칙 이라고 불린다. 이 방법의 가장 큰 장점은 다른 일반적인 방법들보다 오차계수가 작은 것이나, 이 방법은 매 단계마다 조금 더 많은 FLOP([[부동소수점]] 연산)들이 필요하다. 그 Butcher 테이블은 다음과 같다.

:<math>\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}</math>

그러나, 가장 간단한 룽게-쿠타 방법은 (전방)[[오일러 방법]]이다. 공식은 <math>y_{n+1}=y_n+hf(t_n,y_n)</math>이다. 이것은 오직 한 단계 뿐인 explicit 룽게-쿠타 방법이다. 대응되는 테이블은 다음과 같다.

:<math>\begin{array}{c|cc} 0 \\ \hline  &1\end{array}</math>

===두 단계의 이차 방법들===
간단한 예로 두 단계의 2차 방법인 [[중간점 방법]]을 보면 다음과 같다.

:<math>y_{n+1}=y_n+hf\biggl(t_n+\frac12h,y_n+\frac12hf(t_n,y_n)\biggr).</math>

대응되는 테이블은 다음과 같다.

:<math>\begin{array}{c|cc} 0\\ 1/2&1/2\\\hline&0&1\\ \end{array}</math>

두 단계의 2차 룽게-쿠타 방법은 중간점 방법 밖에 없는 것은 아니다. 이런 종류의 방법은 α에 의해 다음과 같이 매개변수화되어 있다.<math>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)</math>

그 Butcher 테이블은 다음과 같다.

:<math>\begin{array}{c|cc}0\\\alpha&\alpha\\ \hline&(1-\frac1{2\alpha})&\frac1{2\alpha}\\ \end{array}</math>

이 때 <math>\alpha=\frac12</math>일 때는 중간점 방법이 되며, <math>\alpha=1</math>일 때는 [[호인의 방법]]이 된다.

== 활용 ==
예를 들어 α=2/3인 두 단계의 2차 룽게 쿠타 방법을 보자, 이것은 라슬톤 방법 이라고 알려져있다. 테이블은 다음과 같이 주어져 있다.

:<math>\begin{array}{c|cc}0 \\ 2/3&2/3\\\hline&1/4&3/4\\ \end{array}</math>

대응되는 등식은  다음과 같다.

:<math>\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}</math>

다음의 초기치 문제를 풀기 위하여 이 방법을 사용한다.

:<math>y^\prime=\tan(y)+1,\quad y_0=1,t\in[1,1.1]</math>

''h'' = 0.025의 크기로 구간을 잡으면 이 방법은 네 단계가 필요하다.

이 방법의 결과는 다음과 같다:

:<math>\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}</math>
수치해석적 결과는 밑줄친 값들이다.

== 적응적(Adaptive) 룽게 - 쿠타 방법 ==
Adaptive한 방법은 단일 룽게-쿠타 방법의 지역적인 절단 오차의 추정치를 찾기 위해 고안되었다. 이것은 테이블에서 차수가 <math>p</math>와 <math>p-1</math>인 두 방법들을 가져서 이루어졌다.

낮은 차수의 단계는 다음과 같이 주어진다.

:<math> y^*_{n+1}=y_n+h\sum_{i=1}^sb^*_ik_i,</math>

<math> k_i</math>가 높은 차원의 방법에서도 같다면 오차는 다음과 같다.

:<math> e_{n+1}=y_{n+1}-y^*_{n+1}=h\sum_{i=1}^s(b_1-b^*_i)k_i,</math>

오차는 <math> O(h^p)</math>이다. 이런 종류의 Butcher 테이블은 <math> b^8_i</math>의 값을 얻기 위하여 확장 되었다.

:<math> \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}</math>

[[룽게-쿠타-펠베르크 방법]]은 차수가 5와 4인 두 방법을 가지고 있다. 그 확장된 [[부처 태블로]]는 다음과 같다:

:<math> \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}</math>

하지만 가장 간단한 adaptive 룽게-쿠타 방법은 2차인 [[호인의 방법]]과 1차인 [[오일러 방법]]과 결합되어 있다. 그 확장된 Butcher 테이블은 다음과 같다:

:<math> \begin{array}{c|cc}0\\1&1\\\hline&1/2&1/2\\&1&0\\ \end{array}</math>

오차 추정은 단계 크기를 제어하기 위해 사용된다.

다른 adaptive 룽게-쿠타 방법은 [[Bogacki–Shampine 방법]](3과 2의 차수), [[Cash–Karp 방법]]과 [[Dormand–Prince 방법]](둘 다 5,4의 차수)이 있다.

== Implicit 룽게 - 쿠타 방법 ==
위에서 언급된 모든 룽게-쿠타 방법들은 모두 explict 방법들이다. explict 룽게-쿠타 방법은 일반적으로 [[:en:Stiff_equation|stiff equation]]에는 맞지 않는다. 왜냐하면 그들의 절대적인 안정성의 범위가 작기 때문이다; 특히, 그 식들은 한정되어 있다. 이 문제는 편미분방정식의 솔루션에서 특히 중요하다.

explict 룽게-쿠타 방법의 불안정성은 음함수적(암시적, [[:en:Explicit_and_implicit_methods|implicit]])한 방법을 만들게 되었다. implicit 룽게-쿠타 방법은 다음의 형태를 가진다.

:<math> y_{n+1}=y_n+h\sum_{i=1}^sb_ik_i</math>

이 때,

:<math> k_i=f\Biggl(t_n+c_ih,y_n+h\sum_{j=1}^sa_{ij}k_j\Biggr),\quad i=1,\dots,s.</math>

explict 방법간의 차이는 explict 방법은 시그마의 위에 ''i'' - 1이 들어간다는 것이다. 이것은 Butcher 테이블에서도 나타난다: explict 방법의 계수들의 행렬은 하 삼각 행렬이다. implicit 방법은 시그마 위에 ''s''가 올라가고 계수들의 행렬은 다음과 같이 삼각행렬이 아니다.

:<math> \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}</math>

이 차이는 매 단계마다, 대수적인 방정식의 시스템을 풀어야 한다는 점이다. 이것은 계산비용을 상당히 증가시킨다. ''m''개의 성분이 있는 미분방정식을 ''s''단계의 방법을 사용하면 대수적인 방정식은 ''ms''개의 성분이 있다.이것은 implict한 [[선형 다단계 방법]](ODE들의 다른 방법)과 대조된다: implicit한 ''s''단계의 선형 다단계 방법은 ''m''개의 성분을 가진 방정식을 풀기 위해서 단계수가 늘어나도 시스템의 크기가 늘어나지 않는다

=== 예시 ===
가장 간단한 implicit 룽게-쿠타 방법은 [[역 오일러 방법]]이다:

:<math> y_{n+1}=y_n+hf(t_n+h,y_{n+1}). </math>

Butcher 테이블은 다음과 같다.

:<math> \begin{array}{c|c}1&1\\\hline&1\\\end{array} </math>

이 Butcher 테이블에 대응하는 식은 다음과 같다.

:<math> k_1=f(t_n+h,y_n+hk_1)\quad and\quad y_{n+1}=y_n+hk_1 </math>

이것은 위에서 나온 역 오일러 방법의 형태로 바꿀 수 있다.

다른 implicit 룽게-쿠타 방법의 예는 [[사다리꼴 법칙]]이다. 그 Butcher 테이블은 다음과 같다.

:<math> \begin{array}{c|cc}0&0&0\\1&\frac12&\frac12\\\hline&\frac12&\frac12\\&1&0\\\end{array} </math><math>  </math>

사다리꼴 법칙은 [[배열 방법]]이다. 모든 배열 방법들은 모두 implicit 룽게-쿠타 방법이지만 모든 implicit룽게-쿠타 방법은 배열 방법이 아니다.
[[Gauss-Legendre 방법]]은 [[가우스 구적법]]에 기반한 배열 방법이다. ''s''단계의 Gauss-Legendre 방법은 차수가 2''s''이다(따라서 임의의 높은 차수의 방법들이 구성될 수 있다).두 단계의 방법(차수는 4이다)의 Buther 테이블은 다음과 같다:

:<math> \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} </math>

===안정도===
explicit한 룽게-쿠타 방법에 비교하여 implicit 방법의 장점은 안정도가 특히 [[딱딱한 방정식]]에서 높다는 것이다. 다음의 선형 방정식 ''y'''=λ''y''를 고려해 보자. 룽게-쿠타 방법을 가용하면 이 식은 다음<math>y_{n+1}=r(h\lambda)y_n</math>의 반복으로 줄며, ''r''은 다음과 같다.

:<math>r(z)=1+zb^T(I-zA)^{-1}e=\frac{\det(I-zA+zeb^T)}{\det(I-zA)},</math>

여기서 ''e''는 1들의 벡터를 나타낸다. 함수 ''r''은 ''안정도 함수''라고 불린다. 만약 이 방법이 ''s''단계라면 ''r''은 차수 ''s''의 두 다항식의 몫의 형태를 띈다. explicit 방법은 순상삼각행렬인 행렬 ''A''를 가진다. 이것은 det(''I''-''zA'') = 1이라는 것을 의미하며 그 안정도 함수는 다항함수이다.

위의 선형방정식은 <math>|r(z)|<1</math>이고 <math>z=h\lambda</math>일 때, 수치적 해법은 0이 된다. 이런''z''들의 집합은 ''절대 안정성의 영역''이라고 부른다. 특히 이 방법은 모든 실수부가 음인 z가 절대 안정성의 영역 안에 있을 때 [[A-안정성|A-안정]]하다고 불린다. explicit 룽게-쿠타 방법의 안정도 함수는 다항함수이기 때문에 explicit 룽게 쿠타 방법은 A-안정할 수 없다.

== 알고리즘 ==

== 코드 ==
* C/C++ 코드<syntaxhighlight lang="cpp">
//C/C++ code of 4th order Runge-Kutta method

# include<stdio.h>
# include<conio.h>
# include<math.h>

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<at);
  getch();
}
</syntaxhighlight>
* Python<syntaxhighlight lang="python" line="1">
# 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()
</syntaxhighlight>

== 외부 링크 ==
* {{springer|title=Runge-Kutta method|id=p/r082810}}
* [http://numericalmethods.eng.usf.edu/topics/runge_kutta_4th_method.html Runge–Kutta 4th-Order Method]
* [https://github.com/USNavalResearchLaboratory/TrackerComponentLibrary/tree/master/Mathematical_Functions/Differential_Equations Tracker Component Library Implementation in Matlab] — Implements 32 embedded Runge Kutta algorithms in <code>RungeKStep</code>, 24 embedded Runge-Kutta Nyström algorithms in <code>RungeKNystroemSStep</code> and 4 general Runge-Kutta Nyström algorithms in <code>RungeKNystroemGStep</code>.

{{수치상미분방정식}}

[[분류:수치해석학]]
[[분류:룽게-쿠타 방법]]