룽게-쿠타 방법

미분방정식
나비에-스토크스 방정식은 장애물 주위의 공기흐름을 시뮬레이션하는데 사용된다.
주제 자연과학 공학 천문학 물리학 화학 생물학 지질학 응용수학 연속체역학 혼돈 이론 동역학계 사회과학 경제학 개체군역학
분류
종류 상미분방정식 편미분방정식 미분대수방정식 적분미분방정식 분수계 미분방정식 선형미분방정식 비선형미분방정식 변수의 종류 독립변수와 종속변수 자율미분방정식 복소미분방정식 연립미분방정식 비연립미분방정식 완전미분방정식 동차미분방정식 비동차미분방정식 특질 차수 미분연산자 미분기호
관련 내용 차분방정식 확률미분방정식 지연미분방정식
풀이
일반론 피카르-린델뢰프 정리 론스키 행렬식 위상모방 위상 공간 랴푸노프 안정성 점근적 안정성 지수적 안정성 수렴속도 미분방정식의 멱급수 해법 미분방정식의 적분 해법 수치적분 디랙 델타 함수
풀이법 검산 특성곡선법 오일러 방법 지수응답공식 유한차분법 (크랭크-니콜슨 방법) 유한요소법 무한요소법 유한체적법 갤러킨 방법 페트로프-갤러킨 방법 적분인자 적분 변환 섭동 이론 룽게-쿠타 방법 변수분리법 미정계수법
인물 아이작 뉴턴 레온하르트 오일러 에밀 피카르 유제프 마리아 호에네브론스키 에른스트 레오나르드 린델뢰프 루돌프 립시츠 오귀스탱 루이 코시 존 크랭크 필리스 니콜슨 카를 다비트 톨메 룽게 마르틴 빌헬름 쿠타
v t e

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

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

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

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

${\displaystyle y'=f(t,y),\quad y(t_{0})=y_{0}.}$

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

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

${\displaystyle {\begin{aligned}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{aligned}}}$

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

${\displaystyle {\begin{aligned}k_{1}&=f(t_{n},y_{n})\\k_{2}&=f(t_{n}+{\tfrac {1}{2}}h,y_{n}+{\tfrac {1}{2}}hk_{1})\\k_{3}&=f(t_{n}+{\tfrac {1}{2}}h,y_{n}+{\tfrac {1}{2}}hk_{2})\\k_{4}&=f(t_{n}+h,y_{n}+hk_{3})\end{aligned}}}$

여기서 ${\displaystyle {y}_{n+1}}$이란 ${\displaystyle y(t_{n+1})}$ 의 RK4 근사 항이며, 다음 단계의 값(${\displaystyle {y}_{n+1}}$)은 현재 단계의 값(${\displaystyle y_{n}}$)에 네 구간의 크기 h의 증분치들의 가중평균을 더한 값이다.

• ${\displaystyle k_{1}}$은 구간의 시작부분의 기울기에 의한 증분값이며, ${\displaystyle y}$를 사용한다.(오일러 방법)
• ${\displaystyle k_{2}}$는 구간의 중간의 기울기에 의한 증분값이며, ${\displaystyle y+{\frac {h}{2}}k_{1}}$을 사용한다.
• ${\displaystyle k_{3}}$는 마찬가지로 구간의 중간의 기울기에 의한 증분값이나, ${\displaystyle y+{\frac {h}{2}}k_{2}}$를 사용한다.
• ${\displaystyle k_{4}}$은 구간의 끝의 기울기에 의한 증분값이며, ${\displaystyle y+hk_{3}}$를 사용한다.

네 증분을 평균을 낼 때, 중간 부분에서 더 많은 무게의 증분이 더해진다. 만약 함수 ${\displaystyle f}$가 ${\displaystyle y}$에 대하여 독립적이라면, 미분방정식은 단순한 적분이 되어 RK4는 심프슨 공식이 된다.

RK4 방법은 4차 해법이다. 이는 각 단계에서의 오류는 점근 표기법으로 ${\displaystyle O(h^{5})}$, 총 누적 오류는 ${\displaystyle O(h^{4})}$라는 것을 의미한다.

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

${\displaystyle y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i}}$

이때 각 ${\displaystyle k_{i}}$는 다음과 같다.

${\displaystyle {\begin{aligned}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{aligned}}}$

특정한 방법을 지정하기 위해서는, 정수 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)이라는 기억장치 내에서 배열되어있다. 존 찰스 부처의 이름을 딴 것이다.

${\displaystyle {\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}}}$

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

${\displaystyle \sum _{j=1}^{i-1}a_{ij}=c_{i}\ for\ i=2,\dots ,s.}$

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

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

${\displaystyle {\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는 이 방법에 맞아들어가며, 그 테이블은 다음과 같다.

${\displaystyle {\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 테이블은 다음과 같다.

${\displaystyle {\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}}}$

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

${\displaystyle {\begin{array}{c|cc}0\\\hline &1\end{array}}}$

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

${\displaystyle y_{n+1}=y_{n}+hf{\biggl (}t_{n}+{\frac {1}{2}}h,y_{n}+{\frac {1}{2}}hf(t_{n},y_{n}){\biggr )}.}$

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

${\displaystyle {\begin{array}{c|cc}0\\1/2&1/2\\\hline &0&1\\\end{array}}}$

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

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

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

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

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

${\displaystyle {\begin{array}{c|cc}0\\2/3&2/3\\\hline &1/4&3/4\\\end{array}}}$

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

${\displaystyle {\begin{aligned}k_{1}&=f(t_{n},y_{n}),\\k_{2}&=f(t_{n}+{\frac {2}{3}}h,y_{n}+{\frac {2}{3}}hk_{1}),\\y_{n+1}&=y_{n}+h{\Bigl (}{\frac {1}{4}}k_{1}+{\frac {3}{4}}k_{2}{\Bigr )}.\end{aligned}}}$

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

${\displaystyle y^{\prime }=\tan(y)+1,\quad y_{0}=1,t\in [1,1.1]}$

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

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

${\displaystyle {\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}+{\frac {2}{3}}h,y_{0}+{\frac {2}{3}}hk_{1})=2.713898184\\&y_{1}=y_{0}+h({\frac {1}{4}}k_{1}+{\frac {3}{4}}k_{2})={\underline {1.066869388}}\\t_{2}=1.05:\\&y_{1}=1.066869388\quad \ k_{1}=2.813524695\quad k_{2}=f(t_{1}+{\frac {2}{3}}h,y_{1}+{\frac {2}{3}}hk_{1})\\&y_{2}=y_{1}+h({\frac {1}{4}}k_{1}+{\frac {3}{4}}k_{2})={\underline {1.141332181}}\\t_{3}=1.075:\\&y_{2}=1.141332181\quad \ k_{1}=3.183536647\quad k_{2}=f(t_{2}+{\frac {2}{3}}h,y_{2}+{\frac {2}{3}}hk_{1})\\&y_{3}=y_{2}+h({\frac {1}{4}}k_{1}+{\frac {3}{4}}k_{2})={\underline {1.227417567}}\\t_{4}=1.1:\\&y_{3}=1.227417567\quad \ k_{1}=3.796866512\quad k_{2}=f(t_{3}+{\frac {2}{3}}h,y_{3}+{\frac {2}{3}}hk_{1})\\&y_{4}=y_{3}+h({\frac {1}{4}}k_{1}+{\frac {3}{4}}k_{2})={\underline {1.335079087}}.\\\end{array}}}$

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

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

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

${\displaystyle y_{n+1}^{*}=y_{n}+h\sum _{i=1}^{s}b_{i}^{*}k_{i},}$

${\displaystyle k_{i}}$가 높은 차원의 방법에서도 같다면 오차는 다음과 같다.

${\displaystyle e_{n+1}=y_{n+1}-y_{n+1}^{*}=h\sum _{i=1}^{s}(b_{1}-b_{i}^{*})k_{i},}$

오차는 ${\displaystyle O(h^{p})}$이다. 이런 종류의 Butcher 테이블은 ${\displaystyle b_{i}^{8}}$의 값을 얻기 위하여 확장 되었다.

${\displaystyle {\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인 두 방법을 가지고 있다. 그 확장된 부처 태블로는 다음과 같다:

${\displaystyle {\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 테이블은 다음과 같다:

${\displaystyle {\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의 차수)이 있다.

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

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

${\displaystyle y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i}}$

이 때,

${\displaystyle k_{i}=f{\Biggl (}t_{n}+c_{i}h,y_{n}+h\sum _{j=1}^{s}a_{ij}k_{j}{\Biggr )},\quad i=1,\dots ,s.}$

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

${\displaystyle {\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 룽게-쿠타 방법은 역 오일러 방법이다:

${\displaystyle y_{n+1}=y_{n}+hf(t_{n}+h,y_{n+1}).}$

Butcher 테이블은 다음과 같다.

${\displaystyle {\begin{array}{c|c}1&1\\\hline &1\\\end{array}}}$

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

${\displaystyle k_{1}=f(t_{n}+h,y_{n}+hk_{1})\quad and\quad y_{n+1}=y_{n}+hk_{1}}$

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

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

${\displaystyle {\begin{array}{c|cc}0&0&0\\1&{\frac {1}{2}}&{\frac {1}{2}}\\\hline &{\frac {1}{2}}&{\frac {1}{2}}\\&1&0\\\end{array}}}$${\displaystyle }$

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

${\displaystyle {\begin{array}{c|cc}{\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}}\\\hline &{\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를 고려해 보자. 룽게-쿠타 방법을 가용하면 이 식은 다음${\displaystyle y_{n+1}=r(h\lambda )y_{n}}$의 반복으로 줄며, r은 다음과 같다.

${\displaystyle 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이라는 것을 의미하며 그 안정도 함수는 다항함수이다.

위의 선형방정식은 ${\displaystyle |r(z)|<1}$이고 ${\displaystyle 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 main():
      t0,y0=input("initial time:"),input("initial value:")
      h=input("interval size:")
      t=input("terminal time:")
      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()
  

• “Runge-Kutta method”. 《Encyclopedia of Mathematics》 (영어). Springer-Verlag. 2001. ISBN 978-1-55608-010-4. 
• 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.