암시적 룽게-쿠타법

수치해석 마지막 수정: 2026-08-12 04:38:09

1. 개요[편집]

암시적 룽게-쿠타법 Implicit Runge–Kutta method
약칭	IRK
특징	부처 표의 계수행렬 A가 꽉 차 있다 (상삼각 성분 ≠ 0)
차수 상한	s단으로 최대 2s차 — 가우스-르장드르가 유일하게 달성
대표 계열	가우스-르장드르 · Radau IIA · Lobatto IIIC · SDIRK/ESDIRK
대가	매 스텝 sn차원 비선형 연립계 → 나이브하게 풀면 (sn)³
강성에서의 함정	차수 감소(order reduction) — 관측 차수가 단 차수로 떨어진다

명시적 RK는 계단을 한 칸씩 밟는다. 암시적 RK는 계단 전체를 한꺼번에 푼다. 밟는 값이 서로를 참조하기 때문이다.

암시적 룽게-쿠타법(implicit Runge–Kutta method, IRK)은 룽게-쿠타법의 부처 표(Butcher tableau)에서 계수행렬 $A = (a_{ij})$ 의 상삼각 부분이 0이 아닌 도식들, 즉 각 스테이지 값이 자기 자신과 뒤쪽 스테이지까지 참조하는 도식들을 말한다. 명시적 도식이 $a_{ij} = 0\ (j \ge i)$ 라는 제약 아래에서 계수를 고르는 것과 달리, 이쪽은 $s^{2}+s$ 개 계수를 전부 자유롭게 쓴다. 그 자유가 사 오는 것이 차수와 안정성이고, 지불하는 것이 연립 비선형계 풀이다.

$s$ 단 IRK는 다음 스테이지 방정식을 푼다.

Y_i = y_n + \Delta t \sum_{j=1}^{s} a_{ij}\, f(t_n + c_j\Delta t,\ Y_j), \qquad i = 1,\dots,s

y_{n+1} = y_n + \Delta t \sum_{i=1}^{s} b_i\, f(t_n + c_i\Delta t,\ Y_i)

$A$ 가 하삼각이 아니면 $Y_1,\dots,Y_s$ 가 서로를 참조하므로 순차적으로 계산할 수 없다. $y \in \mathbb{R}^{n}$ 일 때 미지수가 $sn$ 개인 하나의 비선형 방정식계를 매 스텝 푸는 것이 IRK의 본질이며, 이 문서 나머지의 절반은 그 비용을 어떻게 깎느냐에 대한 이야기다.

존재 이유는 명확하다. 달퀴스트 제2장벽 때문에 A-안정인 선형다단계법은 2차를 넘을 수 없다. 강성 문제에서 고차 정확도와 무조건 안정성을 동시에 원한다면 갈 곳은 여기뿐이다. 안정영역과 A/L-안정성의 정의, 그리고 그 그림 자체는 강성 방정식 문서에 이미 있으니 여기서는 계열별 결론만 정리한다.

2. 배점법으로 보기 — 계수가 하늘에서 오지 않는다[편집]

IRK 계수는 대부분 배점법(collocation)에서 유도된다. 구간 $[t_n, t_n+\Delta t]$ 에서 $s$ 차 다항식 $P(t)$ 를 찾되, $P(t_n) = y_n$ 이고 서로 다른 배점 노드 $c_1,\dots,c_s$ 에서 미분방정식을 정확히 만족하도록 요구한다.

P'(t_n + c_i\Delta t) = f\big(t_n + c_i\Delta t,\ P(t_n + c_i\Delta t)\big)

이러면 계수가 자동으로 결정된다. $\ell_j$ 를 노드에 대한 라그랑주 기저라 할 때

a_{ij} = \int_{0}^{c_i}\ell_j(\tau)\,d\tau, \qquad b_j = \int_{0}^{1}\ell_j(\tau)\,d\tau

즉 $b$ 는 그 노드에 대응하는 수치적분 구적 가중치이고, IRK의 차수는 그 구적법의 차수를 물려받는다. 노드를 어디에 찍느냐가 계열의 이름이 된다.

노드 선택	계열	단수 대비 차수	노드 위치
가우스-르장드르 점	가우스	$2s$	내부에만
Radau 점 ( $c_s = 1$ )	Radau IIA	$2s-1$	오른쪽 끝 포함
Lobatto 점 ( $c_1=0,\ c_s=1$ )	Lobatto IIIA	$2s-2$	양 끝 포함

구적법에서 가우스 점이 $2s$ 차, Radau 점이 $2s-1$ 차, Lobatto 점이 $2s-2$ 차 정확도를 주는 그 순서 그대로다. 그리고 $s$ 단 IRK의 차수 상한이 $2s$ 라는 것이 부처의 정리이며, 그것을 달성하는 도식은 가우스-르장드르 하나뿐이다.¹

3. 세 계열[편집]

3.1. 가우스-르장드르[편집]

노드는 이동 르장드르 다항식 $P_s(2x-1)$ 의 근이다(직교다항식). $s$ 단 $2s$ 차. 1단이 그 유명한 암시적 중점법(implicit midpoint), 2단 4차는 이렇게 생겼다.

\begin{array}{c|cc} \tfrac12-\tfrac{\sqrt3}{6} & \tfrac14 & \tfrac14-\tfrac{\sqrt3}{6}\\[2pt] \tfrac12+\tfrac{\sqrt3}{6} & \tfrac14+\tfrac{\sqrt3}{6} & \tfrac14\\[2pt]\hline & \tfrac12 & \tfrac12 \end{array}

성질이 화려하다. A-안정이고, 대칭(시간 역전 가능)이며, 심플렉틱이다. 심플렉틱 조건 $b_i a_{ij} + b_j a_{ji} = b_i b_j$ 를 정확히 만족하기 때문이고, 그래서 해밀턴계 장시간 적분에서 심플렉틱 적분기로 쓰인다. 명시적 RK4가 에너지를 표류시키는 문제를 근본적으로 피할 수 있는 고차 선택지다.

대신 L-안정이 아니다. 안정함수가 $e^{z}$ 의 대각 파데 근사 $R_{ss}$ 라서 $R(\infty) = (-1)^{s}$ , 즉 절댓값 1이다. 무한히 강성인 모드가 감쇠하지 않고 매 스텝 부호만 뒤집으며 남아 있다 — 크랭크-니콜슨법의 사다리꼴 링잉과 정확히 같은 병이고, 실제로 사다리꼴은 Lobatto IIIA 2단이다. 화학반응처럼 빠른 모드를 확실히 죽여야 하는 문제에는 부적합하다.

3.2. Radau IIA[편집]

노드는 오른쪽 끝점을 강제로 포함한 Radau 구적점. $s$ 단 $2s-1$ 차, L-안정, 그리고 강성 정확(stiffly accurate)하다. 강성 정확이란 $a_{sj} = b_j$ 가 성립해 $y_{n+1} = Y_s$ 가 되는 성질로, 마지막 스테이지 값이 곧 다음 스텝 값이라는 뜻이다. 2단 3차는 다음과 같다.

\begin{array}{c|cc} \tfrac13 & \tfrac{5}{12} & -\tfrac{1}{12}\\[2pt] 1 & \tfrac34 & \tfrac14\\[2pt]\hline & \tfrac34 & \tfrac14 \end{array}

3단 5차 판본이 하이러-반너의 RADAU5 로 구현되어 강성 ODE와 미분대수방정식 솔버의 사실상 표준이 되었다(SciPy의 Radau, Julia의 RadauIIA5 도 같은 계열). 강성 정확성 덕에 지표 1~3의 DAE에서 대수 구속조건이 스텝 끝에서 정확히 만족된다는 점이 결정적이다. 안정함수는 부대각 파데 근사 $R_{s-1,s}$ 이고 $R(\infty)=0$ .

3.3. Lobatto IIIC[편집]

양 끝점을 모두 포함하는 Lobatto 노드. $s$ 단 $2s-2$ 차로 차수는 손해지만 L-안정이고 대수적으로 안정하다. 2단 2차:

\begin{array}{c|cc} 0 & \tfrac12 & -\tfrac12\\[2pt] 1 & \tfrac12 & \tfrac12\\[2pt]\hline & \tfrac12 & \tfrac12 \end{array}

첫 행에 음수 계수가 들어 있다는 것이 눈에 띄는데, 이 때문에 Lobatto IIIC는 배점법이 아니다(같은 노드의 배점법은 Lobatto IIIA, 즉 사다리꼴 계열이다). 노드에서 방정식을 만족시키는 대신 안정성을 사 온 “불연속 배점법”이며, IIIA/IIIB가 못 가진 강한 감쇠를 얻는다. 지표 2 DAE와 강한 감쇠가 필요한 문제에서 쓰인다.

대수적 안정성(algebraic stability)은 $b_i \ge 0$ 이고 $M_{ij} = b_i a_{ij} + b_j a_{ji} - b_i b_j$ 가 준정부호라는 조건으로, 이게 성립하면 단조 감쇠하는 비선형 문제에서 수치해도 반드시 감쇠한다(B-안정성). 가우스·Radau IA·Radau IIA·Lobatto IIIC가 대수적으로 안정하고, Lobatto IIIA/IIIB는 아니다. 선형 A-안정성만으로는 비선형 문제의 안전을 보장하지 못한다는 점이 여기서 갈린다.

4. 비용 — (sn)³ 를 어떻게 깎나[편집]

스테이지 방정식에 뉴턴-랩슨법을 적용하면 매 반복마다 다음을 풀어야 한다.

\big(I_{sn} - \Delta t\,(A \otimes J)\big)\,\Delta = -F, \qquad J = \frac{\partial f}{\partial y}

크기가 $sn \times sn$ 이므로 LU 분해 비용이 $(sn)^{3}/3$ , 즉 같은 크기 후진 오일러의 $s^{3}$ 배다. $s=3$ 이면 27배. 이대로면 아무도 안 쓴다.

실전 코드가 쓰는 핵심 기법은 $A^{-1}$ 의 대각화다. 양변에 $(A^{-1}\otimes I)$ 를 곱하면 계수행렬이

A^{-1}\otimes I - \Delta t\,(I \otimes J)

가 되고, $A^{-1} = T\Lambda T^{-1}$ 로 대각화하면 $s$ 개의 독립적인 $n\times n$ 계 $(\lambda_i I - \Delta t J)$ 로 분리된다. Radau IIA 3단의 경우 $A^{-1}$ 의 고유값이 실수 하나( $\gamma \approx 3.6378$ )와 켤레 복소쌍 하나( $\approx 2.6811 \pm 3.0504i$ )라, 실수 $n\times n$ LU 하나 + 복소 $n\times n$ LU 하나로 끝난다. 복소 LU가 실수의 약 4배 비용이니 총 $\approx 5n^{3}/3$ — 나이브한 $9n^{3}$ 대비 5배 이상 절약이다. RADAU5가 실용적인 이유가 이 한 줄이다.²

그래도 비싸다. 그래서 나온 실무적 타협이 $A$ 를 삼각화하는 것이다.

DIRK(diagonally implicit RK). $a_{ij}=0\ (j>i)$ , 즉 $A$ 를 하삼각으로. 스테이지를 $Y_1, Y_2, \dots$ 순서대로 하나씩 풀 수 있어 매번 $n\times n$ 계만 다룬다. 대가는 차수 — $s$ 단 DIRK의 차수는 $s+1$ 을 넘지 못한다. 가우스가 $2s$ 를 내는 것과 비교하면 자유도를 버린 값을 그대로 치른 셈이다.
SDIRK(singly diagonally implicit). 여기에 $a_{ii} = \gamma$ 로 대각을 전부 같게 만든다. 그러면 $(I - \gamma\Delta t J)$ 하나만 LU 분해해서 모든 스테이지가 재사용한다. $\Delta t$ 를 안 바꾸는 한 스텝 사이에서도 재사용 가능.
ESDIRK(explicit first stage SDIRK). 첫 스테이지를 명시적으로( $a_{11}=0,\ c_1=0$ ) 두고 나머지 대각을 $\gamma$ 로 맞춘다. 첫 스테이지가 공짜인 데다 단 차수를 2로 올릴 수 있어 아래에서 볼 차수 감소에 훨씬 강하다. 여기에 마지막 행을 $b$ 와 같게 잡아 강성 정확성까지 챙긴다. 전산유체역학과 연소 시뮬레이션의 IMEX 도식(케네디-카펜터의 ARK 계열)이 대부분 ESDIRK를 암시적 절반으로 쓴다.

한 걸음 더 나가면 뉴턴 반복 자체를 포기하고 선형 방정식 한 번으로 스테이지를 끝내는 로젠브록 방법이 된다. 수렴 실패가 없다는 장점과 자코비안 정확도에 성능이 매달린다는 단점을 함께 얻는다.

5. 차수 감소 — 논문의 차수와 내 코드의 차수가 다른 이유[편집]

강성 문제에서 5차 방법을 돌렸는데 수렴 그래프 기울기가 3이 나오면, 대개 버그가 아니라 차수 감소(order reduction)다.

원인은 고전적 차수 이론이 $\Delta t \to 0$ 극한의 테일러 전개에 기반한다는 데 있다. 강성 문제에서 실제로 쓰는 것은 $\Delta t\lVert J\rVert \gg 1$ 인 영역이고, 거기서는 그 전개가 무의미하다. 이 영역의 거동을 지배하는 것은 고전적 차수 $p$ 가 아니라 단 차수(stage order) $q$ — 각 스테이지 값 $Y_i$ 자체가 참해 $y(t_n + c_i\Delta t)$ 를 몇 차로 근사하느냐다.

\sum_{j} a_{ij}c_j^{k-1} = \frac{c_i^{k}}{k},\qquad k = 1,\dots,q

계열	고전 차수 $p$	단 차수 $q$
가우스	$2s$	$s$
Radau IIA	$2s-1$	$s$
Lobatto IIIC	$2s-2$	$s-1$
SDIRK	$\le s+1$	1
ESDIRK	$\le s+1$	2

강성 극한에서 관측되는 차수는 대략 $q$ 또는 $q+1$ 수준으로 떨어진다. DIRK 계열의 단 차수가 1이라는 사실이 이 계열의 진짜 약점이고, 3차·4차 SDIRK가 강성 문제에서 1~2차처럼 행동하는 일이 흔하다. 반면 Radau IIA는 단 차수가 $s$ 라 감소 폭이 훨씬 작다 — 비싼 꽉 찬 $A$ 를 감수하는 값이 여기서 회수된다. 시간 의존 경계조건이 있는 PDE의 시간 적분, 그리고 지표 2 이상 DAE에서 이 현상이 특히 두드러진다.³

엄밀한 언어로는 B-수렴(B-convergence) 이론이 이걸 다룬다. 강성 매개변수에 의존하지 않는 오차 상수를 요구하는 수렴 개념이고, 프로테로-로빈슨 문제 $y' = \lambda(y - g(t)) + g'(t)$ 가 표준 시험대다. 해가 $g(t)$ 로 뻔한데도 방법마다 관측 차수가 갈리는 것을 눈으로 볼 수 있다.

6. 스텝 제어 — 오차 추정이 까다롭다[편집]

명시적 RK에서 적응 스텝은 쉬웠다. 차수가 다른 두 해를 같은 스테이지로 뽑아 차이를 재면 그만이었다(내장 쌍, embedded pair). IRK에서는 사정이 나쁘다. 스테이지가 몇 개 안 되고 계수 자유도를 이미 차수에 다 써 버렸기 때문에, 같은 스테이지 위에 얹을 여분의 낮은 차수 해가 잘 안 나온다.

RADAU5가 쓰는 방법은 이렇다. 스테이지 값들의 선형결합에 초기 도함수 $f(t_n,y_n)$ 을 하나 더 섞어 저차 근사를 만들고, 그 차이를 오차 추정으로 삼되 추정값에 $(I - \gamma\Delta t J)^{-1}$ 를 한 번 더 적용한다. 이미 분해해 둔 실수 LU를 재사용하므로 비용이 거의 공짜이고, 이 필터를 거치지 않으면 강성 극한에서 추정치가 터무니없이 커져 스텝이 무한히 줄어드는 사고가 난다. 강성 문제에서는 오차 추정기 자체도 안정해야 한다는 교훈이 여기 박혀 있다.

한 걸음 더 나가면 가변 차수다. Radau IIA는 노드 수만 바꾸면 같은 구조로 5차·9차·13차를 만들 수 있어서, 요구 허용오차가 빡빡하면 단수를 올리고 느슨하면 내리는 전략이 가능하다. 하이러-반너의 RADAU 코드가 정확히 그렇게 동작한다. 스텝 크기와 차수를 함께 조절하는 이 구조는 BDF 코드의 가변 차수 전략과 발상이 같다.

7. 언제 무엇을 쓰나[편집]

화학 반응·연소 시뮬레이션·반응속도론 처럼 강성비가 $10^{9}$ 급이고 $n$ 이 수십~수백: Radau IIA(RADAU5). 비용이 $n^{3}$ 이라도 $n$ 이 작으니 감당되고, L-안정 + 강성 정확 + 낮은 차수 감소를 전부 가져간다.
대형 PDE 반이산화(격자점 수십만~수백만, $J$ 는 희소행렬): ESDIRK + IMEX. 스테이지를 하나씩 풀어야 전처리기와 크리로프 부분공간법을 붙일 수 있고, 꽉 찬 $A$ 의 $sn$ 차원 계는 애초에 크리로프 친화적이지 않다.
해밀턴계 장시간 적분: 가우스-르장드르. 고차 심플렉틱이 필요하고 강성이 아니라면 자연스러운 선택. 다만 스텝마다 비선형계를 풀어야 하므로, 싸고 명시적인 베를레 적분을 이길 만한 이유가 있어야 한다.
DAE: 강성 정확한 방법(Radau IIA, ESDIRK, Lobatto IIIC). 강성 정확하지 않은 방법은 대수 구속조건이 스텝 끝에서 어긋난다.
적당한 강성 + 적당한 정확도: 사실 선형다단계법 계열의 BDF(ode15s, CVODE)가 여전히 가성비 1등인 경우가 많다. 스텝당 자코비안 하나, LU 하나면 끝나기 때문이다. IRK는 BDF가 안 되는 곳(고차 정확도, 진동 모드, 심플렉틱, 심한 차수 감소)에서 값을 한다.

8. 관련 문서[편집]

명시적 쪽에는 “4차를 넘으면 스테이지가 차수보다 많아진다”는 부처 장벽이 있는데, 암시적 쪽은 정반대로 스테이지 하나가 차수 둘을 산다. 자유 매개변수가 $s^2+s$ 개로 늘어난 대가를 정직하게 돌려받는 셈이다. 물론 그 대가는 매 스텝 $sn$ 차원 뉴턴이라는 청구서로 따로 온다. ↩
$A^{-1}$ 이 항상 예쁘게 대각화되는 것은 아니다. 고유값이 겹치면 조르당 블록이 생겨 이 기법이 안 먹히고, 대각화 변환 $T$ 의 조건수가 나쁘면 정확도가 깎인다. Radau IIA 3단은 다행히 둘 다 문제가 없어서 하이러-반너가 상수를 아예 소스코드에 박아 놓았다. 수치해석 코드에서 열몇 자리 상수가 하드코딩되어 있으면 대개 이런 사연이 있다. ↩
“논문에는 5차라던데 왜 3차가 나오죠?” 는 강성 솔버 사용자 게시판의 스테디셀러 질문이다. 답은 대개 “당신 문제가 강성이라서”인데, 정작 강성이라서 이 방법을 골랐다는 게 아이러니다. 고전적 차수는 비강성 극한의 성질이고, 강성 솔버는 정의상 그 극한에서 쓰이지 않는다. ↩