룽게-쿠타 방법

틀:미분방정식 사이드바 수치 해석에서 룽게-쿠타 방법(Runge-Kutta方法, 틀:Llang)은 적분 방정식 중 초기값 문제를 푸는 방법 중 하나이다.

역사

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

일반적인 4차 룽게-쿠타 방법

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

다음과 같이 정의된 초기값 문제를 두자:

$y^{\prime }=f(t,y),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{array}{cc}{y}_{n+1}& =y_n+\frac{1}{6}h(k_1+2k_2+2k_3+k_4)\\ t_{n+1}& =t_n+h\\ \end{array}$

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

$\begin{array}{cc}{k}_1& =f(t_n,y_n)\\ k_2& =f(t_n+\frac{1}{2}h,y_n+\frac{1}{2}h{k}_1)\\ k_3& =f(t_n+\frac{1}{2}h,y_n+\frac{1}{2}h{k}_2)\\ k_4& =f(t_n+h,y_n+h{k}_3)\end{array}$

여기서 $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{array}{cc}{k}_1& =f(t_n,y_n),\\ k_2& =f(t_n+c_{2}h,y_n+h(a_{21}{k}_1)),\\ k_3& =f(t_n+c_{3}h,y_n+h(a_{31}{k}_1+a_{32}{k}_2)),\\ & \cdots \\ k_s& =f(t_n+c_{s}h,y_n+h(a_{s1}{k}_1+a_{s2}{k}_2+\cdots +a_{s,s-1}{k}_{s-1})).\end{array}$

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

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

$\begin{array}{c}0\\ c_2& a_{21}\\ c_3& a_{31}& a_{32}\\ ⋮& ⋮& & \ddots \\ c_s& a_{s1}& a_{s2}& \cdots & a_{s,s-1}\\ & b_1& b_2& \cdots & b_{s-1}& b_s\end{array}$

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

$\sum _{j=1}^{i-1}{a}_{ij}=c_{i}fori=2,\dots ,s.$

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

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

$\begin{array}{ccccccccc}p& 1& 2& 3& 4& 5& 6& 7& 8\\ \min s& 1& 2& 3& 4& 6& 7& 9& 11\end{array}$

예시

RK4는 이 방법에 맞아들어가며, 그 테이블은 다음과 같다.

$\begin{array}{c}0\\ 1/2& 1/2\\ 1/2& 0& 1/2\\ 1& 0& 0& 1\\ & 1/6& 1/3& 1/3& 1/6\end{array}$

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

$\begin{array}{c}0\\ 1/3& 1/3\\ 2/3& -1/3& 1\\ 1& 1& -1& 1\\ & 1/8& 3/8& 3/8& 1/8\end{array}$

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

$\begin{array}{c}0\\ & 1\end{array}$

두 단계의 이차 방법들

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

$y_{n+1}=y_n+hf(t_n+\frac{1}{2}h,y_n+\frac{1}{2}hf(t_n,y_n)).$

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

$\begin{array}{c}0\\ 1/2& 1/2\\ & 0& 1\end{array}$

두 단계의 2차 룽게-쿠타 방법은 중간점 방법 밖에 없는 것은 아니다. 이런 종류의 방법은 α에 의해 다음과 같이 매개변수화되어 있다.$y_{n+1}=y_n+h\left(\right(1-\frac{1}{2\alpha })f(t_n,y_n)+\frac{1}{2\alpha }f(t_n+\alpha h,y_n+\alpha hf(t_n,y_n)\left)\right)$

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

$\begin{array}{c}0\\ \alpha & \alpha \\ & (1-\frac{1}{2\alpha })& \frac{1}{2\alpha }\end{array}$

이 때 $\alpha =\frac{1}{2}$일 때는 중간점 방법이 되며, $\alpha =1$일 때는 호인의 방법이 된다.

활용

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

$\begin{array}{c}0\\ 2/3& 2/3\\ & 1/4& 3/4\end{array}$

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

$\begin{array}{cc}{k}_1& =f(t_n,y_n),\\ k_2& =f(t_n+\frac{2}{3}h,y_n+\frac{2}{3}h{k}_1),\\ y_{n+1}& =y_n+h(\frac{1}{4}{k}_1+\frac{3}{4}{k}_2).\end{array}$

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

$y^{\prime }=\tan (y)+1,y_0=1,t\in [1,1.1]$

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

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

$\begin{array}{c}{t}_0=1:\\ & y_0=1\\ t_1=1.025:\\ & y_0=1k_1=2.557407725k_2=f(t_0+\frac{2}{3}h,y_0+\frac{2}{3}h{k}_1)=2.713898184\\ & y_1=y_0+h(\frac{1}{4}{k}_1+\frac{3}{4}{k}_2)=\underset{\_}{1.066869388}\\ t_2=1.05:\\ & y_1=1.066869388k_1=2.813524695k_2=f(t_1+\frac{2}{3}h,y_1+\frac{2}{3}h{k}_1)\\ & y_2=y_1+h(\frac{1}{4}{k}_1+\frac{3}{4}{k}_2)=\underset{\_}{1.141332181}\\ t_3=1.075:\\ & y_2=1.141332181k_1=3.183536647k_2=f(t_2+\frac{2}{3}h,y_2+\frac{2}{3}h{k}_1)\\ & y_3=y_2+h(\frac{1}{4}{k}_1+\frac{3}{4}{k}_2)=\underset{\_}{1.227417567}\\ t_4=1.1:\\ & y_3=1.227417567k_1=3.796866512k_2=f(t_3+\frac{2}{3}h,y_3+\frac{2}{3}h{k}_1)\\ & y_4=y_3+h(\frac{1}{4}{k}_1+\frac{3}{4}{k}_2)=\underset{\_}{1.335079087}.\end{array}$

수치해석적 결과는 밑줄친 값들이다.

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

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

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

$y_{n+1}^{*}=y_n+h\sum _{i=1}^s{b}_i^{*}{k}_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_i^8$의 값을 얻기 위하여 확장 되었다.

$\begin{array}{c}0\\ c_2& a_{21}\\ c_3& a_{31}& a_{32}\\ ⋮& ⋮& & \ddots \\ c_s& a_{s1}& a_{s2}& \cdots & a_{s,s-1}\\ & 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}{c}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\\ & 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}0\\ 1& 1\\ & 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}^s{b}_i{k}_i$

이 때,

$k_i=f(t_n+c_{i}h,y_n+h\sum _{j=1}^s{a}_{ij}{k}_j),i=1,\dots ,s.$

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

$\begin{array}{ccccc}{c}_1& a_{11}& a_{12}& \cdots & a_{1s}\\ c_2& a_{21}& a_{22}& \cdots & a_{2s}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ c_s& a_{s1}& a_{s2}& \cdots & a_{ss}\\ & b_1& b_2& \cdots & b_s\\ & b_1^{*}& b_2^{*}& \cdots & b_s^{*}\end{array}=\begin{array}{cc}c& A\\ & 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}{cc}1& 1\\ & 1\end{array}$

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

$k_1=f(t_n+h,y_n+h{k}_1)and{y}_{n+1}=y_n+h{k}_1$

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

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

$\begin{array}{ccc}0& 0& 0\\ 1& \frac{1}{2}& \frac{1}{2}\\ & \frac{1}{2}& \frac{1}{2}\\ & 1& 0\end{array}$${\displaystyle }$

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

$\begin{array}{ccc}\frac{1}{2}-\frac{1}{6}\sqrt{3}& \frac{1}{4}& \frac{1}{4}-\frac{1}{6}\sqrt{3}\\ \frac{1}{2}+\frac{1}{6}\sqrt{3}& \frac{1}{4}+\frac{1}{6}\sqrt{3}& \frac{1}{4}\\ & \frac{1}{2}& \frac{1}{2}\\ & \frac{1}{2}+\frac{1}{2}\sqrt{3}& \frac{1}{2}-\frac{1}{2}\sqrt{3}\end{array}$

안정도

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

$r(z)=1+z{b}^T(I-zA{)}^{-1}e=\frac{\det (I-zA+ze{b}^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<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();
  }
  

• 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()
  

외부 링크

• 틀:Springer
• Runge–Kutta 4th-Order Method
• Tracker Component Library Implementation in Matlab — Implements 32 embedded Runge Kutta algorithms in RungeKStep, 24 embedded Runge-Kutta Nyström algorithms in RungeKNystroemSStep and 4 general Runge-Kutta Nyström algorithms in RungeKNystroemGStep.