암시적 룽게-쿠타법

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-안정성만으로는 비선형 문제의 안전을 보장하지 못한다는 점이 여기서 갈린다.

3.4. 검산: 위 세 표와 성질들을 배점법에서 다시 뽑기

노드만 주고 $a_{ij} = \int_0^{c_i}\ell_j$, $b_j = \int_0^1 \ell_j$ 를 유리수 산술로 계산해, 본문에 적은 표와 성질을 전부 재확인했다.

방법	차수 $p$	단 차수 $q$	$R(\infty)$	강성 정확	대수적 안정
가우스 1단 (암시적 중점)	2	1	$-1$	아니다	그렇다
Radau IIA 1단 (후진 오일러)	1	1	$0$	그렇다	그렇다
Radau IIA 2단	3	2	$0$	그렇다	그렇다
Lobatto IIIA 2단 (사다리꼴)	2	2	$-1$	그렇다	아니다
Lobatto IIIC 2단	2	1	$0$	그렇다	그렇다

• Radau IIA 2단의 계수가 유리수로 $A = \begin{pmatrix} 5/12 & -1/12 \\ 3/4 & 1/4\end{pmatrix}$, $b = (3/4,\ 1/4)$ 로 나왔다 — 위 표와 같다. $a_{2j} = b_j$ 가 성립해 강성 정확이고, $p = 3 = 2s-1$, $q = 2 = s$ 다.
• 사다리꼴이 정말 Lobatto IIIA 2단이다. $c=(0,1)$ 의 배점법이 $A = \begin{pmatrix} 0 & 0 \\ 1/2 & 1/2\end{pmatrix}$ 를 주고 $R(\infty) = -1$ 이 나와 링잉 주장과 맞는다. 그리고 대수적으로 안정하지 않다는 것도 확인됐다.
• Lobatto IIIC 2단은 본문의 표(첫 행에 음수)를 그대로 넣었더니 $p=2 = 2s-2$, $q = 1 = s-1$, $R(\infty)=0$ 이고, $M = \begin{pmatrix} 1/4 & -1/4 \\ -1/4 & 1/4\end{pmatrix}$ 가 고유값 $0,\ 1/2$ 의 준정부호라 대수적으로 안정하다.
• 가우스 2단은 계수가 $1/4$, $1/4 \mp \sqrt{3}/6$, $b=(1/2,1/2)$ 로 본문 표와 일치했고, 심플렉틱 조건 $b_i a_{ij} + b_j a_{ji} - b_i b_j$ 의 최댓값이 $1.1\times10^{-16}$ — 반올림 한계다. 안정함수는 $e^z$ 의 대각 파데 근사 $\frac{1 + z/2 + z^2/12}{1 - z/2 + z^2/12}$ 와 시험한 모든 $z$ 에서 14자리까지 같고, $R(\infty) = +1 = (-1)^s$ 라 A-안정이지만 L-안정이 아니다는 서술이 확인된다.

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가 실용적인 이유가 이 한 줄이다.²

4.1. 검산: Radau IIA 3단의 $A^{-1}$ 고유값

위 세 숫자가 이 절의 전부이므로 직접 뽑았다. 노드 $c = \left(\frac{4-\sqrt6}{10},\ \frac{4+\sqrt6}{10},\ 1\right)$ 로 배점법 계수를 만들고, $A^{-1}$ 의 특성다항식을 카르다노 공식으로 풀었다.

고유값	값	절댓값
실수 $\gamma$	$3.6378342527$	$3.6378342527$
켤레쌍	$2.6810828736 \pm 3.0504301992\,i$	$4.0611980715$

본문의 $\gamma \approx 3.6378$, $\approx 2.6811 \pm 3.0504i$ 와 자리까지 맞는다. 덤으로 $A^{-1}$ 의 대각합이 정확히 9, 행렬식이 정확히 60 이었다 — 세 고유값이 무리수인데 대칭함수는 정수로 떨어진다. 구적조건 $B(k)$ 는 $k \le 5$ 에서 잔차 $10^{-16}$, $k=6$ 에서 $1.7\times10^{-3}$ 이라 $p = 5 = 2s-1$ 이고, 단조건 $C(k)$ 는 $k \le 3$ 까지만 성립해 $q = 3 = s$ 다.

비용 산술도 맞는다. 실수 LU $n^3/3$ + 복소 LU $4n^3/3 = 5n^3/3$ 이고, 나이브한 $(3n)^3/3 = 9n^3$ 과의 비는 $27/5 = 5.4$ 배다.

그래도 비싸다. 그래서 나온 실무적 타협이 $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)$ 로 뻔한데도 방법마다 관측 차수가 갈리는 것을 눈으로 볼 수 있다.

5.1. 검산: 차수 감소를 프로테로-로빈슨으로 직접 재기

“눈으로 볼 수 있다”고 적어 놓고 안 보여 주면 안 되니 실제로 쟀다. $g(t) = \cos t$, $y(0)=g(0)$, $T=1$ 이고 스테이지 방정식이 선형이라 뉴턴 없이 정확히 풀었다(뉴턴 수렴 오차가 섞이지 않는다). 아래는 $N = 40 \to 320$ 구간의 관측 차수다.

방법	$p$	$q$	비강성 $\lambda=-1$	강성 $\lambda=-10^{4}$
가우스 2단	4	2	4.00	2.59 → 3.98
Radau IIA 2단	3	2	3.00	2.01
Radau IIA 3단	5	3	4.99	3.02
SDIRK 2단 (L-안정)	2	1	2.00	1.08
SDIRK 3차 (Crouzeix)	3	1	3.00	2.00
ESDIRK (TR-BDF2)	2	2	2.00	2.00
Lobatto IIIC 2단	2	1	2.00	1.02

비강성 열에서는 모두 고전 차수가 그대로 나오고, 강성 열에서는 전부 $q$ 또는 $q+1$ 로 내려앉는다. 본문의 ”$q$ 또는 $q+1$ 수준”이 실측으로 확인되는 셈이다. 읽을 거리가 셋 있다.

• 각주의 그 질문이 표에 있다. Radau IIA 3단은 고전 차수 5차인데 강성에서 3.02 다. “논문에는 5차라던데 왜 3차가 나오죠”의 답이 $q=3$ 이다.
• DIRK 의 약점이 그대로 보인다. SDIRK 2단은 2차 → 1.08, Crouzeix 3차는 3차 → 2.00 으로, 단 차수 1 인 계열이 강성에서 1~2차처럼 행동한다는 서술과 맞는다.
• ESDIRK 가 왜 IMEX 의 암시적 절반인지도 보인다. 단 차수가 2 라 강성에서도 2.00 을 그대로 유지한다 — 표에서 유일하게 차수를 안 잃은 방법이다.

$\lambda = -10^{6}$ 까지 내리면 여러 방법의 오차가 $10^{-14}$ 대의 반올림 바닥에 닿아 관측 차수가 의미를 잃는다. 강성 극한의 차수를 잴 때는 오차가 바닥에 붙지 않는 구간을 골라야 한다는 것도 이 실험의 부산물이다.

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등인 경우가 많다(강성 방정식 이 달퀴스트 장벽과 BDF1~6 의 A($\alpha$) 각을 본문에서 다룬다). 스텝당 자코비안 하나, LU 하나면 끝나기 때문이다. IRK는 BDF가 안 되는 곳(고차 정확도, 진동 모드, 심플렉틱, 심한 차수 감소)에서 값을 한다.

8. 관련 문서

• 룽게-쿠타법 · 강성 방정식 · 크랭크-니콜슨법
• 심플렉틱 적분기 · 베를레 적분
• 뉴턴-랩슨법 · 자코비안 행렬 · LU 분해 · 희소행렬
• 미분대수방정식 · 지수 적분기 · 연산자 분리
• 파데 근사 · 수치적분 · 직교다항식
• 연소 시뮬레이션 · 반응속도론 · 전처리기

9. Footnotes

1. 명시적 쪽에는 “4차를 넘으면 스테이지가 차수보다 많아진다”는 부처 장벽이 있는데, 암시적 쪽은 정반대로 스테이지 하나가 차수 둘을 산다. 자유 매개변수가 $s^2+s$ 개로 늘어난 대가를 정직하게 돌려받는 셈이다. 물론 그 대가는 매 스텝 $sn$ 차원 뉴턴이라는 청구서로 따로 온다. ↩

2. $A^{-1}$ 이 항상 예쁘게 대각화되는 것은 아니다. 고유값이 겹치면 조르당 블록이 생겨 이 기법이 안 먹히고, 대각화 변환 $T$ 의 조건수가 나쁘면 정확도가 깎인다. Radau IIA 3단은 다행히 둘 다 문제가 없어서 하이러-반너가 상수를 아예 소스코드에 박아 놓았다. 수치해석 코드에서 열몇 자리 상수가 하드코딩되어 있으면 대개 이런 사연이 있다. ↩

3. “논문에는 5차라던데 왜 3차가 나오죠?” 는 강성 솔버 사용자 게시판의 스테디셀러 질문이다. 답은 대개 “당신 문제가 강성이라서”인데, 정작 강성이라서 이 방법을 골랐다는 게 아이러니다. 고전적 차수는 비강성 극한의 성질이고, 강성 솔버는 정의상 그 극한에서 쓰이지 않는다. ↩