7.48 4차 룽게-쿠타법(RK4)

1. 고전적 4차 룽게-쿠타법의 공식

4차 룽게-쿠타법(fourth-order Runge-Kutta method, RK4)은 수치 적분에서 가장 널리 사용되는 방법이다. 초기값 문제 $\dot{y} = f(t, y)$ , $y(t_0) = y_0$ 에 대하여, 고전적 RK4의 갱신 공식은 다음과 같다.

$y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$

여기서 네 개의 단계 기울기(stage slope)는 순차적으로 계산된다.

$k_1 = f(t_n, y_n)$

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

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

$k_4 = f(t_n + h, \, y_n + hk_3)$

최종 갱신은 이 네 기울기의 가중 평균이며, 가중치 $1/6, 2/6, 2/6, 1/6$ 은 심프슨 적분 공식(Simpson’s rule)의 가중치와 일치한다.

2. 부처 배열

고전적 RK4의 부처 배열(Butcher tableau)은 다음과 같다.

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

이 배열에서 룽게-쿠타 행렬 $A$ 가 강하게 하삼각(strictly lower triangular)이므로 RK4는 명시적(explicit) 방법이다. 각 단계 $k_i$ 는 이전 단계의 결과만을 참조하여 순차적으로 계산할 수 있다.

3. 절단 오차 해석

3.1 차수 조건

$s$ 단계 룽게-쿠타 방법이 $p$ 차 정확도를 갖기 위해서는 테일러 급수 전개에서 $p$ 차항까지의 일치를 보장하는 차수 조건(order conditions)을 만족하여야 한다. $p$ 차 방법의 차수 조건 수는 $p$ 가 증가함에 따라 급격히 증가한다.

차수 $p$	조건의 수	최소 단계 수 $s$
1	1	1
2	2	2
3	4	3
4	8	4
5	17	6

4차 방법은 8개의 차수 조건을 4단계로 만족시킬 수 있으므로, 단계 수와 차수가 일치하는 마지막 경우이다. 5차 이상에서는 최소 단계 수가 차수보다 크며, 이를 부처 장벽(Butcher barrier)이라 한다.

3.2 국소 절단 오차와 전역 절단 오차

RK4의 국소 절단 오차는 $O(h^5)$ 이고, 전역 절단 오차는 $O(h^4)$ 이다. 시간 간격을 절반으로 줄이면 전역 오차는 대략 $1/16$ 로 감소한다.

$\frac{e(h)}{e(h/2)} \approx 2^4 = 16$

이 높은 차수의 정확도가 RK4가 실용적으로 널리 채택되는 핵심 이유이다.

3.3 오차 상수

고전적 RK4의 국소 절단 오차의 주 항(leading term)은 다음과 같이 표현된다.

$\tau_{n+1} = \sum_{i} c_i h^5 D_i f + O(h^6)$

여기서 $D_i f$ 는 $f$ 의 다양한 편도함수의 조합이고, $c_i$ 는 방법에 의존하는 상수이다. 오차 상수(error constant)는 이 주 항의 계수의 2-노름으로 정의되며, 고전적 RK4의 오차 상수는 약 $1/120$ 이다.

4. 안정성 해석

4.1 안정성 함수

달퀴스트 시험 방정식 $\dot{y} = \lambda y$ 에 RK4를 적용하면, 안정성 함수는 $e^z$ 의 4차 테일러 다항식이 된다.

$R(z) = 1 + z + \frac{z^2}{2} + \frac{z^3}{6} + \frac{z^4}{24}$

여기서 $z = h\lambda$ 이다.

4.2 안정성 영역

안정성 영역 $\{z \in \mathbb{C} : \lvert R(z) \rvert \leq 1\}$ 은 실수 축을 따라 대략 $-2.785 \leq \text{Re}(z) \leq 0$ 의 범위를 갖는다. 허수 축 방향으로는 대략 $\lvert\text{Im}(z)\rvert \leq 2.828$ 까지 안정하다.

실수 고유값 $\lambda < 0$ 에 대한 안정성 조건은 다음과 같다.

$0 < h < \frac{2.785}{\lvert\lambda\rvert}$

2차 방법의 안정성 구간( $h < 2/\lvert\lambda\rvert$ )에 비하여 약 39% 넓은 안정성 범위를 갖는다.

4.3 순수 허수 고유값에 대한 안정성

순수 허수 고유값 $\lambda = j\omega$ 에 대하여 $z = jh\omega$ 를 대입하면 다음과 같다.

$\lvert R(jh\omega)\rvert^2 = 1 + \frac{(h\omega)^8}{576}$

이 값은 항상 1보다 크므로, RK4는 순수 허수 고유값에 대하여 불안정하다. 이는 보존계(conservative system)의 장시간 시뮬레이션에서 에너지가 점진적으로 증가하는 수치적 불안정성을 초래할 수 있다.

5. 연립 방정식에의 적용

벡터 형태의 연립 방정식 $\dot{\mathbf{y}} = \mathbf{f}(t, \mathbf{y})$ 에 대하여 RK4는 다음과 같이 적용된다.

$\mathbf{k}_1 = \mathbf{f}(t_n, \mathbf{y}_n)$

$\mathbf{k}_2 = \mathbf{f}\left(t_n + \frac{h}{2}, \, \mathbf{y}_n + \frac{h}{2}\mathbf{k}_1\right)$

$\mathbf{k}_3 = \mathbf{f}\left(t_n + \frac{h}{2}, \, \mathbf{y}_n + \frac{h}{2}\mathbf{k}_2\right)$

$\mathbf{k}_4 = \mathbf{f}(t_n + h, \, \mathbf{y}_n + h\mathbf{k}_3)$

$\mathbf{y}_{n+1} = \mathbf{y}_n + \frac{h}{6}(\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4)$

각 단계에서 벡터 $\mathbf{f}$ 의 모든 성분이 동시에 평가되므로, 변수 간 결합(coupling)이 자연스럽게 처리된다.

6. 계산 비용 분석

RK4는 각 시간 단계에서 함수 $\mathbf{f}$ 를 4회 평가한다. $n$ 차원 상태 벡터에 대한 함수 평가 비용을 $C_f$ 라 하면, 한 단계의 총 비용은 $4C_f$ 이다.

동일 정밀도를 달성하기 위한 총 비용을 비교하면 다음과 같다. 전역 오차를 $\epsilon$ 이하로 유지하기 위하여 RK4는 $N = O(\epsilon^{-1/4})$ 단계가 필요하므로, 총 함수 평가 횟수는 $4N = O(\epsilon^{-1/4})$ 이다. 이에 비하여 오일러 방법은 $O(\epsilon^{-1})$ 회의 평가가 필요하다. 따라서 높은 정밀도가 요구될수록 RK4의 계산 효율이 오일러 방법에 비하여 월등히 우수하다.

7. 로봇 동역학 시뮬레이션에서의 적용

7.1 다자유도 로봇의 시뮬레이션

$n$ 자유도 로봇의 상태 벡터 $\mathbf{x} = [\mathbf{q}^T, \dot{\mathbf{q}}^T]^T \in \mathbb{R}^{2n}$ 에 대한 상태 방정식

$\dot{\mathbf{x}} = \begin{bmatrix} \dot{\mathbf{q}} \\ M(\mathbf{q})^{-1}[\boldsymbol{\tau} - C(\mathbf{q}, \dot{\mathbf{q}})\dot{\mathbf{q}} - \mathbf{g}(\mathbf{q})] \end{bmatrix}$

에 RK4를 적용할 때, 각 단계 $\mathbf{k}_i$ 의 계산에서 관성 행렬 $M(\mathbf{q})$ 의 역행렬 연산 또는 이와 동등한 선형 시스템 풀이가 필요하다. $\mathbf{k}_2$ 와 $\mathbf{k}_3$ 계산에서는 중간 시각의 관절 위치에 대한 관성 행렬을 새로이 구성하여야 하므로, 함수 평가 비용이 상당하다.

7.2 시간 간격 선택

RK4의 높은 정확도 덕분에 오일러 방법에 비하여 훨씬 큰 시간 간격을 사용할 수 있다. 일반적인 로봇 동역학 시뮬레이션에서 RK4는 1~10 ms의 시간 간격으로 충분한 정밀도를 제공한다. 다만, 시간 간격의 상한은 안정성 조건에 의해 제한되며, 시스템의 최대 고유 진동수에 의해 결정된다.

7.3 에너지 보존의 수치적 검증

보존계에서의 장시간 시뮬레이션 시, 시스템의 총 에너지를 매 시간 단계에서 계산하여 수치적 에너지 편이(drift)를 모니터링하는 것이 바람직하다. RK4는 비심플렉틱(non-symplectic) 방법이므로 에너지가 엄밀하게 보존되지 않으나, 짧은 시간 구간에서는 에너지 오차가 $O(h^4)$ 로 매우 작다.

8. 다른 4차 방법과의 비교

고전적 RK4 외에도 다양한 4단계 4차 룽게-쿠타 방법이 존재한다. 대표적인 변형으로 3/8 규칙(3/8-rule) RK4가 있으며, 그 부처 배열은 다음과 같다.

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

이 방법은 고전적 RK4와 동일한 차수와 유사한 오차 상수를 가지나, 고전적 RK4가 약간 더 작은 오차 상수를 갖는다. 따라서 실용적으로는 고전적 RK4가 더 널리 사용된다.

9. RK4의 한계와 적응적 방법의 필요성

RK4는 고정 시간 간격 방법으로, 해의 변화가 급격한 구간과 완만한 구간에서 동일한 시간 간격을 사용한다. 이는 급격한 변화 구간에서는 정밀도가 부족하고, 완만한 구간에서는 불필요하게 많은 계산을 수행하는 비효율을 초래한다.

또한 RK4는 국소 절단 오차의 추정치를 자체적으로 제공하지 않으므로, 시간 간격의 적절성을 판단할 수 있는 내재적 오차 지표가 없다. 이러한 한계를 극복하기 위하여 내장형(embedded) 룽게-쿠타 쌍을 이용한 적응적 시간 간격 제어 방법이 개발되었다.

10. 요약

4차 룽게-쿠타법(RK4)은 4단계 함수 평가로 4차 정확도를 달성하는 효율적인 수치 적분 방법이다. 국소 절단 오차 $O(h^5)$ 와 전역 절단 오차 $O(h^4)$ 를 가지며, 안정성 영역은 하위 차수 방법에 비하여 넓다. 부처 배열에 의한 체계적 표현과 벡터 형태로의 자연스러운 확장이 가능하여, 다자유도 로봇 동역학 시뮬레이션에 널리 사용된다. 다만 고정 시간 간격 사용의 비효율과 자체 오차 추정 부재라는 한계가 있으며, 이는 적응적 방법으로 보완할 수 있다.

참고 문헌

Butcher, J. C. (2016). Numerical Methods for Ordinary Differential Equations (3rd ed.). Wiley.
Hairer, E., Nørsett, S. P., & Wanner, G. (1993). Solving Ordinary Differential Equations I: Nonstiff Problems (2nd ed.). Springer.
Kutta, W. (1901). “Beitrag zur näherungsweisen Integration totaler Differentialgleichungen.” Zeitschrift für Mathematik und Physik, 46, 435–453.
Dormand, J. R. (1996). Numerical Methods for Differential Equations: A Computational Approach. CRC Press.
Featherstone, R. (2008). Rigid Body Dynamics Algorithms. Springer.

version: 0.1.0