7.50 강성 방정식과 암시적 해법

1. 강성의 정의와 특성

1.1 강성 시스템의 개념

강성(stiffness)은 상미분 방정식 시스템에서 서로 극단적으로 다른 시간 규모(time scale)가 공존하는 현상을 지칭한다. 선형 시스템 \dot{\mathbf{y}} = A\mathbf{y}에서 시스템 행렬 A의 고유값이 \lambda_1, \lambda_2, \ldots, \lambda_n일 때, 강성비(stiffness ratio)는 다음과 같이 정의된다.

\rho = \frac{\max_i \lvert\text{Re}(\lambda_i)\rvert}{\min_i \lvert\text{Re}(\lambda_i)\rvert}

강성비가 매우 큰 경우(\rho \gg 1), 시스템은 빠르게 감쇠하는 모드와 느리게 변화하는 모드를 동시에 포함한다. 명시적 방법의 안정성 조건은 가장 빠른 모드에 의해 결정되므로, 느린 모드를 추적하기 위한 긴 시뮬레이션 구간에서도 매우 작은 시간 간격을 사용하여야 한다. 이는 계산적으로 극히 비효율적이다.

1.2 강성의 정성적 기준

Hairer와 Wanner(1996)에 따르면, 강성 문제는 다음의 조건을 동시에 만족하는 경우이다.

1. 명시적 방법의 안정성 조건이 정확도가 아닌 안정성에 의해 시간 간격을 제한한다.
2. 암시적 방법이 명시적 방법보다 훨씬 효율적으로 문제를 풀 수 있다.

즉, 해 자체는 완만하게 변화하지만, 명시적 방법이 요구하는 시간 간격이 해의 변화 규모에 비하여 불합리하게 작은 경우가 강성 문제의 본질이다.

2. 로봇 시스템에서의 강성 발생

2.1 유연 관절 로봇

유연 관절(flexible joint) 모델에서 관절 강성 K가 매우 큰 경우, 시스템의 고유 진동수 중 일부는 \sqrt{K/J}에 비례하여 매우 높아진다. 반면 링크의 강체 운동 모드는 상대적으로 낮은 주파수를 갖는다. 이로 인하여 관절 강성이 증가할수록 강성비가 증가하며, 명시적 방법의 시간 간격은 고주파 모드에 의해 심하게 제한된다.

2.2 접촉 역학

로봇이 환경과 접촉하는 경우, 접촉 강성(contact stiffness)이 매우 높은 접촉력 모델을 사용하면 강성 방정식이 발생한다. 패널티 기반(penalty-based) 접촉 모델에서 접촉 강성 계수 k_c가 크면, 접촉 모드의 고유 진동수가 매우 높아져 시스템이 극도로 강성이 된다.

2.3 구동기 동역학

DC 모터의 전기적 시정수(L/R, 여기서 L은 인덕턴스, R은 저항)는 기계적 시정수에 비하여 수십에서 수백 배 작다. 전기-기계 결합 시스템의 동역학을 함께 모델링하면 필연적으로 강성 시스템이 형성된다.

3. 명시적 방법의 한계

3.1 안정성 제한의 수학적 분석

고전적 RK4의 안정성 영역은 실수 축을 따라 대략 -2.785 \leq h\lambda \leq 0의 범위로 제한된다. 강성비가 \rho = 10^6인 시스템에서 최대 고유값 크기가 \lvert\lambda_{max}\rvert = 10^6이면, 시간 간격은 h < 2.785 \times 10^{-6}으로 제한된다. 그러나 느린 모드의 해를 정확히 추적하기 위해서는 h \approx 10^{-2} 정도면 충분하다. 이 경우 안정성이 요구하는 시간 간격은 정확도가 요구하는 시간 간격의 약 10^{-4}배에 불과하며, 이는 막대한 계산 낭비이다.

4. A-안정성

4.1 정의

수치 방법이 A-안정(A-stable)하다 함은, 달퀴스트 시험 방정식 \dot{y} = \lambda y에 대하여 \text{Re}(\lambda) \leq 0인 모든 \lambda와 모든 양의 시간 간격 h > 0에 대하여 수치 해가 유계(bounded)인 것을 의미한다. 즉, 안정성 영역이 복소 평면의 왼쪽 반평면(left half-plane) 전체를 포함하여야 한다.

\{z \in \mathbb{C} : \text{Re}(z) \leq 0\} \subseteq \{z \in \mathbb{C} : \lvert R(z) \rvert \leq 1\}

4.2 달퀴스트 제2 장벽

달퀴스트(Dahlquist, 1963)는 다음의 중요한 정리를 증명하였다.

달퀴스트 제2 장벽(Dahlquist’s second barrier). 명시적 선형 다단계 방법(explicit linear multistep method)은 A-안정할 수 없다. A-안정한 선형 다단계 방법 중 최고 차수는 2차이다.

이 정리는 강성 문제를 풀기 위해서는 암시적 방법이 필수적임을 이론적으로 확립한다.

5. 암시적 방법

5.1 후진 오일러 방법

가장 단순한 A-안정 방법은 후진 오일러 방법(backward Euler method)이다.

y_{n+1} = y_n + hf(t_{n+1}, y_{n+1})

안정성 함수는 R(z) = 1/(1 - z)이며, \text{Re}(z) \leq 0에서 \lvert R(z) \rvert \leq 1이 항상 성립한다. 또한 \lvert z \rvert \to \infty일 때 \lvert R(z) \rvert \to 0이므로, 매우 빠른 모드를 강하게 감쇠시킨다. 이 성질을 L-안정성(L-stability)이라 한다.

5.2 사다리꼴 규칙

암시적 사다리꼴 규칙(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]

안정성 함수는 R(z) = (1 + z/2)/(1 - z/2)이며, 이는 A-안정하다. 그러나 \lvert z \rvert \to \infty일 때 \lvert R(z) \rvert \to 1이므로 L-안정하지 않다. 이로 인하여 매우 큰 고유값에 대응하는 모드가 충분히 감쇠되지 않아 수치적 진동이 발생할 수 있다.

5.3 암시적 룽게-쿠타 방법

암시적 룽게-쿠타(implicit Runge-Kutta, IRK) 방법은 부처 배열의 행렬 A가 일반적인(상삼각 원소도 포함하는) 형태를 갖는다. 대표적인 예로 가우스-르장드르(Gauss-Legendre) 방법이 있으며, s단계에서 2s차 정확도를 달성한다.

2단계 가우스-르장드르 방법의 부처 배열은 다음과 같다.

\begin{array}{c|cc} \frac{1}{2} - \frac{\sqrt{3}}{6} & \frac{1}{4} & \frac{1}{4} - \frac{\sqrt{3}}{6} \\[4pt] \frac{1}{2} + \frac{\sqrt{3}}{6} & \frac{1}{4} + \frac{\sqrt{3}}{6} & \frac{1}{4} \\[4pt] \hline & \frac{1}{2} & \frac{1}{2} \end{array}

이 방법은 4차 정확도를 가지면서 A-안정하고 심플렉틱(symplectic)이다.

6. 비선형 방정식의 풀이

6.1 뉴턴-랩슨 반복법

암시적 방법에서는 각 시간 단계에서 비선형 방정식을 풀어야 한다. 후진 오일러 방법의 경우, 다음 방정식의 해 y_{n+1}을 구하여야 한다.

G(y_{n+1}) \equiv y_{n+1} - y_n - hf(t_{n+1}, y_{n+1}) = 0

뉴턴-랩슨(Newton-Raphson) 반복법을 적용하면 다음과 같다.

y_{n+1}^{(k+1)} = y_{n+1}^{(k)} - \left[G'(y_{n+1}^{(k)})\right]^{-1} G(y_{n+1}^{(k)})

여기서 야코비안 행렬은 다음과 같다.

G' = I - h\frac{\partial f}{\partial y}

초기 추정값으로는 명시적 오일러 예측값 y_{n+1}^{(0)} = y_n + hf(t_n, y_n)을 사용하는 것이 일반적이다.

6.2 연립 방정식에서의 야코비안

n차원 연립 방정식 \dot{\mathbf{y}} = \mathbf{f}(t, \mathbf{y})에 대한 후진 오일러 방법에서, 뉴턴 반복의 각 단계는 다음의 n \times n 선형 시스템을 풀어야 한다.

(I - hJ)\Delta\mathbf{y}^{(k)} = -\mathbf{G}(\mathbf{y}_{n+1}^{(k)})

여기서 J = \partial\mathbf{f}/\partial\mathbf{y}는 야코비안 행렬이다. 이 선형 시스템의 풀이가 암시적 방법의 주된 계산 비용을 구성한다.

6.3 야코비안 갱신 전략

야코비안 행렬의 계산과 LU 분해는 비용이 크므로, 매 뉴턴 반복마다 야코비안을 갱신하지 않고 이전에 계산한 야코비안을 재사용하는 전략이 사용된다. 이를 수정 뉴턴법(modified Newton method)이라 한다. 야코비안은 수렴이 느려지거나, 시간 간격이 변경되거나, 일정 단계 수가 경과한 경우에만 갱신한다.

7. SDIRK 방법

7.1 정의와 구조

단일 대각 암시적 룽게-쿠타(singly diagonally implicit Runge-Kutta, SDIRK) 방법은 부처 배열의 행렬 A가 하삼각이고 대각 원소가 모두 동일한 값 \gamma를 갖는 방법이다.

A = \begin{bmatrix} \gamma & 0 & \cdots & 0 \\ a_{21} & \gamma & \cdots & 0 \\ \vdots & \vdots & \ddots & 0 \\ a_{s1} & a_{s2} & \cdots & \gamma \end{bmatrix}

7.2 계산 효율

SDIRK 방법의 핵심 장점은 모든 단계에서 동일한 행렬 (I - h\gamma J)의 분해를 재사용할 수 있다는 점이다. 완전 암시적 룽게-쿠타 방법에서는 sn \times sn 크기의 선형 시스템을 풀어야 하는 반면, SDIRK에서는 n \times n 크기의 시스템을 s회 순차적으로 풀면 된다. 이는 계산 비용을 대폭 절감한다.

7.3 L-안정한 SDIRK 방법

\gamma = (2 - \sqrt{2})/2 \approx 0.2929를 선택하면 2단계 3차 SDIRK 방법이 L-안정하게 된다. 이 방법은 강성 로봇 동역학 시뮬레이션에서 널리 사용되는 적분기 중 하나이다.

8. BDF 방법

8.1 후진 차분 공식

후진 차분 공식(backward differentiation formula, BDF)은 이전 여러 시간 단계의 해를 이용하여 도함수를 근사하는 다단계 암시적 방법이다. k단계 BDF의 일반 형태는 다음과 같다.

\sum_{j=0}^{k} \alpha_j y_{n+1-j} = h\beta_0 f(t_{n+1}, y_{n+1})

여기서 \beta_0 \neq 0이고 \beta_j = 0 (j \geq 1)이다. 즉, 우변에 f의 현재 시각 값만 포함되므로, 뉴턴 반복에서 n \times n 선형 시스템만을 풀면 된다.

8.2 차수별 BDF 공식

1차 및 2차 BDF는 A-안정하며, 3~6차 BDF는 A(\alpha)-안정하다. 7차 이상의 BDF는 영안정(zero-stable)하지 않으므로 사용할 수 없다.

BDF 방법은 MATLAB의 ode15s, Sundials의 CVODE 등 주요 강성 적분기에 채택되어 있다.

9. 요약

강성 방정식은 서로 다른 시간 규모의 동적 모드가 공존하는 시스템에서 발생하며, 유연 관절, 접촉 역학, 구동기 전기 동역학 등 로봇공학의 다양한 영역에서 나타난다. 명시적 방법은 안정성 제한에 의해 비효율적이므로, 강성 문제에서는 A-안정한 암시적 방법이 필수적이다. 후진 오일러, 사다리꼴 규칙, SDIRK, BDF 등의 암시적 방법은 각 시간 단계에서 비선형 방정식 풀이를 필요로 하며, 뉴턴 반복과 야코비안 분해가 주된 계산 비용이다. SDIRK 방법은 단일 행렬 분해의 재사용으로 계산 효율을 높이며, BDF 방법은 다단계 구조를 통해 높은 차수의 정확도를 제공한다.

참고 문헌

• Hairer, E. & Wanner, G. (1996). Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems (2nd ed.). Springer.
• Ascher, U. M. & Petzold, L. R. (1998). Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations. SIAM.
• Butcher, J. C. (2016). Numerical Methods for Ordinary Differential Equations (3rd ed.). Wiley.
• Dahlquist, G. (1963). “A special stability problem for linear multistep methods.” BIT Numerical Mathematics, 3, 27–43.
• Shampine, L. F. & Reichelt, M. W. (1997). “The MATLAB ODE Suite.” SIAM Journal on Scientific Computing, 18(1), 1–22.
• Featherstone, R. (2008). Rigid Body Dynamics Algorithms. Springer.

version: 0.1.0