7.47 개선 오일러 방법과 2차 룽게-쿠타법

1. 개선 오일러 방법의 동기

전진 오일러 방법은 시간 구간 [t_n, t_{n+1}]의 시작점에서의 기울기만을 사용하여 해를 외삽한다. 이로 인하여 구간 내에서 기울기가 변화하는 경우 상당한 오차가 발생한다. 개선 오일러 방법(improved Euler method)은 구간의 시작점과 끝점에서의 기울기 정보를 함께 활용하여 정확도를 향상시킨다.

2. 호인 방법

2.1 유도와 공식

호인 방법(Heun’s method)은 개선 오일러 방법의 대표적 형태이다. 초기값 문제 \dot{y} = f(t, y)에 대하여 다음 두 단계로 구성된다.

예측 단계(predictor): 전진 오일러 방법으로 다음 시각의 예비값을 계산한다.

\tilde{y}_{n+1} = y_n + hf(t_n, y_n)

수정 단계(corrector): 시작점과 예측된 끝점에서의 기울기 평균을 사용하여 수정한다.

y_{n+1} = y_n + \frac{h}{2}\left[f(t_n, y_n) + f(t_{n+1}, \tilde{y}_{n+1})\right]

이 방법은 사다리꼴 공식(trapezoidal rule)의 명시적 근사로 해석할 수 있다. 정적분 \int_{t_n}^{t_{n+1}} f(t, y(t))\,dt를 사다리꼴 규칙으로 근사하되, 끝점에서의 미지 함수값을 오일러 방법으로 예측하는 것이다.

2.2 국소 절단 오차

호인 방법의 국소 절단 오차를 유도하기 위하여 y(t_{n+1})의 테일러 전개를 2차항까지 수행한다.

y(t_{n+1}) = y(t_n) + hy'(t_n) + \frac{h^2}{2}y''(t_n) + O(h^3)

여기서 y' = f이고, y'' = f_t + ff_y (f_t = \partial f/\partial t, f_y = \partial f/\partial y)이다.

호인 방법의 갱신 공식을 테일러 전개하면 다음과 같다.

y_{n+1} = y_n + hf + \frac{h^2}{2}(f_t + ff_y) + O(h^3)

이는 정확한 테일러 전개와 2차항까지 일치하므로, 국소 절단 오차는 O(h^3)이고 전역 절단 오차는 O(h^2)이다. 따라서 호인 방법은 2차 방법이다.

3. 수정 오일러 방법

3.1 중점 방법

수정 오일러 방법(modified Euler method) 또는 중점 방법(midpoint method)은 구간의 중점에서의 기울기를 사용하는 방법이다.

1단계: 구간 중점에서의 값을 오일러 방법으로 예측한다.

y_{n+1/2} = y_n + \frac{h}{2}f(t_n, y_n)

2단계: 중점에서의 기울기로 전체 구간을 진행한다.

y_{n+1} = y_n + hf\left(t_n + \frac{h}{2}, \, y_{n+1/2}\right)

이 방법도 2차 정확도를 가지며, 국소 절단 오차는 O(h^3)이다.

4. 차 룽게-쿠타법의 일반 형태

4.1 일반 공식

2차 룽게-쿠타(Runge-Kutta) 방법의 일반적인 2단계(two-stage) 형태는 다음과 같다.

k_1 = f(t_n, y_n)

k_2 = f(t_n + \alpha h, \, y_n + \alpha h k_1)

y_{n+1} = y_n + h(w_1 k_1 + w_2 k_2)

여기서 \alpha, w_1, w_2는 방법을 정의하는 매개변수이다.

4.2 차수 조건의 유도

이 일반 공식이 2차 정확도를 갖기 위하여, 테일러 전개를 통해 다음의 조건이 유도된다.

w_1 + w_2 = 1

w_2 \alpha = \frac{1}{2}

이 조건은 2개의 방정식에 3개의 미지수(w_1, w_2, \alpha)를 포함하므로, 1개의 자유 매개변수가 존재한다. 즉, 2차 룽게-쿠타 방법은 무한히 많은 변형이 가능하다.

4.3 주요 변형

자유 매개변수의 선택에 따른 대표적 변형은 다음과 같다.

랄스턴(Ralston) 방법은 국소 절단 오차의 상계를 최소화하는 매개변수 선택에 해당한다.

5. 부처 배열

5.1 표기법

룽게-쿠타 방법의 계수를 체계적으로 표현하기 위하여 부처 배열(Butcher tableau)을 사용한다. s단계 룽게-쿠타 방법의 부처 배열은 다음과 같은 형태이다.

\begin{array}{c|c} \mathbf{c} & A \\ \hline & \mathbf{b}^T \end{array}

여기서 A = [a_{ij}] \in \mathbb{R}^{s \times s}는 룽게-쿠타 행렬, \mathbf{b} = [b_1, \ldots, b_s]^T는 가중치 벡터, \mathbf{c} = [c_1, \ldots, c_s]^T는 노드 벡터이다.

5.2 호인 방법의 부처 배열

호인 방법의 부처 배열은 다음과 같다.

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

5.3 중점 방법의 부처 배열

중점 방법의 부처 배열은 다음과 같다.

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

명시적 방법에서는 행렬 A가 강하게 하삼각(strictly lower triangular)이며, 이는 k_i의 계산이 이전에 구한 k_1, \ldots, k_{i-1}에만 의존함을 의미한다.

6. 안정성 해석

6.1 시험 방정식에 대한 안정성 함수

달퀴스트 시험 방정식 \dot{y} = \lambda y에 2차 룽게-쿠타 방법을 적용하면, 호인 방법의 경우 안정성 함수는 다음과 같다.

R(z) = 1 + z + \frac{z^2}{2}

여기서 z = h\lambda이다. 이는 e^z의 2차 테일러 다항식에 해당한다. 안정성 영역 \{z \in \mathbb{C} : \lvert R(z) \rvert \leq 1\}은 전진 오일러 방법의 안정성 영역보다 넓다.

6.2 실수 축에서의 안정성 구간

실수 고유값 \lambda < 0에 대하여 \lvert R(z) \rvert \leq 1을 만족하는 z = h\lambda의 범위는 -2 \leq z \leq 0이다. 따라서 안정성 조건은 다음과 같다.

0 < h < \frac{2}{\lvert\lambda\rvert}

이는 전진 오일러 방법과 동일한 실수 축 안정성 구간이다. 그러나 허수 축 방향으로는 2차 방법의 안정성 영역이 더 넓으므로, 진동성 시스템에서 상대적으로 유리하다.

7. 계산 비용과 정확도의 비교

2차 룽게-쿠타 방법은 각 시간 단계에서 함수 f를 2회 평가한다. 전진 오일러 방법은 1회 평가로 1차 정확도를 달성하고, 2차 룽게-쿠타 방법은 2회 평가로 2차 정확도를 달성한다.

동일한 전역 오차 수준 \epsilon을 달성하기 위해 필요한 총 함수 평가 횟수를 비교하면 다음과 같다.

• 오일러 방법: N = O(1/\epsilon)이므로 총 함수 평가 O(1/\epsilon)회
• 2차 룽게-쿠타: N = O(1/\sqrt{\epsilon})이므로 총 함수 평가 O(1/\sqrt{\epsilon}) \times 2 = O(1/\sqrt{\epsilon})회

따라서 높은 정밀도가 요구될수록(\epsilon이 작을수록) 2차 방법의 계산 효율이 오일러 방법에 비하여 현저히 우수하다.

8. 연립 방정식에의 적용

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

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

\mathbf{k}_2 = \mathbf{f}(t_n + h, \, \mathbf{y}_n + h\mathbf{k}_1)

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

로봇 동역학의 상태 벡터 \mathbf{x} = [\mathbf{q}^T, \dot{\mathbf{q}}^T]^T에 대하여 각 \mathbf{k}_i를 계산할 때, 상태 벡터의 모든 성분이 동시에 갱신되어야 한다. 위치 성분과 속도 성분을 분리하여 순차적으로 갱신하면 일관성이 깨질 수 있으므로 주의하여야 한다.

9. 예측자-수정자 해석

호인 방법은 예측자-수정자(predictor-corrector) 쌍으로 해석할 수 있다. 예측자는 전진 오일러 방법이고, 수정자는 암시적 사다리꼴 규칙의 1회 반복 근사이다.

암시적 사다리꼴 규칙(implicit trapezoidal rule)은 다음과 같다.

y_{n+1} = y_n + \frac{h}{2}\left[f(t_n, y_n) + f(t_{n+1}, y_{n+1})\right]

이 공식은 우변에 y_{n+1}이 포함되어 있으므로 암시적이다. 호인 방법은 이 암시적 방정식의 우변에서 y_{n+1}을 오일러 예측값 \tilde{y}_{n+1}로 대체한 명시적 근사이다. 수정 단계를 반복적으로 적용하면 암시적 사다리꼴 규칙의 해에 수렴하며, 이를 반복 예측자-수정자(iterated predictor-corrector) 방법이라 한다.

10. 요약

개선 오일러 방법과 2차 룽게-쿠타법은 구간 내 복수 지점의 기울기 정보를 활용하여 전진 오일러 방법의 정확도를 1차에서 2차로 향상시킨다. 호인 방법, 중점 방법, 랄스턴 방법은 모두 동일한 2차 차수 조건을 만족하는 변형이며, 부처 배열로 체계적으로 분류된다. 함수 평가 횟수가 2배로 증가하나, 시간 간격을 더 크게 설정할 수 있으므로 동일 정밀도에서의 총 계산 비용은 오히려 감소한다. 이러한 계산 효율의 향상은 로봇 동역학 시뮬레이션의 실시간 구현에서 중요한 의미를 갖는다.

참고 문헌

• 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.
• Süli, E. & Mayers, D. F. (2003). An Introduction to Numerical Analysis. Cambridge University Press.
• Ascher, U. M. & Petzold, L. R. (1998). Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations. SIAM.
• Ralston, A. (1962). “Runge-Kutta methods with minimum error bounds.” Mathematics of Computation, 16(80), 431–437.

version: 0.1.0