수치 해석 에서 룽게-쿠타 방법 (Runge-Kutta方法, 영어 : Runge–Kutta method )은 적분 방정식 중 초기값 문제 를 푸는 방법 중 하나이다.
1900년 경 독일의 수학자 카를 룽게 와 마르틴 빌헬름 쿠타 가 개발하였다.
일반적으로 사용하는 룽게-쿠타 방법은 4차 룽게-쿠타 방법으로 보통 "RK4" 라고 쓴다. 별다른 수식어 없이 룽게-쿠타 방법을 쓴다고 말한다면 대체로 4차 룽게-쿠타 방법을 쓴다는 뜻이다.
다음과 같이 정의된 초기값 문제 를 두자:
y
′
=
f
(
t
,
y
)
,
y
(
t
0
)
=
y
0
.
{\displaystyle y'=f(t,y),\quad y(t_{0})=y_{0}.}
y는 시간 t에 대한 미지의 함수이며, 우리가 근사하려는 것이다. y의 변화인
y
˙
{\displaystyle {\dot {y}}}
t
0
{\displaystyle t_{0}}
y
0
{\displaystyle y_{0}}
t
0
{\displaystyle t_{0}}
y
0
{\displaystyle y_{0}}
이제 h > 0 인 단계 크기로 다음을 정의한다.
y
n
+
1
=
y
n
+
1
6
h
(
k
1
+
2
k
2
+
2
k
3
+
k
4
)
t
n
+
1
=
t
n
+
h
{\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, ... 에서, 다음을 사용한다.
k
1
=
f
(
t
n
,
y
n
)
k
2
=
f
(
t
n
+
1
2
h
,
y
n
+
1
2
h
k
1
)
k
3
=
f
(
t
n
+
1
2
h
,
y
n
+
1
2
h
k
2
)
k
4
=
f
(
t
n
+
h
,
y
n
+
h
k
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}}}
여기서
y
n
+
1
{\displaystyle {y}_{n+1}}
y
(
t
n
+
1
)
{\displaystyle y(t_{n+1})}
y
n
+
1
{\displaystyle {y}_{n+1}}
y
n
{\displaystyle y_{n}}
k
1
{\displaystyle k_{1}}
y
{\displaystyle y}
오일러 방법 )
k
2
{\displaystyle k_{2}}
y
+
h
2
k
1
{\displaystyle y+{\frac {h}{2}}k_{1}}
k
3
{\displaystyle k_{3}}
y
+
h
2
k
2
{\displaystyle y+{\frac {h}{2}}k_{2}}
k
4
{\displaystyle k_{4}}
y
+
h
k
3
{\displaystyle y+hk_{3}}
네 증분을 평균을 낼 때, 중간 부분에서 더 많은 무게의 증분이 더해진다. 만약 함수
f
{\displaystyle f}
y
{\displaystyle y}
심프슨 공식 이 된다.
일반적인 미분방정식에서 RK4 방법과 다른 낮은 차원의 방법을 비교한 것이다. RK4 방법은 4차 해법이다. 이는 각 단계에서의 오류 는 점근 표기법 으로
O
(
h
5
)
{\displaystyle O(h^{5})}
O
(
h
4
)
{\displaystyle O(h^{4})}
양함수적(명시적, explicit ) 룽게-쿠타 방법은 위에서 설명된 RK4 방법의 일반화이다. 이것은 다음과 같이 주어진다.
y
n
+
1
=
y
n
+
h
∑
i
=
1
s
b
i
k
i
{\displaystyle y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i}}
이때 각
k
i
{\displaystyle k_{i}}
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
)
)
,
⋯
k
s
=
f
(
t
n
+
c
s
h
,
y
n
+
h
(
a
s
1
k
1
+
a
s
2
k
2
+
⋯
+
a
s
,
s
−
1
k
s
−
1
)
)
.
{\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 (각 단계들의 수) 가 주어져야 하고, 각각의 계수들 aij (1 ≤ i < j ≤ s ), bi (i = 1, 2, ..., s ), ci (i =2, 3, ..., s ) 가 필요하다. 행렬 [aij ]는 룽게-쿠타 행렬 이며, bi 와ci 는 알려져 있듯이 무게 와 마디 이다.
이 값들은 부처 태블로 존 찰스 부처 의 이름을 딴 것이다.
0
c
2
a
21
c
3
a
31
a
32
⋮
⋮
⋱
c
s
a
s
1
a
s
2
⋯
a
s
,
s
−
1
b
1
b
2
⋯
b
s
−
1
b
s
{\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}}}
룽게-쿠타 방법은 다음을 만족할 때 일관성이 있다.
∑
j
=
1
i
−
1
a
i
j
=
c
i
f
o
r
i
=
2
,
…
,
s
.
{\displaystyle \sum _{j=1}^{i-1}a_{ij}=c_{i}\ for\ i=2,\dots ,s.}
만약 하나가 특정 차수 p ,즉 지역 절단 오차가
O
(
h
p
+
1
)
{\displaystyle O(h^{p+1})}
b
1
+
b
2
=
1
{\displaystyle b_{1}+b_{2}=1}
b
2
c
2
=
1
/
2
{\displaystyle b_{2}c_{2}=1/2}
a
21
=
c
2
{\displaystyle a_{21}=c_{2}}
일반적으로, s단계의 explicit 룽게-쿠타 방법이 p의 차수를 가진다면,
s
≥
p
{\displaystyle s\geq p}
p
≥
5
{\displaystyle p\geq 5}
s
>
p
{\displaystyle s>p}
p
1
2
3
4
5
6
7
8
min
s
1
2
3
4
6
7
9
11
{\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는 이 방법에 맞아들어가며, 그 테이블은 다음과 같다.
0
1
/
2
1
/
2
1
/
2
0
1
/
2
1
0
0
1
1
/
6
1
/
3
1
/
3
1
/
6
{\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 테이블은 다음과 같다.
0
1
/
3
1
/
3
2
/
3
−
1
/
3
1
1
1
−
1
1
1
/
8
3
/
8
3
/
8
1
/
8
{\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}}}
그러나, 가장 간단한 룽게-쿠타 방법은 (전방)오일러 방법 이다. 공식은
y
n
+
1
=
y
n
+
h
f
(
t
n
,
y
n
)
{\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n})}
0
1
{\displaystyle {\begin{array}{c|cc}0\\\hline &1\end{array}}}
간단한 예로 두 단계의 2차 방법인 중간점 방법 을 보면 다음과 같다.
y
n
+
1
=
y
n
+
h
f
(
t
n
+
1
2
h
,
y
n
+
1
2
h
f
(
t
n
,
y
n
)
)
.
{\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 )}.}
대응되는 테이블은 다음과 같다.
0
1
/
2
1
/
2
0
1
{\displaystyle {\begin{array}{c|cc}0\\1/2&1/2\\\hline &0&1\\\end{array}}}
두 단계의 2차 룽게-쿠타 방법은 중간점 방법 밖에 없는 것은 아니다. 이런 종류의 방법은 α에 의해 다음과 같이 매개변수화되어 있다.
y
n
+
1
=
y
n
+
h
(
(
1
−
1
2
α
)
f
(
t
n
,
y
n
)
+
1
2
α
f
(
t
n
+
α
h
,
y
n
+
α
h
f
(
t
n
,
y
n
)
)
)
{\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 테이블은 다음과 같다.
0
α
α
(
1
−
1
2
α
)
1
2
α
{\displaystyle {\begin{array}{c|cc}0\\\alpha &\alpha \\\hline &(1-{\frac {1}{2\alpha }})&{\frac {1}{2\alpha }}\\\end{array}}}
이 때
α
=
1
2
{\displaystyle \alpha ={\frac {1}{2}}}
α
=
1
{\displaystyle \alpha =1}
호인의 방법 이 된다.
예를 들어 α=2/3인 두 단계의 2차 룽게 쿠타 방법을 보자, 이것은 라슬톤 방법 이라고 알려져있다. 테이블은 다음과 같이 주어져 있다.
0
2
/
3
2
/
3
1
/
4
3
/
4
{\displaystyle {\begin{array}{c|cc}0\\2/3&2/3\\\hline &1/4&3/4\\\end{array}}}
대응되는 등식은 다음과 같다.
k
1
=
f
(
t
n
,
y
n
)
,
k
2
=
f
(
t
n
+
2
3
h
,
y
n
+
2
3
h
k
1
)
,
y
n
+
1
=
y
n
+
h
(
1
4
k
1
+
3
4
k
2
)
.
{\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}}}
다음의 초기치 문제를 풀기 위하여 이 방법을 사용한다.
y
′
=
tan
(
y
)
+
1
,
y
0
=
1
,
t
∈
[
1
,
1.1
]
{\displaystyle y^{\prime }=\tan(y)+1,\quad y_{0}=1,t\in [1,1.1]}
h = 0.025의 크기로 구간을 잡으면 이 방법은 네 단계가 필요하다.
이 방법의 결과는 다음과 같다:
t
0
=
1
:
y
0
=
1
t
1
=
1.025
:
y
0
=
1
k
1
=
2.557407725
k
2
=
f
(
t
0
+
2
3
h
,
y
0
+
2
3
h
k
1
)
=
2.713898184
y
1
=
y
0
+
h
(
1
4
k
1
+
3
4
k
2
)
=
1.066869388
_
t
2
=
1.05
:
y
1
=
1.066869388
k
1
=
2.813524695
k
2
=
f
(
t
1
+
2
3
h
,
y
1
+
2
3
h
k
1
)
y
2
=
y
1
+
h
(
1
4
k
1
+
3
4
k
2
)
=
1.141332181
_
t
3
=
1.075
:
y
2
=
1.141332181
k
1
=
3.183536647
k
2
=
f
(
t
2
+
2
3
h
,
y
2
+
2
3
h
k
1
)
y
3
=
y
2
+
h
(
1
4
k
1
+
3
4
k
2
)
=
1.227417567
_
t
4
=
1.1
:
y
3
=
1.227417567
k
1
=
3.796866512
k
2
=
f
(
t
3
+
2
3
h
,
y
3
+
2
3
h
k
1
)
y
4
=
y
3
+
h
(
1
4
k
1
+
3
4
k
2
)
=
1.335079087
_
.
{\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한 방법은 단일 룽게-쿠타 방법의 지역적인 절단 오차의 추정치를 찾기 위해 고안되었다. 이것은 테이블에서 차수가
p
{\displaystyle p}
p
−
1
{\displaystyle p-1}
낮은 차수의 단계는 다음과 같이 주어진다.
y
n
+
1
∗
=
y
n
+
h
∑
i
=
1
s
b
i
∗
k
i
,
{\displaystyle y_{n+1}^{*}=y_{n}+h\sum _{i=1}^{s}b_{i}^{*}k_{i},}
k
i
{\displaystyle k_{i}}
e
n
+
1
=
y
n
+
1
−
y
n
+
1
∗
=
h
∑
i
=
1
s
(
b
1
−
b
i
∗
)
k
i
,
{\displaystyle e_{n+1}=y_{n+1}-y_{n+1}^{*}=h\sum _{i=1}^{s}(b_{1}-b_{i}^{*})k_{i},}
오차는
O
(
h
p
)
{\displaystyle O(h^{p})}
b
i
8
{\displaystyle b_{i}^{8}}
0
c
2
a
21
c
3
a
31
a
32
⋮
⋮
⋱
c
s
a
s
1
a
s
2
⋯
a
s
,
s
−
1
b
1
b
2
⋯
b
s
−
1
b
s
b
1
∗
b
2
∗
⋯
b
s
−
1
∗
b
s
∗
{\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인 두 방법을 가지고 있다. 그 확장된 부처 태블로 는 다음과 같다:
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
{\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 테이블은 다음과 같다:
0
1
1
1
/
2
1
/
2
1
0
{\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 룽게-쿠타 방법은 다음의 형태를 가진다.
y
n
+
1
=
y
n
+
h
∑
i
=
1
s
b
i
k
i
{\displaystyle 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
∑
j
=
1
s
a
i
j
k
j
)
,
i
=
1
,
…
,
s
.
{\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 가 올라가고 계수들의 행렬은 다음과 같이 삼각행렬이 아니다.
c
1
a
11
a
12
⋯
a
1
s
c
2
a
21
a
22
⋯
a
2
s
⋮
⋮
⋮
⋱
⋮
c
s
a
s
1
a
s
2
⋯
a
s
s
b
1
b
2
⋯
b
s
b
1
∗
b
2
∗
⋯
b
s
∗
=
c
A
b
T
{\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 룽게-쿠타 방법은 역 오일러 방법 이다:
y
n
+
1
=
y
n
+
h
f
(
t
n
+
h
,
y
n
+
1
)
.
{\displaystyle y_{n+1}=y_{n}+hf(t_{n}+h,y_{n+1}).}
Butcher 테이블은 다음과 같다.
1
1
1
{\displaystyle {\begin{array}{c|c}1&1\\\hline &1\\\end{array}}}
이 Butcher 테이블에 대응하는 식은 다음과 같다.
k
1
=
f
(
t
n
+
h
,
y
n
+
h
k
1
)
a
n
d
y
n
+
1
=
y
n
+
h
k
1
{\displaystyle k_{1}=f(t_{n}+h,y_{n}+hk_{1})\quad and\quad y_{n+1}=y_{n}+hk_{1}}
이것은 위에서 나온 역 오일러 방법의 형태로 바꿀 수 있다.
다른 implicit 룽게-쿠타 방법의 예는 사다리꼴 법칙 이다. 그 Butcher 테이블은 다음과 같다.
0
0
0
1
1
2
1
2
1
2
1
2
1
0
{\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 테이블은 다음과 같다:
1
2
−
1
6
3
1
4
1
4
−
1
6
3
1
2
+
1
6
3
1
4
+
1
6
3
1
4
1
2
1
2
1
2
+
1
2
3
1
2
−
1
2
3
{\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 를 고려해 보자. 룽게-쿠타 방법을 가용하면 이 식은 다음
y
n
+
1
=
r
(
h
λ
)
y
n
{\displaystyle y_{n+1}=r(h\lambda )y_{n}}
r 은 다음과 같다.
r
(
z
)
=
1
+
z
b
T
(
I
−
z
A
)
−
1
e
=
det
(
I
−
z
A
+
z
e
b
T
)
det
(
I
−
z
A
)
,
{\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이라는 것을 의미하며 그 안정도 함수는 다항함수이다.
위의 선형방정식은
|
r
(
z
)
|
<
1
{\displaystyle |r(z)|<1}
z
=
h
λ
{\displaystyle z=h\lambda }
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 ( " \n t \t y" );
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 \t y \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 ()
![{\displaystyle y^{\prime }=\tan(y)+1,\quad y_{0}=1,t\in [1,1.1]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/dafd44e47f94754f8638c04ed5a91b2cb5485aff)
