강성 방정식

수치해석 마지막 수정: 2026-08-02 04:47:19

1. 개요[편집]

강성 방정식 Stiff equation
용어 등장	Curtiss & Hirschfelder (1952)
증상	명시적 적분기가 정확도가 아니라 안정성 때문에 스텝을 줄인다
처방	암시적 방법(후진 오일러 · BDF · Radau IIA · Rosenbrock)
대가	매 스텝 뉴턴 반복 + 자코비안 + LU 분해
이론적 상한	달퀴스트 제2장벽 — A-안정 선형다단계법은 2차가 최대

해는 100초 동안 얌전한 코사인인데, 적분기는 $10^{-6}$ 초짜리 스텝을 100만 번 밟는다. 그리고 아무도 그게 왜 필요한지 물어보지 않는다.

강성 방정식(stiff equation)은 해 자체는 매끄럽게 변하는데도 명시적 시간 적분기가 정확도 요구가 아니라 수치적 안정성 요구 때문에 터무니없이 작은 시간 스텝을 강요당하는 미분방정식을 말한다. 1952년 커티스와 히르슈펠더가 화학반응 계산 도중 이 현상에 이름을 붙였고¹, 이후 70년 동안 수치 상미분방정식 이론의 절반은 사실상 이 한 단어를 중심으로 굴러갔다.

주의할 점은 강성이 방정식만의 고유한 성질이 아니라는 것이다. 강성은 (문제) × (적분기) × (요구 정확도) × (적분 구간)의 관계에서 나타나는 성질이다. 같은 방정식이 명시적 RK4에게는 지옥이고 후진 오일러에게는 아무 일도 아니며, $10^{-3}$ 의 정확도를 원할 때는 강성이지만 $10^{-9}$ 을 원할 때는 어차피 스텝을 잘게 써야 하므로 강성이 문제되지 않는다.²

2. 정의가 미묘한 이유[편집]

교과서가 흔히 드는 정의는 강성비(stiffness ratio)다. 선형화된 계 $y' = Jy$ 의 자코비안 행렬 고유값 $\lambda_i$ 가 모두 좌반평면에 있을 때

S = \frac{\max_i \lvert \mathrm{Re}\,\lambda_i \rvert}{\min_i \lvert \mathrm{Re}\,\lambda_i \rvert} \gg 1

이면 강성이라는 것이다. 직관은 맞다 — 빠른 모드가 이미 소멸했는데도 그 시간척도가 스텝 크기를 계속 인질로 잡는 상황이니까. 문제는 이 정의가 필요조건도 충분조건도 아니라는 데 있다.

스칼라도 강성일 수 있다. $y' = -\lambda(y - \cos t) - \sin t$ 는 고유값이 $-\lambda$ 하나뿐이라 강성비가 정의되지 않지만, 해석해 $y = \cos t$ 는 $\lambda$ 와 무관하게 매끄러운데도 $\lambda = 10^4$ 이면 명시적 방법은 죽는다.
강성비가 커도 강성이 아닐 수 있다. 적분 구간이 빠른 시간척도 수준으로 짧으면(초기 과도 구간만 보는 경우) 어차피 그 스케일을 분해해야 하므로 문제될 게 없다.
비선형계에서는 고유값 자체가 시간에 따라 변한다. 어떤 구간에서만 강성이 나타나는 문제가 오히려 일반적이다.

그래서 실무적으로 살아남은 것은 커티스-히르슈펠더의 서술적 정의다. “명시적 방법이 안정성 때문에 극히 작은 스텝을 쓰도록 강요당하는 문제.” 정량적이지는 않지만 현장에서 판별은 쉽다. 스텝을 절반으로 줄였을 때 해가 눈에 띄게 안 바뀌는데도 그렇게 안 하면 발산한다면, 그게 강성이다.

3. 시험 방정식과 안정영역[편집]

이 모든 논의의 기준점은 달퀴스트의 시험 방정식 $y' = \lambda y$ ( $\mathrm{Re}\,\lambda < 0$ )다. 어떤 일단계 방법을 여기에 적용하면 반드시 $y_{n+1} = R(z)\,y_n$ ( $z = \lambda\Delta t$ ) 꼴이 되고, $R(z)$ 를 안정함수(stability function)라 부른다. 정확한 해는 $e^{z}$ 이므로 $R(z)$ 는 $e^z$ 의 유리근사인 셈이다.

방법	$R(z)$	$R(\infty)$	성질
전진 오일러	$1 + z$	발산	조건부 안정
RK4	$1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24}$	발산	실축 안정한계 $z \approx -2.785$
후진 오일러	$(1-z)^{-1}$	$0$	A-안정, L-안정
사다리꼴	$\frac{1+z/2}{1-z/2}$	$-1$	A-안정, L-안정 아님

A-안정성은 좌반평면 전체( $\mathrm{Re}\,z < 0$ )에서 $\lvert R(z) \rvert \le 1$ 인 성질이다. 즉 물리적으로 감쇠하는 모드는 스텝 크기와 무관하게 수치적으로도 감쇠한다. L-안정성은 여기에 $R(\infty) = 0$ 을 추가로 요구한다. 이 차이가 실무에서 뼈아프게 드러나는 곳이 사다리꼴(크랭크-니콜슨법)이다. $\lambda\Delta t \to -\infty$ 에서 $R \to -1$ 이므로, 빠른 모드가 감쇠하는 대신 매 스텝 부호를 바꾸며 진동한다. 감쇠는 하지만 아주 천천히 한다. SPICE 회로 해석에서 스위칭 직후 전압 파형에 톱니가 끼는 그 유명한 “사다리꼴 링잉”이 정확히 이것이고, 그래서 SPICE 계열은 대개 기어(BDF) 적분기를 대안으로 준비해 둔다.

y' = -λ(y - cos t) - sin t 를 같은 Δt로 명시적 RK4와 사다리꼴로 동시에 적분한다. 해석해는 λ와 무관하게 y = cos t 인데도 RK4는 λΔt가 안정영역을 벗어나는 순간 터진다. 아래 패널의 안정영역은 그려둔 그림이 아니라 복소평면 격자에서 |R(z)|를 실제로 평가해 채운 것이고, 빨간 점이 현재 λΔt다.

여기에 두 개의 유명한 벽이 있다. 달퀴스트 제1장벽은 영-안정(zero-stable)인 $k$ 단계 선형다단계법의 차수가 $k$ 가 홀수면 $k+1$ , 짝수면 $k+2$ 를 넘을 수 없다고 말한다. 달퀴스트 제2장벽은 훨씬 잔인하다. A-안정인 선형다단계법은 2차 정확도를 넘을 수 없고, 그 2차 중에서 오차상수가 가장 작은 것이 사다리꼴이다. 강성 문제에서 고차 정확도와 무조건 안정성을 동시에 원한다면 다단계법을 버리고 암시적 룽게-쿠타법으로 가야 한다는 뜻이다.

4. 방법들[편집]

후진 오일러. 1차지만 L-안정하고 무엇에도 안 터진다. 정확도를 버리고 안정성을 사는 극단. 수치 감쇠가 커서 진동 문제에는 못 쓴다.
사다리꼴 / 크랭크-니콜슨법. 2차 A-안정. 열전달 해석과 대류-확산 방정식의 국룰이지만 위에서 본 진동 때문에 급격한 초기조건에는 주의가 필요하다.
BDF (후진차분식). 기어(Gear)가 대중화한 다단계 계열로 1~~6차까지 존재한다. 7차부터는 영-안정성을 잃어서 아예 쓸 수 없다.³ BDF1(=후진 오일러)과 BDF2만 A-안정이고, BDF3~~6은 A( $\alpha$ )-안정으로 각각 $\alpha \approx 86.0^\circ,\ 73.4^\circ,\ 51.8^\circ,\ 17.8^\circ$ 까지만 커버한다. 차수를 올릴수록 안정영역의 부채꼴이 좁아지므로, 허축 근처 고유값(진동 모드)이 있으면 저차로 내려가야 한다.
암시적 룽게-쿠타. Radau IIA는 $s$ 단으로 $2s-1$ 차이면서 L-안정이고 강성 정확(stiffly accurate)하다. 3단 5차짜리가 유명한 RADAU5. SDIRK/ESDIRK는 계수행렬을 하삼각으로 잡아 단마다 순차적으로 풀 수 있게 만든 절충안으로, 큰 계에서 비용이 훨씬 싸다.
Rosenbrock(선형 암시적). 뉴턴 반복을 아예 포기하고, 자코비안을 한 번 만들어 각 단을 선형 방정식 풀이 한 번으로 끝낸다. 반복 수렴 실패가 없다는 것이 최대 장점이고, 대신 자코비안 정확도에 성능이 민감하다. 부정확한 자코비안을 허용하도록 일반화한 것이 W-방법.
지수 적분기. 강성의 원인인 선형부 $L$ 을 $e^{\Delta t L}$ 과 $\varphi$ -함수로 정확히 처리하고 비선형부만 명시적으로 다룬다. 행렬 지수를 크리로프 부분공간법으로 근사하는 것이 실용 열쇠다.

5. 암시적 방법의 청구서[편집]

암시적이라는 말은 매 스텝 $F(y_{n+1}) = 0$ 을 푼다는 뜻이고, 이걸 푸는 것은 결국 뉴턴-랩슨법이다. 각 뉴턴 반복은 $\left(I - \gamma\Delta t\,J\right)\Delta y = -F$ 를 풀어야 하므로, 자코비안 행렬 구성 + LU 분해(혹은 희소행렬 분해) 비용이 스텝마다 붙는다. 대형 계에서는 이 한 줄이 전체 계산 시간의 80~90%를 먹는다. 그래서 실전 코드는 온갖 꼼수를 쓴다.

자코비안을 매 스텝이 아니라 수렴이 나빠질 때만 재계산한다(modified Newton).
같은 $\Delta t$ 와 $\gamma$ 를 유지하는 동안 LU 인수를 재사용한다. BDF 코드가 스텝 크기를 잘 안 바꾸려 하는 진짜 이유가 정확도가 아니라 이것이다.
직접 분해 대신 크리로프 반복 + 전처리기를 쓴다(Newton-Krylov). 이때 조건수가 곧 반복 횟수를 결정한다.
자코비안이 나쁘면 뉴턴이 안 도는데, 그건 곧 스텝 축소로 이어진다. 무조건 안정하다고 무조건 큰 스텝을 쓸 수 있는 게 아니다. 안정성은 공짜여도 수렴성은 아니다.

6. 강성은 어디서 나오나[편집]

반응속도론과 연소 시뮬레이션. 강성의 원조. 라디칼 재결합은 $10^{-9}$ s, 전체 연소는 $10^{0}$ s 스케일이라 강성비가 $10^{9}$ 을 우습게 넘는다. 표준 벤치마크인 로버트슨 문제(속도상수 $0.04,\ 3\times10^{7},\ 10^{4}$ 인 3화학종 계)는 모든 강성 솔버가 통과해야 하는 신고식이다.⁴
확산항. 대류-확산 방정식이나 열전달 해석에서 확산을 유한차분법으로 이산화하면 최대 고유값이 $-4\alpha/\Delta x^2$ 수준이라, 격자를 반으로 줄일 때마다 강성비가 4배로 뛴다. 명시적 스텝 제한 $\Delta t \le \Delta x^2/(2\alpha)$ 가 폰 노이만 안정성 해석에서 나오는 그 조건이다.
구조의 고주파 모드. 명시적 동해석이 마이크로초 스텝을 쓰는 것은 CFL 조건 — 정확히는 가장 작은 요소를 응력파가 가로지르는 시간 — 때문이다. 강판 1 mm 요소면 $\Delta t \approx 0.2\ \mu\mathrm{s}$ . 충돌 해석이 아니라 준정적 성형 해석을 명시적으로 돌릴 때 질량 스케일링이라는 편법이 등장하는 이유가 이것이다.
SPICE 회로. 기생 커패시턴스의 시상수는 ps, 관심 있는 스위칭 파형은 ns~μs.
흥분성 매질과 심장/신경 모형. 활동전위의 빠른 탈분극과 느린 회복이 한 계에 공존한다.

전부 암시적으로 풀 필요는 없다. IMEX(implicit-explicit) 전략은 연산자 분리로 강성항(확산, 화학원)만 암시적으로, 비강성항(대류)은 명시적으로 처리한다. 확산은 무조건 안정하게, 대류는 싼 명시적 계산으로 — 이 조합이 대류-확산-반응 문제의 현실적인 타협점이다. 분할오차가 추가로 생긴다는 것이 대가이고, 그 크기가 두 연산자의 교환자에 비례한다는 점은 연산자 분리 문서에 있다.

7. 관련 문서[편집]

커티스와 히르슈펠더의 1952년 PNAS 논문 제목이 아예 Integration of stiff equations다. “stiff”라는 단어 선택은 구조역학의 강성(stiffness)과는 무관하며, 굳이 따지면 “뻣뻣해서 잘 안 움직이는” 쪽의 어감이다. 한국어 번역이 강성행렬의 그것과 겹쳐 버린 건 순전한 사고다. ↩
그래서 “이 방정식은 강성인가요?”라는 질문에 정색하고 답하려면 최소한 “어떤 적분기로, 어느 정확도로, 얼마나 긴 구간을”까지 되물어야 한다. 물론 현장에서 이렇게 되묻는 순간 상대는 그냥 암시적 솔버를 켜고 자리를 뜬다. ↩
BDF7 이상이 못 쓰이는 것은 안정영역이 좁아서가 아니라 영-안정성 자체가 깨지기 때문이다. 즉 $\Delta t \to 0$ 에서도 수렴하지 않는다. 수치해석에서 가장 기본적인 요구조건이 무너지는 것이라, 아무리 급해도 방법이 없다. ↩
로버트슨 문제를 명시적 RK4로 풀어 보는 것은 수치해석 수업의 통과의례다. 결과는 대개 “노트북 팬이 이륙 준비를 하고, 30초 뒤 NaN이 뜬다”로 요약된다. ↩