지수 적분기

1. 개요[편집]

지수 적분기 Exponential integrator
대상	반선형 강성계 $u' = Lu + N(u)$
핵심 도구	$\varphi$-함수, $\varphi_1(z) = (e^z-1)/z$
대표 기법	지수 오일러, ETDRK4, 지수 로젠브록
병목	행렬 지수-벡터 곱 $\varphi_k(hL)\,v$
최적 조건	$L$ 이 푸리에 공간에서 대각인 주기 문제

강성의 원인이 선형항이라면, 선형항을 근사하지 말고 그냥 정확히 풀면 되지 않나? — 사실 이 생각은 1960년대부터 있었다. 문제는 $e^{hL}$ 을 계산할 방법이 없었다는 것이다.

지수 적분기(exponential integrator)는 미분방정식의 선형 강성항을 행렬 지수로 정확히 적분하고 비선형항만 다항식으로 근사하는 시간 적분 계열이다. 대상은 반선형(semilinear) 형태

$$u'(t) = L\,u(t) + N(u(t))$$

이며, $L$ 이 큰 음의 고유값을 갖는 강성 부분(확산, 초확산, 분산항)이고 $N$ 은 비교적 온순한 비선형 부분(대류, 반응)인 상황을 전제한다. 강성 방정식의 표준 처방인 암시적 방법이 매 스텝 뉴턴-랩슨법과 LU 분해를 청구하는 반면, 지수 적분기는 뉴턴 반복 없이 무조건 안정성을 얻는다. 대신 청구서가 다른 곳으로 온다 — 행렬 지수-벡터 곱.

출발점은 상수변화법(variation of constants)의 정확한 항등식이다.

$$u(t_n + h) = e^{hL} u_n + \int_0^h e^{(h-\tau)L}\, N\!\left(u(t_n+\tau)\right) d\tau$$

여기까지는 근사가 하나도 없다. 남은 자유도는 적분 안의 $N(u(t_n+\tau))$ 를 어떻게 다항식으로 보간할 것인가 하나뿐이고, 그 선택이 곧 방법의 이름이 된다.

2. φ 함수[편집]

적분 안의 $N$ 을 $\tau$ 에 대한 다항식으로 두면 $\int_0^h e^{(h-\tau)L}\tau^{k}d\tau$ 꼴의 적분이 나오고, 이것들이 $\varphi$ 함수로 정리된다. 정의는

$$\varphi_0(z) = e^{z}, \qquad \varphi_{k+1}(z) = \frac{\varphi_k(z) - 1/k!}{z}$$

이고 처음 몇 개는 $\varphi_1(z) = (e^z - 1)/z$, $\varphi_2(z) = (e^z - 1 - z)/z^2$ 다. 적분 표현 $\varphi_k(z) = \frac{1}{(k-1)!}\int_0^1 e^{(1-\theta)z}\theta^{k-1}d\theta$ 에서 보이듯 $z = 0$ 에서 특이점이 없고 $\varphi_k(0) = 1/k!$ 다. 그런데 부동소수점에서는 $|z|$ 가 작을 때 정의식이 재앙이다. $e^z - 1$ 을 $z$ 로 나누는 순간 자릿수 소거가 일어나 $|z| \sim 10^{-8}$ 근처에서 유효숫자가 반토막 난다.¹

$N$ 을 $\tau$ 에 대해 상수 $N(u_n)$ 으로 두는 가장 단순한 선택이 지수 오일러다.

$$u_{n+1} = e^{hL} u_n + h\,\varphi_1(hL)\,N(u_n)$$

일반적으로 1차 정확도지만, $N$ 이 상수인 비동차 선형계에서는 정확한 해를 준다. 이 성질은 생각보다 중요하다 — 정상해(steady state)를 정확히 재현하므로 장시간 적분에서 인위적 표류가 없다. 여기서 $L$ 을 고정 행렬이 아니라 매 스텝의 야코비안 $J_n = \partial(Lu+N)/\partial u|_{u_n}$ 으로 바꾼 것이 지수 로젠브록-오일러이고, 이 경우 차수가 2로 올라간다. 자코비안을 매 스텝 새로 쓴다는 점에서 준-뉴턴법이 아니라 로젠브록 계열의 사고방식이다.

3. ETDRK4와 차수 조건[편집]

실무에서 가장 널리 쓰이는 4차 방법이 콕스와 매튜스(2002)의 ETDRK4(Exponential Time Differencing, 4th-order Runge-Kutta)다. 고전 RK4와 같은 4단 구조를 갖되, 각 단의 계수가 스칼라가 아니라 $\varphi_k(hL)$, $\varphi_k(hL/2)$ 행렬 함수다. 카삼과 트레페텐(2005)이 이 방법을 스펙트럴 공간 구현과 함께 정리하면서 사실상 표준이 됐다.

그들의 기여 중 하나가 위에서 말한 소거 문제의 해법이다. $\varphi_k$ 를 정의식으로 계산하는 대신, $z$ 를 중심으로 하는 복소평면상의 작은 원 위에서 등간격 표본(예: 32점)의 평균을 취해 코시 적분으로 값을 얻는다. 실수축상의 소거가 복소 등고선 위에서는 일어나지 않으므로 기계 정밀도 근처의 정확도가 회복된다. 코드 열 줄짜리 요령인데 이게 없으면 ETDRK4는 소형 $|hL|$ 영역에서 4차 수렴을 보여주지 못한다.

또 하나 짚어야 할 것이 강성 차수 조건(stiff order conditions)이다. 고전 룽게-쿠타의 차수 조건은 $hL$ 이 유계라는 가정 아래 유도된 것이라, $\|hL\| \to \infty$ 인 강성 극한에서는 실제 차수가 떨어진다. 호흐브루크와 오스터만이 2010년 Acta Numerica 리뷰에서 정리한 강성 차수 조건은 이 극한에서도 차수를 보증하도록 훨씬 빡빡한 조건을 요구하며, 이 관점에서 보면 ETDRK4 계열의 일부는 문제에 따라 3차 정도로 저하된다. 크로그스타드의 변형(ETDRK4-B)이 이 저하를 완화하려는 시도 중 하나다.

지수 적분기의 사촌으로 적분인자법(integrating factor, IF)이 있다. $v = e^{-tL}u$ 로 치환한 뒤 $v$ 에 보통의 RK를 돌리는 방식이라 구현이 더 쉽다. 하지만 치환된 방정식의 해가 시간에 따라 격렬하게 진동하는 성분을 갖게 되어 오차상수가 훨씬 크고, 정상해를 정확히 재현하지 못한다. 같은 문제에서 IF-RK4와 ETDRK4를 비교하면 후자가 수십 배 정확한 것이 보통이다.

4. 행렬 지수-벡터 곱[편집]

지수 적분기의 성패는 결국 $\varphi_k(hL)\,v$ 를 얼마나 싸게 구하느냐로 결정된다. $L$ 을 명시적으로 지수화하는 것($\exp(hL)$ 행렬 자체를 구하는 것)은 $N \times N$ 밀집 행렬이 되므로 큰 문제에서는 논외다. 실전 경로는 셋이다.

• 대각화가 공짜인 경우. 주기 경계 + 상수계수 미분연산자면 고속 푸리에 변환 뒤 $L$ 이 대각이므로 $\varphi_k(hL)$ 은 성분별 스칼라 함수 평가다. 계수를 시작 전에 한 번 계산해 캐싱하면 이후 스텝 비용은 FFT 두 번뿐. 스펙트럴 방법과 지수 적분기가 세트로 다니는 이유이며, ETDRK4의 화려한 성능 그래프는 대부분 이 조건에서 찍힌 것이다.
• 크릴로프 부분공간. 일반적인 희소 $L$ 에는 아놀디 알고리즘으로 $m$ 차원 크리로프 부분공간법 $\mathcal{K}_m(L, v)$ 를 만들고 $\varphi_k(hL)v \approx \|v\|\,V_m \varphi_k(hH_m) e_1$ 로 근사한다. $H_m$ 은 $m \times m$ ($m \sim 30$) 소형 행렬이라 밀집 알고리즘으로 처리 가능하다. 여러 $\varphi_k$ 를 한꺼번에 얻는 확대행렬 요령(블록 상삼각 행렬 하나의 지수로 $\varphi_0 \dots \varphi_p$ 를 동시에 뽑는 트릭)이 표준으로 쓰이며, EXPOKIT·phipm·KIOPS 같은 구현이 이를 담고 있다.
• 유리근사. 파데 근사 + 스케일링·제곱법(expm의 고전 알고리즘), 체비셰프 전개, 음의 실축 위 최량 유리근사인 카라테오도리-페예르(CF) 근사 등. CF의 (14,14)차 근사는 음의 실축에서 오차가 $10^{-14}$ 수준이라 부분분수 전개 형태로 쓰면 선형계 14개 풀이로 $\varphi$ 를 얻는다 — 실질적으로 암시적 방법과 같은 작업량이 된다. 유리 크릴로프와 이동-역변환 아놀디도 같은 계열의 절충이다.

주의할 점은 크릴로프 수렴 자체가 $L$ 의 스펙트럼 분포에 강하게 의존한다는 것이다. $L$ 이 대칭이면 예측 가능하지만, 대류가 섞여 비정규가 되면 의사스펙트럼이 부풀어 오르고 수렴 예측이 무의미해진다.

5. 이득과 한계[편집]

지수 적분기가 빛나는 표준 시험대가 쿠라모토-시바신스키 방정식 $u_t = -u u_x - u_{xx} - u_{xxxx}$ 다. 푸리에 공간에서 $L$ 의 기호가 $k^2 - k^4$ 이라 최고파수 성분의 고유값이 $-k_{\max}^4$ 로 격자수의 4제곱에 비례한다. 명시적 방법은 $\Delta t \sim \Delta x^4$ 라는 절망적 제한을 받고, 크랭크-니콜슨류는 안정하지만 4차 정확도를 못 낸다. ETDRK4는 이 항을 정확히 적분하므로 $\Delta t$ 를 물리적 시간척도에 맞춰 잡을 수 있고, 같은 정확도에서 명시적 RK4 대비 수십~수백 배 빠르다. 앨런-칸 방정식 $u_t = \epsilon^2 \Delta u + u - u^3$, 그레이-스콧류 반응-확산(튜링 패턴), 복소 긴즈부르크-란다우, 주기 상자 안의 비압축성 유동도 같은 구조다.

한계도 분명하다.

• $L$ 이 대각화하기 어려우면 이점이 급격히 준다. 복잡한 형상의 유한요소법 이산화, 변수계수, 비주기 경계 — 이 경우 매 스텝 크릴로프를 돌려야 하고 그 비용은 암시적 솔버의 전처리기가 붙은 선형계 풀이와 비슷해진다. “지수 적분기가 빠르다”는 벤치마크가 대부분 스펙트럴 설정에서 나온다는 사실은 항상 기억해야 한다.
• 비선형항이 강성이면 소용없다. 강성이 $N$ 쪽(예: 빠른 화학반응)에 있으면 분해 자체가 무의미하다. 이때는 지수 로젠브록처럼 야코비안 전체를 $L$ 로 삼는 변형이 필요하고, 그러면 $L$ 이 매 스텝 바뀌어 캐싱 이점이 사라진다.
• 메모리. 크릴로프 기저 $V_m$ 을 저장해야 하고, 다단계 지수 방법은 여러 $\varphi$ 벡터를 동시에 들고 있어야 한다.
• 차수 저하와 오차 추정. 위에서 말한 강성 차수 조건 문제 때문에, 관측 차수가 이론값보다 낮게 나오는 일이 흔하다. 적응 스텝 제어를 붙이려면 임베디드 쌍이 필요한데 선택지가 암시적 RK만큼 풍부하지 않다.

요약하면 주기 경계 + 스펙트럴 + 반선형 강성 PDE 라는 삼박자가 맞을 때는 사실상 대체재가 없는 방법이고, 그 밖에서는 연산자 분리나 IMEX와 정직하게 성능을 겨뤄야 하는 후보 중 하나다.²³

6. 관련 문서[편집]

• 강성 방정식 · 룽게-쿠타법 · 크랭크-니콜슨법 · 연산자 분리
• 크리로프 부분공간법 · 아놀디 알고리즘 · 의사스펙트럼
• 스펙트럴 방법 · 고속 푸리에 변환 · 편미분방정식
• 튜링 패턴 · 벨루소프-자보틴스키 반응 · 심플렉틱 적분기

7. Footnotes[편집]

1. 이 소거 문제는 지수 적분기 구현자가 반드시 한 번은 밟는 지뢰다. “이론상 4차인데 실측이 2차가 나온다”는 버그 리포트의 절대다수가 여기서 나온다. 알고리즘이 틀린 게 아니라 $\varphi_1$ 이 틀린 것이다. ↩

2. 지수 적분기의 아이디어 자체는 1960년대(Certaine, Lawson)까지 거슬러 올라가지만 30년 넘게 잊혀 있었다. 이유는 단순하다 — 그때는 $e^{hL}v$ 를 실용적으로 계산할 방법이 없었다. 사드의 크릴로프 근사(1992)와 FFT 하드웨어가 갖춰지고 나서야 되살아났으니, 알고리즘도 인프라가 준비돼야 데뷔한다. ↩

3. 여담으로, “지수 적분기”라는 이름 때문에 “그럼 지수적으로 빠른가요”라는 질문을 받는다고 한다. 아니다. 지수함수를 쓴다는 뜻이다. 굳이 따지면 수렴은 여전히 다항식 차수다. ↩