미분대수방정식

1. 개요[편집]

미분대수방정식 Differential-Algebraic Equation (DAE)
일반형	$F(t, y, y') = 0$,   $\partial F/\partial y'$ 특이
난이도 지표	미분지수 (differential index) $\nu$
대표 사례	다물체 구속 동역학($\nu=3$), 비압축성 유동($\nu=2$), 회로 MNA($\nu \le 2$)
해법	BDF(DASSL·IDA), Radau IIA, 지수 축소 + 안정화
함정	드리프트, 일관되지 않은 초기값, 차수 저하

ODE는 “다음에 어디로 갈지”를 알려준다. DAE는 거기에 더해 “애초에 여기 있으면 안 된다”고 말한다.

미분대수방정식(Differential-Algebraic Equation, DAE)은 미분방정식과 대수 구속조건이 한 계 안에 섞여 있어, 일반형 $F(t, y, y') = 0$ 에서 야코비안 $\partial F / \partial y'$ 가 특이한 방정식계를 말한다. 즉 $y'$ 에 대해 풀 수 없는 성분이 존재한다. 가장 흔한 반명시적(semi-explicit) 형태는

$$x' = f(t, x, z), \qquad 0 = g(t, x, z)$$

로, $x$ 는 미분 변수, $z$ 는 대수 변수다. 두 번째 줄이 있다는 것만으로 룽게-쿠타법 교과서의 모든 정리가 통째로 무효화된다. 초기값을 아무렇게나 줄 수 없고, 국소 오차 추정이 안 먹고, 차수가 저하되고, 심지어 해가 존재하지 않을 수도 있다.

DAE가 튀어나오는 이유는 단순하다. 물리 모형을 만들 때 구속조건을 없애는 것보다 남겨두는 게 훨씬 쉽기 때문이다. 로봇 팔의 관절, 회로의 키르히호프 법칙, 비압축성 조건 — 전부 “이 상태공간 안의 어떤 부분다양체 위에만 존재하라”는 대수식이고, 이걸 최소 좌표계로 치환해 없애는 것은 대개 불가능하거나 지저분하다.

2. 미분지수[편집]

DAE의 난이도를 재는 표준 척도는 미분지수(differential index) $\nu$ 다. 정의는 “이 계를 시간에 대해 몇 번 미분해야 순수 ODE로 바뀌는가”의 최소 횟수. 지수가 1 올라갈 때마다 수치적 고통은 대략 한 등급씩 뛴다.

• 지수 0 — 그냥 ODE.
• 지수 1 — 반명시적 형태에서 $\partial g / \partial z$ 가 정칙인 경우. $g$ 를 한 번 미분하면 $z'$ 에 대해 풀린다. 실무적으로는 “거의 ODE”라서 뉴턴-랩슨법이 붙은 암시적 솔버면 대개 잘 돌아간다.
• 지수 2 — $\partial g/\partial z$ 가 아예 0인 헤센베르크 형 $x' = f(x,z),\ 0 = g(x)$ 에서 $g_x f_z$ 가 정칙인 경우. 비압축성 나비에-스토크스 방정식을 반이산화한 계가 여기 속하며, 압력이 비발산 조건의 라그랑주 승수 역할을 한다.¹
• 지수 3 — 구속된 역학계의 표준 위치. 아래 진자가 교과서 예제다.

진자를 데카르트 좌표로 쓰면 지수 3이 된다. 질량 $m$, 길이 $L$ 인 진자를

$$x' = u,\quad y' = v,\quad m u' = -\lambda x,\quad m v' = -\lambda y - mg,\quad 0 = x^2 + y^2 - L^2$$

로 쓰면, 마지막 구속식에는 $\lambda$ 가 아예 안 나온다. 한 번 미분하면 속도 수준 구속 $xu + yv = 0$ 이 나오고(여기도 $\lambda$ 없음), 두 번 미분하면 가속도 수준에서

$$u^2 + v^2 - \frac{\lambda}{m}\,(x^2+y^2) - g y = 0$$

이 되어 비로소 $\lambda$ 가 등장한다. 여기서 한 번 더 미분해야 $\lambda'$ 를 얻어 완전한 ODE가 되므로 미분지수는 3이다. 각도 $\theta$ 하나로 쓰면 2계 ODE 한 줄로 끝나는 문제인데, 좌표계 선택 하나 때문에 지수 3 DAE가 된 것이다.²

지수의 사촌으로 섭동지수(perturbation index)가 있다. 우변에 작은 잡음 $\delta(t)$ 를 넣었을 때 해의 오차가 $\delta$ 의 몇 계 도함수까지 의존하는가로 정의하며, 수치해석적으로는 이쪽이 더 정직한 지표다. 지수 3이라는 말은 곧 입력 오차가 두 번 미분되어 증폭된다는 뜻이고, 그래서 지수가 높을수록 반올림 오차에 취약하다.

3. 지수 축소와 드리프트[편집]

높은 지수를 그대로 두면 솔버가 고생하니, 구속식을 미분해 지수를 낮추는 것이 지수 축소(index reduction)다. 진자를 가속도 수준 구속으로 바꾸면 지수 1이 되고, BDF든 뭐든 잘 돈다. 문제는 그 대가다.

가속도 수준 구속만 만족시키면 원래의 위치 구속 $x^2+y^2-L^2 = 0$ 은 적분상수만큼 어긋날 자유를 얻는다. 수치오차가 두 번 적분되므로 위반량이 시간에 대해 대략 $t^2$ 로 자란다. 이것이 드리프트(drift)다. 물리 엔진에서 관절이 슬금슬금 늘어나 캐릭터가 분해되는 그 현상이 정확히 이것이고, 강체 동역학·구속 동역학 구현자들의 영원한 숙적이다.

대응책은 셋이다.

• 바움가르테 안정화(Baumgarte, 1972). 가속도 구속 $\ddot g = 0$ 대신 $\ddot g + 2\alpha \dot g + \beta^2 g = 0$ 을 강제한다. 구속 위반에 감쇠 스프링을 매다는 셈이라 위반이 지수적으로 소멸한다. 구현이 세 줄이라 물리 엔진에서 압도적으로 많이 쓰이지만, $\alpha, \beta$ 를 고르는 원리적 기준이 없고 크게 잡으면 계가 강성 방정식이 되며 시간 스텝에 민감해진다. 튜닝 상수 두 개로 바꿔치기한 것이라는 비판이 정당하다.
• GGL 형식(Gear-Gupta-Leimkuhler, 1985). 위치 구속과 속도 구속을 둘 다 남기고, 추가 승수 $\mu$ 를 도입해 $q' = v + G^{\mathsf{T}}\mu$, $Mv' = f - G^{\mathsf{T}}\lambda$, $g(q) = 0$, $G(q)\,v = 0$ 으로 쓴다($G = \partial g/\partial q$). 해석적으로는 $\mu \equiv 0$ 이지만 수치적으로는 $\mu$ 가 위치 구속을 매 스텝 되돌리는 사영 역할을 한다. 지수는 2로 내려가고 드리프트는 원리적으로 사라진다. 튜닝 상수가 없다는 점에서 바움가르테보다 우월하다.
• 좌표 사영. 매 스텝 끝에서 계산된 상태를 구속 다양체 위로 최근접 사영한다. 개념이 명확하고 구속을 기계 정밀도로 만족시키지만, 심플렉틱 적분기의 구조 보존성을 깨뜨릴 수 있다는 것이 대가다.

4. 해법과 초기값[편집]

DAE는 본질적으로 암시적이라 명시적 방법이 존재하지 않는다. 실전 코드는 두 계보다.

• BDF 계열. Petzold의 DASSL과 그 후예인 SUNDIALS의 IDA가 표준이다. 지수 1까지는 이론과 실무가 모두 튼튼하고, 지수 2 헤센베르크 형까지는 가변 스텝에서 제약을 두면 수렴한다. 지수 3부터는 보증이 없어 사실상 지수 축소를 전제한다.
• 암시적 룽게-쿠타. Radau IIA(RADAU5)는 강성 정확(stiffly accurate)해서 대수 구속을 각 스텝 끝에서 정확히 만족시킨다. 지수 3 문제를 직접 다루는 몇 안 되는 범용 코드이며, 대신 오차 추정을 지수별로 스케일링한다 — 지수 2 변수는 $h^{-1}$, 지수 3 변수(승수 $\lambda$)는 $h^{-2}$ 로 오차가 증폭되므로, 승수 성분을 오차 제어에서 빼거나 가중치를 낮추지 않으면 솔버가 스텝을 0으로 줄이며 자살한다.³

차수 저하(order reduction)도 상수처럼 따라온다. 고차 방법이 ODE에서 보이던 차수를 DAE에서 그대로 내주는 경우는 드물고, 특히 대수 변수 쪽에서 한두 차수를 잃는 것이 보통이다.

또 하나의 고전적 함정이 일관된 초기값(consistent initial condition)이다. ODE에서는 $y(0)$ 을 마음대로 줘도 되지만, DAE에서는 $y(0)$ 이 $g = 0$ 뿐 아니라 숨은 구속($\dot g = 0$, $\ddot g = 0$ …)까지 전부 만족해야 한다. 진자를 지수 3으로 풀려면 초기 위치가 원 위에 있어야 하고(위치 구속), 초기 속도가 원의 접선 방향이어야 하며(속도 구속), 초기 $\lambda$ 마저 가속도 구속과 일관돼야 한다. IDA의 IDACalcIC 같은 루틴이 하는 일이 정확히 이 비선형 대수계를 푸는 것이고, 여기서 실패하는 것이 DAE 초보자의 첫 번째 벽이다.

5. 왜 필연적으로 나오는가[편집]

• 회로 해석. SPICE가 쓰는 수정 절점 해석(MNA)은 절점 전압에 KCL을 세우고 전압원·인덕터 전류를 미지수로 추가한다. 커패시터가 붙은 절점은 미분식을, 저항만 붙은 절점은 대수식을 준다 — 정의상 DAE다. 지수는 회로 위상에 따라 결정되며, 전압원과 커패시터로만 이루어진 루프나 전류원과 인덕터로만 이루어진 컷셋이 있으면 지수 2로 올라가고, 그렇지 않으면 지수 1이다. 넷리스트 하나 잘못 쓰면 수렴 실패가 나는 진짜 이유가 이 위상 조건이다.
• 다물체 동역학. 관절로 이어진 강체계를 데카르트 기술자(descriptor) 형식 $M\ddot q = f - G^{\mathsf{T}}\lambda,\ g(q) = 0$ 로 쓰면 무조건 지수 3이다. 최소 좌표계로 가면 DAE를 피할 수 있지만 폐루프 기구에서는 그 좌표계 자체가 존재하지 않는다. 제약 해결기와 위치 기반 동역학은 이 문제를 실시간 예산 안에서 근사적으로 푸는 공학적 타협의 산물이다.
• 비압축성 유동. 압력은 상태변수가 아니라 $\nabla \cdot \mathbf{u} = 0$ 의 라그랑주 승수다. 압력의 시간 미분식이 없다는 사실이 SIMPLE 알고리즘류 압력-속도 연성 절차 전체의 존재 이유이며, 엇갈림 격자와 체커보드 문제도 결국 이 대수 구속의 이산 안정성 조건에서 나온다.

6. 관련 문서[편집]

• 강성 방정식 · 룽게-쿠타법 · 뉴턴-랩슨법 · 자코비안 행렬
• 구속 동역학 · 강체 동역학 · 제약 해결기 · 위치 기반 동역학
• 라그랑주 승수법 · 라그랑주 역학 · 심플렉틱 적분기 · 기하 수치적분
• SPICE · 반도체 소자 시뮬레이션 · SIMPLE 알고리즘 · 엇갈림 격자

7. Footnotes[편집]

1. 그래서 CFD 하는 사람들은 자기가 평생 DAE를 풀고 있었다는 사실을 대개 모른다. “압력 포아송 방정식”이라는 이름 뒤에 숨어 있어서 그렇지, 그건 지수 축소를 손으로 한 결과물이다. ↩

2. 이 예제가 유명한 이유는 “좌표계를 잘못 고르면 문제 난이도가 올라간다”를 세 줄로 보여주기 때문이다. 반대로 말하면 지수 3 DAE를 붙잡고 밤을 새우기 전에 최소 좌표계가 존재하는지 30초만 고민해 보라는 교훈이기도 하다. 물론 폐루프 기구에서는 그 30초가 헛수고다. ↩

3. RADAU5 문서에 “지수 2·3 변수는 오차 제어에서 제외하라”는 취지의 안내가 대문짝만하게 붙어 있는데, 이걸 안 읽고 돌리면 스텝 크기가 $10^{-14}$ 로 수렴하는 장엄한 광경을 보게 된다. 발산이 아니라 수렴이라 더 억울하다. ↩