강성 방정식

1. 개요

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

해는 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$이면 명시적 방법은 죽는다. 이 예제에는 이름이 있다 — 프로테로-로빈슨 문제(1974)로, 강성 정확도와 차수 저하를 논할 때의 표준 시험대다.
• 강성비가 커도 강성이 아닐 수 있다. 적분 구간이 빠른 시간척도 수준으로 짧으면(초기 과도 구간만 보는 경우) 어차피 그 스케일을 분해해야 하므로 문제될 게 없다.
• 비선형계에서는 고유값 자체가 시간에 따라 변한다. 어떤 구간에서만 강성이 나타나는 문제가 오히려 일반적이다.
• 자코비안이 특이하면 분모가 0이다. 보존량이 있는 계에서는 이게 예외가 아니라 기본값이다(바로 아래 검산).

그래서 실무적으로 살아남은 것은 커티스-히르슈펠더의 서술적 정의다. “명시적 방법이 안정성 때문에 극히 작은 스텝을 쓰도록 강요당하는 문제.” 정량적이지는 않지만 현장에서 판별은 쉽다. 스텝을 절반으로 줄였을 때 해가 눈에 띄게 안 바뀌는데도 그렇게 안 하면 발산한다면, 그게 강성이다. 2015년에 쇠데르린드·제이·칼보가 “Stiffness 1952–2012: 정의를 찾아 60년”이라는 제목으로 이 분야를 한 번 정리했는데, 결론은 아직 모두가 동의하는 정의는 없다는 것이었다.³

2.1. 검산: 로버트슨 문제의 강성비는 아예 정의되지 않는다

강성 솔버의 표준 신고식인 로버트슨 문제(Robertson, 1966)로 위 주장을 직접 확인했다. 3화학종 자기촉매 반응

$$\mathrm{A} \xrightarrow{k_1} \mathrm{B},\qquad \mathrm{B} + \mathrm{B} \xrightarrow{k_2} \mathrm{B} + \mathrm{C},\qquad \mathrm{B} + \mathrm{C} \xrightarrow{k_3} \mathrm{A} + \mathrm{C}$$

에 $k_1 = 0.04$, $k_2 = 3\times10^7$, $k_3 = 10^4$, $y(0) = (1, 0, 0)$을 넣고 $t \in [0,\ 4\times10^{10}]$을 후진 오일러 + 뉴턴으로 적분하면서 매 스텝 자코비안의 고유값을 뽑았다.

$t$	$y_1$	$y_2$	$y_3$	$\lvert\lambda\rvert$ 세 개
$10^{-5}$	1.0000	4.08e-7	1.79e-11	0, 4.41e-2, 2.45e+1
$1$	0.9654	3.06e-5	3.46e-2	0, 2.91e-1, 2.18e+3
$10^{3}$	0.3382	2.03e-6	6.62e-1	0, 1.09e-3, 6.74e+3
$4\times10^{10}$	5.3e-8	2.10e-13	1.0000	0, 5.05e-11, 1.00e+4

고유값 하나가 모든 시각에서 정확히 0이다. 질량보존 $y_1+y_2+y_3 = 1$ 때문에 자코비안의 열 합이 0이고(실측 열 합 $= 0,\ 5.6\times10^{-13},\ 0$), 따라서 $J$는 구조적으로 특이하다. 위 정의의 분모가 0이니 강성비는 무한이거나 정의 불가다. 이 바닥에서 가장 유명한 강성 문제가 가장 유명한 강성 정의를 통과하지 못한다.

0이 아닌 고유값만 보면 궤적 전체에서 $\lvert\lambda\rvert$가 $5.05\times10^{-11}$부터 $1.00\times10^{4}$까지 퍼져 비가 $2.0\times10^{14}$다. 같은 구간을 안정성을 유지하며 명시적 RK4로 밟으려면 $\Delta t \le 2.785/10^4$이므로 약 $1.4\times10^{14}$ 스텝이 필요하다. 반면 스텝을 매번 1.05배로만 키우는 조악한 후진 오일러 코드가 1136 스텝으로 끝냈고 $\sum y_i = 1$이 소수점 12자리까지 유지됐다. 이 두 숫자의 거리가 강성이라는 단어의 전부다.

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) 적분기를 대안으로 준비해 둔다.

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

3.1. 검산: 안정한계와 BDF 의 A($\alpha$) 각을 직접 재보기

교과서 숫자를 그대로 믿지 않고 안정함수를 복소평면에서 평가해 다시 뽑았다. RK4는 이분법으로 $\lvert R \rvert = 1$이 되는 점을 찾았고, BDF는 $t=0,\dots,k$ 노드의 보간 도함수로 계수를 만든 뒤 경계 궤적 $z(\theta) = \rho(e^{i\theta})/e^{ik\theta}$ 를 200만 점에서 훑어 음실축과의 최소각을 구했다.

방법	실측값	비고
전진 오일러 실축 한계	$-2$	$\lvert 1+z\rvert \le 1$ 원
RK4 실축 한계	$-2.7852935634$	흔히 인용되는 $-2.785$ 확인
RK4 허축 한계	$\pm 2.82842712\,i$	$2\sqrt2 = 2.82842712$ 와 일치
BDF1 · BDF2	$180^\circ$ · $90^\circ$	A-안정
BDF3	$86.03^\circ$	A($\alpha$)-안정
BDF4	$73.35^\circ$
BDF5	$51.84^\circ$
BDF6	$17.84^\circ$	실용 한계

RK4의 허축 한계가 정확히 $2\sqrt2$라는 점이 재밌다 — 감쇠 없는 순수 진동 모드에 대해서는 $\Delta t \le 2\sqrt2/\omega$가 상한이고, 이것이 명시적 동해석이 CFL 조건에서 만나는 그 숫자와 같은 종류의 제약이다.

BDF7이 못 쓰이는 이유도 같은 코드로 확인했다. $\rho(w)$의 근 중 최대 모듈러스가 BDF6에서는 정확히 $1.000000$인데 BDF7에서는 $1.022218$, BDF8에서는 $1.183870$이다. 단위원을 벗어나는 순간 영-안정성이 깨진다 — 안정영역이 좁아지는 문제가 아니라 $\Delta t \to 0$에서도 수렴하지 않는다는 뜻이다.⁵

4. 방법들

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

4.1. 검산: 사다리꼴 링잉과 L-안정성의 차이

$y' = \lambda y$에 $\lambda\Delta t = -200$을 넣고 두 방법을 돌렸다. 사다리꼴의 증폭계수는 $R = -0.980198$이므로 매 스텝 부호가 뒤집힌다. 여섯 스텝 값이 $-0.980,\ +0.961,\ -0.942,\ +0.923,\ -0.905,\ +0.887$로 진동하고, $\lvert y\rvert < 0.01$까지 줄어드는 데 231 스텝이 걸린다. 정확해는 한 스텝에 $e^{-200}$으로 사라지는 모드이고, 후진 오일러는 $R = 4.975\times10^{-3}$으로 한 스텝에 그 자리를 떠난다. A-안정이라는 도장만으로는 이 차이를 구분할 수 없다는 것이 L-안정성 개념이 필요한 이유다.

4.2. 검산: 안정성이 정확도를 인질로 잡는 순간

프로테로-로빈슨 문제 $y' = -\lambda(y - \cos t) - \sin t$, $y(0) = 1$, $\lambda = 10^4$, $t \in [0, 10]$. 해석해는 $y = \cos t$로 $\lambda$와 무관하다.

방법	$\Delta t$	$\lambda\Delta t$	최대 오차	스텝 수
RK4	2.757e-4	2.757	3.99e-7	36,270
RK4	2.813e-4	2.813	발산 (17,216스텝에서 NaN)	—
RK4	1e-2	100	발산 (48스텝에서 NaN)	—
사다리꼴	1e-2	100	8.33e-10	1,000
사다리꼴	1e-1	1,000	8.34e-8	100
후진 오일러	1e-1	1,000	5.00e-6	100

RK4의 생사가 오차가 아니라 $\lambda\Delta t$가 $2.785$를 넘는지에만 달려 있다. $\lambda\Delta t = 2.757$에서 오차가 이미 $4\times10^{-7}$인데도 스텝을 2% 키우면 그냥 죽는다. 같은 정확도를 사다리꼴은 $\Delta t = 0.1$, 즉 362배 적은 스텝으로 낸다. “정확도 때문이 아니라 안정성 때문에 스텝을 줄인다”는 문장의 실측판이 이 표다.

5. 암시적 방법의 청구서

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

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

암시적 방법에는 정확도 쪽 청구서도 따로 온다. 강성 구간에서 룽게-쿠타법의 실제 관측 차수가 이론 차수보다 낮게 나오는 차수 저하(order reduction) 현상이 그것이다. 프로테로와 로빈슨이 1974년에 이 문제를 정면으로 다루면서 강성 차수(stiff order)라는 개념을 도입했고, Radau IIA 계열이 “강성 정확”하다고 불리는 근거도 여기에 있다. 미분대수방정식은 이 이야기를 $\lambda \to \infty$ 극한까지 밀고 간 경우다.

6. 강성은 어디서 나오나

• 반응속도론과 연소 시뮬레이션. 강성의 원조. 라디칼 재결합은 $10^{-9}$ s, 전체 연소는 $10^{0}$ s 스케일이라 시간척도 비가 $10^{9}$을 우습게 넘는다. 위에서 검산한 로버트슨 문제가 이 계열의 최소 재현 예제다.⁸
• 확산항. 대류-확산 방정식이나 열전달 해석에서 확산을 유한차분법으로 이산화하면 최대 고유값이 $-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. 관련 문서

• 크랭크-니콜슨법 · 룽게-쿠타법 · 암시적 룽게-쿠타법 · 지수 적분기 · ADI 법
• 폰 노이만 안정성 해석 · CFL 조건 · 수렴성 · 절단오차
• 뉴턴-랩슨법 · 자코비안 행렬 · LU 분해 · 희소행렬 · 조건수 · 전처리기 · 크리로프 부분공간법
• 연산자 분리 · 대류-확산 방정식 · 유한차분법 · 미분대수방정식
• 반응속도론 · 연소 시뮬레이션 · 흥분성 매질
• 명시적 동해석 · 열전달 해석 · SPICE

8. Footnotes

1. 커티스와 히르슈펠더의 1952년 PNAS 논문 제목이 아예 Integration of stiff equations다. “stiff”라는 단어 선택은 구조역학의 강성(stiffness)과는 무관하며, 굳이 따지면 “뻣뻣해서 잘 안 움직이는” 쪽의 어감이다. 한국어 번역이 강성행렬의 그것과 겹쳐 버린 건 순전한 사고다. ↩

2. 그래서 “이 방정식은 강성인가요?”라는 질문에 정색하고 답하려면 최소한 “어떤 적분기로, 어느 정확도로, 얼마나 긴 구간을”까지 되물어야 한다. 물론 현장에서 이렇게 되묻는 순간 상대는 그냥 암시적 솔버를 켜고 자리를 뜬다. ↩

3. 논문의 결론 중 하나가 “강성은 문제의 성질이 아니라 문제-방법-정확도 쌍의 성질”이라는, 이미 1952년에 커티스와 히르슈펠더가 서술적으로 말한 것이었다. 60년 동안 정의를 찾다가 원점 확인으로 끝나는 전개는 수치해석 특유의 정직함이다. ↩

4. 이 두 장벽이 다 달퀴스트 한 사람에게서 나왔다는 사실이 이 분야의 지분 구조를 요약한다. 제2장벽의 원 논문 제목은 겸손하게도 “선형다단계법의 특수한 안정성 문제”다. ↩

5. 위 검산의 $\rho(w)$ 근 모듈러스가 그 증거다. BDF7의 $1.0222$ 는 “조금 넘었다”가 아니라 스텝마다 2.2%씩 증폭되는 기생 모드가 있다는 뜻이고, 수치해석에서 가장 기본적인 요구조건이 무너지는 것이라 아무리 급해도 방법이 없다. ↩

6. 그 코드가 CACM Algorithm 407 DIFSUB(1971)이고, 오늘날 MATLAB ode15s·SUNDIALS CVODE·scipy.integrate.BDF 로 이어지는 계보의 출발점이다. Fortran 서브루틴 한 개가 반세기짜리 표준이 된 드문 사례. ↩

7. 그래서 암시적 솔버의 로그를 보면 “스텝 크기를 바꿀 이유가 충분히 크지 않으면 안 바꾼다”는 히스테리시스가 들어 있다. 정확도 제어가 아니라 LU 인수를 아끼려는 구두쇠 로직이다. ↩

8. 로버트슨 문제를 명시적 RK4로 풀어 보는 것은 수치해석 수업의 통과의례다. 결과는 대개 “노트북 팬이 이륙 준비를 하고, 30초 뒤 NaN이 뜬다”로 요약된다. 위 검산의 $1.4\times10^{14}$ 스텝이 그 팬 소리의 정량적 근거다. ↩