룽게-쿠타 방법

3. 룽게-쿠타 방법의 종류

일반적으로 사용되는 룽게-쿠타 방법은 4차 룽게-쿠타 방법(RK4)이며, 특별한 언급이 없으면 대체로 이 방법을 의미한다.

룽게-쿠타 방법은 크게 양함수적(명시적) 방법과 음함수적(암시적) 방법으로 나뉜다.

• 양함수적(Explicit) 룽게-쿠타 방법: RK4 방법을 일반화한 형태로, 미래 값을 계산하기 위해 현재 및 과거 단계의 값을 명시적으로 사용한다.
• 음함수적(Implicit) 룽게-쿠타 방법: 강성 방정식에 적합하며, 미래 값을 계산하기 위해 현재 및 미래 단계의 값을 암시적으로 사용한다. 안정성이 높지만, 매 단계마다 대수 방정식을 풀어야 하므로 계산 비용이 더 크다.

3. 1. 고전적 룽게-쿠타 방법 (RK4)

• $k\_1$은 구간 시작점의 기울기 (오일러 방법)
• $k\_2$는 구간 중간점의 기울기 ($y$와 $k\_1$ 사용)
• $k\_3$는 다시 중간점의 기울기 ($y$와 $k\_2$ 사용)
• $k\_4$는 구간 끝점의 기울기 ($y$와 $k\_3$ 사용)

3. 2. 양함수적(명시적) 룽게-쿠타 방법

3. 2. 1. 예시

2단계 2차 방법의 예시로는 중간점 방법이 있다. 공식은 다음과 같다.

:$y\_\{n+1\}\; =\; y\_n\; +\; hf\backslash left(t\_n+\backslash frac\{1\}\{2\}h,\; y\_n+\backslash frac\{1\}\{2\}hf(t\_n,\backslash \; y\_n)\backslash right).$

해당하는 부처 태블로는 다음과 같다.

중간점 방법은 2단계 2차 룽게-쿠타 방법의 유일한 형태는 아니다. 다음과 같은 공식으로 α로 매개변수화되는 방법들이 있다.

:$y\_\{n+1\}\; =\; y\_n\; +\; h\backslash bigl(\; (1-\backslash tfrac1\{2\backslash alpha\})\; f(t\_n,\; y\_n)\; +\; \backslash tfrac1\{2\backslash alpha\}\; f(t\_n\; +\; \backslash alpha\; h,\; y\_n\; +\; \backslash alpha\; h\; f(t\_n,\; y\_n))\; \backslash bigr).$

이에 대한 부처 태블로는 다음과 같다.

이때, $\backslash alpha=\backslash tfrac12$는 중간점 방법을, $\backslash alpha=1$는 호인의 방법을,^[5] $\backslash alpha=\backslash tfrac23$는 Ralston's method를 나타낸다.

α = 2/3인 2단계 2차 룽게-쿠타 방법은 Ralston 방법으로도 알려져 있으며, 부처 태블로는 다음과 같다.

해당하는 방정식은 다음과 같다.

:$\backslash begin\{align\}$

k_1 &= f(t_n,\ y_n), \\

k_2 &= f(t_n + \tfrac{2}{3}h,\ y_n + \tfrac{2}{3}h k_1), \\

y_{n+1} &= y_n + h\left(\tfrac{1}{4}k_1+\tfrac{3}{4}k_2\right).

\end{align}

이 방법은 초기값 문제 $\backslash frac\{dy\}\{dt\}\; =\; \backslash tan(y)+1,\backslash quad\; y\_0=1,\backslash \; t\backslash in\; [1,\; 1.1]$를 해결하는 데 사용될 수 있다. 단계 크기 ''h'' = 0.025를 사용하면, 네 단계를 거쳐야 한다.

계산 과정은 다음과 같다.

수치 해는 밑줄 친 값에 해당한다.

RK4로 불리는 4차 룽게-쿠타 방법의 부처 태블로는 다음과 같다.^[18]

쿠타의 3/8 법칙의 부처 태블로는 다음과 같다.^[19]

3. 3. 적응적(Adaptive) 룽게-쿠타 방법

3. 4. 음함수적(Implicit) 룽게-쿠타 방법

$b^*$ 행은 적응형 룽게-쿠타 방법에 해당한다.

이러한 차이로 인해 음함수적 방법은 매 단계마다 대수 방정식 시스템을 풀어야 한다.^[27] 이는 계산 비용을 상당히 증가시키는데, ''s'' 단계 방법을 사용하여 ''m''개의 구성 요소가 있는 미분 방정식을 풀 경우, 대수 방정식 시스템은 ''ms''개의 구성 요소를 갖는다. 반면, 음적 선형 다단계 방법은 단계 수가 증가해도 시스템의 크기가 증가하지 않는다.^[27]

4. 룽게-쿠타-뉘스트룀 방법

룽게-쿠타-뉘스트룀 방법은 2차 미분 방정식을 풀기 위해 최적화된 특수한 룽게-쿠타 방법이다.^[22][23]

차수가 $s$인 2차 상미분 방정식(ODE) 시스템에 대한 일반적인 룽게-쿠타-뉘스트룀 방법의 형태는 다음과 같다.

$$$$

\ddot y_i = f_i (y_1, y_2, \ldots, y_n )



$$\backslash begin\{cases\}$$

g_i = y_m + c_i h {\dot y}_m + h^2 \sum_{j=1}^s a_{ij} f(g_j), & i=1,2,\ldots, s\\

y_{m+1} = y_{m} + h {\dot y}_m + h^2 \sum_{j=1}^s \bar{b}_j f(g_j) \\

{\dot y}_{m+1} = {\dot y}_m + h \sum_{j=1}^s b_j f(g_j)

\end{cases}

이는 다음과 같은 형태의 부처 테이블(Butcher tableau)로 나타낼 수 있다.

$$$$

\begin{array}{c|cccc}

c_1 & a_{11} & a_{12}& \dots & a_{1s}\\

c_2 & a_{21} & a_{22}& \dots & a_{2s}\\

\vdots & \vdots & \vdots& \ddots& \vdots\\

c_s & a_{s1} & a_{s2}& \dots & a_{ss} \\

\hline

& {\bar b}_1 & {\bar b}_2 & \dots & {\bar b}_s\\

& b_1 & b_2 & \dots & b_s

\end{array} =

\begin{array}{c|c}

\mathbf{c}& \mathbf{A}\\

\hline

& \mathbf^\top \\

& \mathbf{b}^\top

\end{array}



다음은 두 가지 4차 명시적 RKN 방법의 부처 테이블 예시이다.

이 두 방식은 원래 방정식이 보존적 고전 역학 시스템에서 파생된 경우, 즉

$$$$

f_i(x_1, \ldots, x_n ) = \frac{\partial V}{\partial x_i} (x_1, \ldots, x_n )



어떤 스칼라 함수 $V$에 대해, 심플렉틱 보존(symplectic preserving) 속성을 갖는다.^[24]

5. 안정성

Runge–Kutta^영어 방법의 안정성은 A-안정성과 B-안정성으로 설명할 수 있다. 음함수적 방법은 양함수적 방법에 비해 안정성이 높다.

지금까지 언급된 모든 룽게-쿠타 방법은 명시적 방법이다. 명시적 룽게-쿠타 방법은 일반적으로 강성 방정식의 해를 구하는 데 적합하지 않은데, 절대 안정 영역이 작기 때문이다.^[25] 특히 그 영역은 제한적이다. 이 문제는 특히 수치 편미분 방정식의 해를 구할 때 중요하다.

명시적 룽게-쿠타 방법의 불안정성은 음적 방법의 개발을 촉진했다. 묵시적 룽게-쿠타 방법의 가장 간단한 예는 후진 오일러 방법이다.

:$y\_\{n\; +\; 1\}\; =\; y\_n\; +\; h\; f(t\_n\; +\; h,\backslash \; y\_\{n\; +\; 1\}).\; \backslash ,$

묵시적 방법은 명시적 방법에 비해 안정성이 높다.

선형 테스트 방정식 $y\text{'}\; =\; \backslash lambda\; y$를 생각해보자. 룽게-쿠타 방법을 사용하여 이 방정식에 적용하면 $y\_\{n+1\}\; =\; r(h\backslash lambda)y\_n$ 이 된다. 여기서,

: $r(z)\; =\; 1\; +\; zb^T(I-zA)^\{-1\}e\; =\; \backslash frac\{\backslash det\; (I-zA+zeb^T)\}\{\backslash det(I-zA)\}$^[40]

는 안정성 함수라고 불린다.

룽게-쿠타 방법에 의한 위의 방정식의 수치해가 0으로 감쇠하는 조건이 ($z\; =\; h\backslash lambda$) 이다. 위의 조건을 만족하는 z의 집합은 방법에 대한 '''절대 안정 영역''' (region of absolute stability) 이다. 특히, 절대 안정 영역이 왼쪽 복소수 평면을 포함할 때, 그 방법은 '''A-안정''' (A-stable) 이라고 한다. 양적 룽게-쿠타 방법은, 안정성 함수가 다항식이므로, A-안정한 방법이 될 수 없다.^[49]

s 단 가우스-르장드르 방법의 차수는 앞서 언급했듯이 2s 이다. 따라서 안정성 함수의 분자와 분모는 마찬가지로 s 차 다항식이 된다. 즉, $m\; =\; n\; =\; s$ 이다. 위의 조건이 충족되므로, 가우스-르장드르 방법은 A-안정이라는 것을 알 수 있다. 따라서 임의로 높은 차수를 갖는, A-안정한 룽게-쿠타 방법이 존재한다. 이에 비해, A-안정성을 갖는 선형 다단법의 차수는 2 이하이다.^[51]

이상의 내용으로부터, 음적 룽게-쿠타 방법은 양적 방법보다 뛰어난 안정성을 갖는다는 것을 알 수 있다.

A-안정성이라는 개념은 선형 자려 방정식 $y\text{'}\; =\; \backslash lambda\; y$를 고려한 결과이다. 일반적인 테스트 방정식 $y\text{'}\; =\; f(x,y)$와 단조성 조건 $\backslash langle\; f(x,y)\; -\; f(x,z),\; y-z\; \backslash rangle\; \backslash le\; 0$을 만족하는 f를 고려한다. 만약 모든 $h\; \backslash ge\; 0$에 대해, 부등식 $\backslash |y\_\{n+1\}\; -\; z\_\{n+1\}\backslash |\; \backslash le\; \backslash |y\_n\; -\; z\_n\backslash |$이 성립할 때, 해당 룽게-쿠타법은 '''B-안정'''이라고 한다. 여기서 $y\_n$과 $z\_n$은 각각의 초기값에 대한 수치해이다. $f(x,y)\; =\; \backslash lambda\; y$에 방법을 적용함으로써, B-안정성은 A-안정성보다 더 강한 조건임을 알 수 있다.

6. 활용

룽게-쿠타 방법은 과학기술 분야에서 시뮬레이션 및 수치 모델링에 널리 활용된다. 특히 공학, 물리학, 전산 유체 역학 등 다양한 분야에서 사용되고 있다.^[39]

예를 들어, 고전적 룽게-쿠타 방법(RK4)은 다음과 같은 초기값 문제를 해결하는 데 사용될 수 있다.

:$y\text{'}\; =\; f(t,\; y),\; \backslash quad\; y(t\_0)\; =\; y\_0.$

여기서 ''y''(''t'')는 구하고자 하는 미지 함수이며, ''t''에서의 기울기는 ''f''(''t'', ''y'')로 주어진다. 이때, RK4 방법은 다음과 같은 공식을 통해 ''y''(''t'')의 근사값을 계산한다.

:$\backslash begin\{align\}$

y_{n+1} &= y_n + {h \over 6} (k_1 + 2k_2 + 2k_3 + k_4), \\

t_{n+1} &= t_n + h.

\end{align}

:

k_2 &= f \left( t_n + {h \over 2}, y_n + {h \over 2} k_1 \right), \\

k_3 &= f \left( t_n + {h \over 2}, y_n + {h \over 2} k_2 \right), \\

k_4 &= f \left( t_n + h, y_n + hk_3 \right) \end{align}

여기서 ''h''는 시간 간격을 나타낸다. 위 식에서 ''k''₁, ''k''₂, ''k''₃, ''k''₄는 각각 구간의 시작, 중간, 끝점에서의 기울기를 나타내며, 이들의 가중 평균을 통해 다음 단계의 ''y'' 값을 예측한다.

RK4는 4차 방법으로, 테일러 전개와 비교했을 때 4차 항까지 일치하며, 한 스텝의 오차는 $O(h^5)$, 전체 오차는 $O(h^4)$ 수준이다.^[39]

이 외에도 다양한 룽게-쿠타 방법들이 존재하며, 문제의 특성과 요구되는 정확도에 따라 적절한 방법을 선택하여 활용할 수 있다. 예를 들어, 후진 오일러 방법이나 사다리꼴 공식과 같은 음적 룽게-쿠타 방법도 활용될 수 있다.

7. 알고리즘

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

다음과 같이 정의된 초기값 문제를 보자.

:$y\text{'}\; =\; f(t,\; y),\; \backslash quad\; y(t\_0)\; =\; y\_0.$

여기서 y는 시간 t에 대한 미지의 함수이며, 우리가 근사하려는 값이다. y의 변화율인 $\backslash dot\{y\}$는 t와 y 자신으로 이루어진 함수이고, 초기 시간 $t\_0$에 대응하는 y의 초기값은 $y\_0$이다. 함수 f와 $t\_0$, $y\_0$의 값은 주어져 있다.

이제 h > 0 인 단계 크기(step size)를 정하고 다음을 정의한다.

:$\backslash 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, ... 에 대해, 다음을 계산한다.^[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+\backslash frac\{h\}\{2\}k\_1$을 사용한다.
• $k\_3$는 마찬가지로 구간의 중간 지점의 기울기에 의한 증분값이나, $y+\backslash frac\{h\}\{2\}k\_2$를 사용한다.
• $k\_4$는 구간의 끝 지점의 기울기에 의한 증분값이며, $y+h\; k\_3$를 사용한다.

8. 코드

cpp



// C/C++ 4차 룽게-쿠타 방법 코드

#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 t0와 y0 값을 입력하세요:");

scanf("%f %f", &t0, &y0);

printf("\n h 값을 입력하세요:");

scanf("%f", &h);

printf("\n 최종 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();

}



```

```python



# 룽게-쿠타 4차 방법 프로그램

from math import *

def f(t, y):

return sin(y**t)

def main():

t0, y0 = input("초기 시간:"), input("초기 값:")

h = input("간격 크기:")

t = input("종료 시간:")

print("t\ty\n", '_' * 16)

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("%.4f\t%.4f" % (t1, y1))

t0, y0 = t1, y1

main()

참조

<sub>[1]</sub> 웹사이트 Runge-Kutta method https://www.dictiona[...] 2021-04-04
<sub>[2]</sub> 서적 A first course in computational physics Jones and Bartlett Publishers 2011
<sub>[3]</sub> 간행물 2007
<sub>[4]</sub> 간행물 1989
<sub>[5]</sub> 간행물 2003
<sub>[6]</sub> 간행물 2007
<sub>[7]</sub> 간행물 1996
<sub>[8]</sub> 간행물 1996
<sub>[9]</sub> 문서
<sub>[10]</sub> 간행물 2008
<sub>[11]</sub> 간행물 1965
<sub>[12]</sub> 간행물 1985
<sub>[13]</sub> 간행물 2008
<sub>[14]</sub> 간행물 1964
<sub>[15]</sub> 간행물 1970
<sub>[16]</sub> 간행물 1993
<sub>[17]</sub> 간행물 1996
<sub>[18]</sub> 간행물 2003
<sub>[19]</sub> 간행물 1993
<sub>[20]</sub> 간행물 2003
<sub>[21]</sub> 간행물 1991
<sub>[22]</sub> 논문 New Runge–Kutta Algorithms for Numerical Simulation in Dynamical Astronomy 1978-10
<sub>[23]</sub> 간행물 Classical seventh-, sixth-, and fifth-order Runge–Kutta–Nyström formulas with stepsize control for general second-order differential equations National Aeronautics and Space Administration 1974-10
<sub>[24]</sub> 논문 Canonical Runge-Kutta-Nyström (RKN) methods for second order ordinary differential equations https://dx.doi.org/1[...] 1991-01-01
<sub>[25]</sub> 간행물 2003
<sub>[26]</sub> 간행물 1996
<sub>[27]</sub> 간행물 2003
<sub>[28]</sub> 간행물 1996
<sub>[29]</sub> 간행물 1996
<sub>[30]</sub> 간행물 1996
<sub>[31]</sub> 간행물 1996
<sub>[32]</sub> 간행물 1996
<sub>[33]</sub> 간행물 1996
<sub>[34]</sub> 간행물 1996
<sub>[35]</sub> 간행물 1963
<sub>[36]</sub> 서적 1991
<sub>[37]</sub> 서적 1991
<sub>[38]</sub> 웹사이트 Appendix C. Derivation of the Numerical Integration Formulae http://www.ss.ncu.ed[...] Institute of Space Science, National Central University 2022-04-17
<sub>[39]</sub> 서적 2003
<sub>[40]</sub> 서적 2008
<sub>[41]</sub> 서적 1989
<sub>[42]</sub> 서적 2008
<sub>[43]</sub> 서적 2008
<sub>[44]</sub> 서적 2003
<sub>[45]</sub> 서적 数値解析 （共立講座 数学探検 17） 共立出版 2017
<sub>[46]</sub> 서적 1993
<sub>[47]</sub> 서적 2003
<sub>[48]</sub> 서적 2008
<sub>[49]</sub> 서적 2008
<sub>[50]</sub> 서적 2008
<sub>[51]</sub> 서적 1963