수치해석

1부의 앞 열한 장은 전부 닫힌 해를 찾는 이야기였다. 적분을 손으로 하고, 급수를 전개하고, 변수를 분리해 고유함수를 얻었다. 그런데 반도체 소자 하나를 앞에 놓으면 그 방법이 거의 다 막힌다. 도핑이 위치마다 다르고, 이동도가 전기장에 따라 변하고, 캐리어 농도가 퍼텐셜의 지수함수라서, 푸아송 방정식 하나조차 닫힌 해가 없다. 그 자리에서 쓰는 것이 이 장이다.

수치해석이 하는 일은 딱 하나다 — 연속을 이산으로 바꾸고, 그때 무엇을 잃는지 잰다. 함수를 점 몇 개로, 미분을 차분으로, 적분을 합으로, 무한 반복을 유한 반복으로 바꾼다. 바꾸는 것 자체는 쉽다. 어려운 것은 바꾼 뒤에 남는 두 질문이다.

• 얼마나 틀리는가. 이산화가 만든 오차(절단 오차)와 컴퓨터가 만든 오차(반올림 오차)가 따로 있고, 둘은 반대 방향으로 움직인다. 한쪽을 줄이면 다른 쪽이 커진다.
• 믿어도 되는가. 답이 나왔다고 해서 맞는 답이 아니다. 반복이 엉뚱한 곳으로 수렴하거나, 시간 전진이 폭발하거나, 조건이 나쁜 문제라 입력의 미세한 흔들림이 답을 통째로 바꿀 수 있다.

이 장의 열네 절은 모두 이 두 질문의 변주다. 1절에서 오차의 종류를 나눈 뒤, 2절부터 9절까지는 정적인 문제 — 방정식 풀기, 함수 근사, 적분, 선형계, 고유값 — 를, 10절부터 12절까지는 시간이 흐르는 문제 — ODE와 PDE — 를 다룬다. 13절은 소자·전자기 시뮬레이터가 실제로 쓰는 두 이산화 방식(유한요소법·모멘트법)의 뼈대이고, 14절은 “가장 좋은 파라미터를 찾는” 최적화다.

이 장이 반도체 어디로 이어지는지 미리 적어 둔다. 소자 시뮬레이터(TCAD)는 푸아송·연속 방정식을 격자에서 이산화하고(11·12절) 뉴턴법으로 비선형계를 풀며(2절) 그때 나오는 희소 선형계를 LU나 반복법으로 처리한다(5·6절). 회로 시뮬레이터(SPICE)는 강성 ODE를 암시적으로 적분하고(10절), 소자 모델의 파라미터는 최소제곱으로 추출하며(8절·14절), 밴드 구조는 큰 에르미트 행렬의 고유값을 반복법으로 구한다(9절). 공정 시뮬레이션의 확산은 포물형 PDE(12절), 배선 기생 커패시턴스는 모멘트법(13절), 기계적 응력은 유한요소법(13절)이다. 뒤의 장들이 “수치적으로 푼다”고 한 줄로 넘기는 자리마다 이 장의 절 하나가 들어 있다.

기호 — 이 장의 약속

수치해석은 다른 장보다 기호가 재활용된다. 같은 글자가 절마다 다른 뜻으로 쓰이는 자리가 많아서, 처음에 한 번 정리해 두고 절이 바뀔 때마다 다시 못박는다.

기호	읽는 법	무엇을 가리키는가	절
$a,\ \tilde a$	참값 · 근삿값	물결표(틸데)가 붙으면 계산된 값	전체
$\varepsilon$	오차	$\varepsilon = a - \tilde a$ (참값 − 근삿값)	1
$\varepsilon_r$	상대 오차	$\varepsilon / a$	1
$u$	반올림 단위	부동소수점 한 번의 상대 오차 상한. 배정도에서 $1.1 \times 10^{-16}$	1
$h$	걸음 폭	독립 변수(대개 $x$)의 격자 간격	3–12
$k$		시간 걸음 폭 $\Delta t$ (11·12절) · 스플라인의 기울기 $k_j$ (3절) · 룽게–쿠타의 기울기 $k_1 \ldots k_4$ (10절)	3·10–12
$x_n,\ y_n$	아래첨자 $n$	스칼라 반복·걸음의 $n$번째 값	2·10
$\mathbf{x}^{(m)}$	위첨자 $(m)$	벡터 반복의 $m$번째 값. 성분 인덱스와 겹치지 않게 위로 올린다	6·9
$u_{ij}$		격자점 $(x_i,\ y_j)$ 또는 $(x_i,\ t_j)$에서의 근사값	11·12
$p_n(x)$		$n$차 보간 다항식	3
$L_k(x)$	라그랑주 기저	$k$번째 마디에서 1, 나머지에서 0인 다항식	3
$f[x_0, x_1]$	분할 차분	$\bigl(f(x_1) - f(x_0)\bigr)/(x_1 - x_0)$	3
$\xi$	크시	오차 공식에 나오는 위치를 모르는 점. 구간 안 어딘가	3·4
$L,\ U$		하삼각 · 상삼각 행렬. $A = LU$	5
$\lVert \cdot \rVert$	노름	벡터·행렬의 크기. 아래첨자 $1,\ 2,\ \infty$가 종류	7
$\kappa(A)$	조건수	$\lVert A \rVert\, \lVert A^{-1} \rVert$	7
$\mathbf{r}$	잔차	$\mathbf{b} - A\tilde{\mathbf{x}}$. 굵은 글자	7
$r,\ r^{*}$	격자비	포물형 $c^2 k / h^2$ · 쌍곡형 $c^2 k^2 / h^2$. 굵지 않은 스칼라	12
$\rho(C)$	스펙트럼 반지름	반복 행렬 $C$의 고유값 절댓값 최대. 수렴의 척도	6·9
$\omega$	오메가	과완화 계수(SOR). 이 장에서는 각진동수가 아니다	6
$\lambda$	람다	고유값 (1.3과 같다)	9·10
$q$		레일리 몫 $\mathbf{x}^{\mathsf{T}}A\mathbf{x} / \mathbf{x}^{\mathsf{T}}\mathbf{x}$. 이 장에서는 전하가 아니다	9
$J$	야코비 행렬	$\partial F_i / \partial x_j$. 벡터 뉴턴법의 도함수	2·8·14
$\nabla f$	기울기	1.1과 같다 — 후진 차분 기호로는 쓰지 않는다	14

주의

짚고 갈 함정이 넷이다.

① $\varepsilon$의 부호. 이 백과는 Kreyszig를 따라 $\varepsilon = a - \tilde a$, 즉 참값에서 근삿값을 뺀다. 그래서 $a = \tilde a + \varepsilon$이다 — “참값은 근삿값에 오차를 더한 것”. 반대로 $\tilde a - a$로 정의하는 책도 많다. 부호만 다르므로 크기 논의에는 영향이 없지만, 오차 전파 식의 부호를 따라갈 때는 어느 쪽인지 확인해야 한다.

② $k$가 세 가지다. 3절의 $k_j$는 스플라인의 마디 기울기, 10절의 $k_1 \ldots k_4$는 룽게–쿠타의 중간 기울기, 11·12절의 $k$는 시간 걸음 $\Delta t$다. 셋 다 Kreyszig의 표기이고, 서로 다른 절에서만 나온다. 볼츠만 상수와 파수도 이 장에는 없다.

③ $r$과 $\mathbf{r}$. 굵은 $\mathbf{r}$은 잔차 벡터(7절), 굵지 않은 $r$은 격자비(12절)다. Kreyszig는 격자비를 $r$로, 잔차를 $\mathbf{r}$로 쓰므로 같은 표기다. 3절에서 Newton 전진 공식의 매개변수 $r = (x - x_0)/h$도 잠깐 나오는데, 그 절 안에서만 쓴다.

④ $\nabla$는 기울기다. Kreyszig는 후진 차분을 $\nabla f_j = f_j - f_{j-1}$로 적는다. 1.1의 기울기 기호와 부딪치므로 이 백과는 후진 차분 기호를 쓰지 않는다. 필요하면 전진 차분 $\Delta$로 옮겨 적거나 성분을 풀어 쓴다.

오차와 부동소수점

세 가지 오차

수치 계산의 답이 참값과 다른 이유는 세 층으로 나뉜다.

1. 모델링 오차. 물리를 수식으로 옮길 때 이미 버린 것. 이동도를 상수로 놓았다면 그 순간부터 답은 실제 소자와 다르다. 수치해석의 대상이 아니다 — 이 층은 물리가 책임진다.
2. 절단 오차(truncation error). 무한을 유한으로 자를 때 생긴다. 테일러 급수를 몇 항에서 끊고, 적분을 유한 합으로 바꾸고, 미분을 차분으로 바꾸는 순간 생긴다. 이산화 오차라고도 한다. 걸음 폭 $h$를 줄이면 준다.
3. 반올림 오차(round-off error). 컴퓨터가 실수를 유한 자릿수로만 저장하기 때문에 생긴다. 연산 한 번마다 조금씩 쌓인다. $h$를 줄이면 연산 횟수가 늘어 오히려 커진다.

이 장의 방법들은 전부 2번을 줄이는 기술이고, 그 기술이 3번을 얼마나 키우는지가 각 방법의 한계다. 둘 사이의 줄다리기를 먼저 보고 시작한다.

부동소수점 — 숫자를 어떻게 저장하는가

컴퓨터는 실수를 가수 × 밑의 거듭제곱 꼴로 저장한다. Kreyszig의 표기로 십진 $t$자리 부동소수점 수는

$$\tilde a = \pm\, 0.d_1 d_2 \cdots d_t \times 10^{e}, \qquad d_1 \neq 0$$

이다. $d_1 \neq 0$이 정규화 조건이다 — 앞자리에 0을 낭비하지 않는다. 실제 컴퓨터는 밑이 2이고, 표준(IEEE 754) 배정도(double)에서는 가수가 53비트, 지수가 11비트다. 십진으로 옮기면 유효숫자 약 16자리, 크기 범위 약 $10^{-308}$에서 $10^{308}$이다.

직관

부동소수점 수는 수직선 위에 균일하게 놓여 있지 않다. 지수가 하나 오를 때마다 이웃한 두 수의 간격이 두 배로 벌어진다. 1 근처에서는 이웃까지의 거리가 $2^{-52} \approx 2.2 \times 10^{-16}$이고, $10^{6}$ 근처에서는 그 $10^{6}$배인 $2.2 \times 10^{-10}$이며, $10^{16}$을 넘으면 간격이 1보다 커져서 정수조차 전부 표현되지 않는다.

그래서 부동소수점의 정밀도는 절대량이 아니라 비율로 말해야 한다. “소수점 아래 16자리”가 아니라 “유효숫자 16자리”다. 이것이 아래의 모든 오차 논의가 상대 오차로 진행되는 이유다.

실수 $a$를 가장 가까운 부동소수점 수 $\text{fl}(a)$로 반올림하면, 상대 오차는 한 값을 넘지 않는다.

$$\text{fl}(a) = a\,(1 + \delta), \qquad \lvert \delta \rvert \le u$$

$u$를 반올림 단위(rounding unit, 단위 반올림 오차)라 한다. 십진 $t$자리에서는 $u = \tfrac12 \times 10^{1-t}$, 배정도에서는 $u = 2^{-53} \approx 1.1 \times 10^{-16}$이다. 1 다음의 부동소수점 수까지의 간격 $2^{-52}$를 머신 엡실론이라 부르는데, $u$의 두 배다. 책과 라이브러리에 따라 둘 중 아무거나 “엡실론”이라 부르므로 2배 차이는 신경 쓰지 않아도 되지만, 이 절에서는 Kreyszig를 따라 $u$를 쓴다.

사칙연산 하나하나도 같은 성질을 갖는다 — 정확히 계산한 뒤 반올림한 것과 같다.

$$\text{fl}(a \circ b) = (a \circ b)(1 + \delta), \qquad \lvert \delta \rvert \le u, \qquad \circ \in \{+,\ -,\ \times,\ \div\}$$

연산 한 번의 오차는 $10^{-16}$ 수준이다. 문제는 그것이 어떻게 쌓이는가다.

오차의 정의와 전파

참값 $a$, 근삿값 $\tilde a$에 대해

$$\varepsilon = a - \tilde a \quad(\text{오차}), \qquad \varepsilon_r = \frac{\varepsilon}{a} = \frac{a - \tilde a}{a} \quad(\text{상대 오차})$$

그리고 $\lvert \varepsilon \rvert \le \beta$인 $\beta$를 오차 한계라 한다. 참값을 모르니 오차 자체는 알 수 없고, 알 수 있는 것은 한계뿐이다.

두 근삿값을 더하거나 곱하면 오차가 어떻게 되는가. $a = \tilde a + \varepsilon_a$, $b = \tilde b + \varepsilon_b$를 넣어 보면

$$\begin{aligned} (a \pm b) - (\tilde a \pm \tilde b) &= \varepsilon_a \pm \varepsilon_b &&\Longrightarrow\quad \lvert \varepsilon_{a \pm b} \rvert \le \beta_a + \beta_b \\[4pt] \frac{ab - \tilde a \tilde b}{ab} &\approx \frac{\varepsilon_a}{a} + \frac{\varepsilon_b}{b} &&\Longrightarrow\quad \lvert \varepsilon_{r,\,ab} \rvert \lesssim \beta_{r,a} + \beta_{r,b} \end{aligned}$$

덧셈·뺄셈에서는 절대 오차 한계가 더해지고, 곱셈·나눗셈에서는 상대 오차 한계가 더해진다. 곱셈의 결과는 $\varepsilon_a \varepsilon_b$ 항을 버린 1차 근사다. 나눗셈도 같은 꼴이다.

이 두 줄이 뜻하는 바는 대칭적이지 않다. 곱셈은 상대 오차가 그저 더해지므로 얌전하다. 그러나 뺄셈은 절대 오차가 더해지는데 결과의 크기는 작아질 수 있으므로, 상대 오차가 폭발할 수 있다.

상쇄 — 뺄셈이 자릿수를 잡아먹는다

$a \approx b$인 두 수를 뺀다고 하자. 각각 유효숫자 16자리로 정확하더라도, 앞의 열 자리가 같다면 뺀 결과에는 여섯 자리만 남는다. 뒤의 열 자리가 0으로 채워져 정확해 보이지만 정보가 없는 자리다. 이것을 상쇄(cancellation)라 한다.

Kreyszig가 드는 예가 이차방정식이다. $x^2 - 40x + 2 = 0$의 두 근을 근의 공식으로 구하면

$$x = 20 \pm \sqrt{398} = 20 \pm 19.94994\ldots$$

작은 근 $x_2 = 20 - 19.94994 = 0.05006$은 유효숫자 4자리 계산에서 $20.00 - 19.95 = 0.05$로 한 자리만 남는다. 참값은 $0.050063\ldots$이다. 반면 큰 근 $x_1 = 39.95$는 멀쩡하다. 해결은 공식을 바꾸는 것이다 — 두 근의 곱이 상수항이므로

$$x_2 = \frac{c}{x_1} = \frac{2}{39.95} = 0.05006$$

뺄셈을 곱셈·나눗셈으로 바꾸면 상쇄가 사라진다. 같은 요령이 $\sqrt{x+1} - \sqrt{x}$를 $1/(\sqrt{x+1} + \sqrt{x})$로, $1 - \cos\theta$를 $2\sin^2(\theta/2)$로 바꾸는 데 쓰인다. 수학적으로 같은 식이 수치적으로는 전혀 다르다.

주의

반도체 계산에서 상쇄가 숨어 있는 자리.

① 순 재결합률. $R = \dfrac{np - n_i^2}{\tau(n + p + 2n_i)}$에서 평형 근처는 $np \approx n_i^2$이라 분자가 상쇄다. $np = 10^{20}$, $n_i^2 = 10^{20}$ 급의 두 수를 빼서 $10^{5}$ 급의 차이를 읽으려 하면 유효숫자가 열다섯 자리에서 다섯 자리로 준다.

② 전류 밀도. 드리프트 항과 확산 항은 평형 근처에서 크기가 같고 부호가 반대다. $J_n = q\mu_n n E + qD_n\, dn/dx$를 두 항으로 따로 계산해 더하면 상쇄가 난다. 12절에서 보는 Scharfetter–Gummel 이산화는 이 상쇄를 피하려고 두 항을 하나의 지수 표현으로 합쳐 버린 것이다.

③ 큰 수에 작은 수 더하기. 누설 전류 $10^{-15}\ \mathrm{A}$를 구동 전류 $10^{-3}\ \mathrm{A}$에 더하면 상대 크기가 $10^{-12}$라 $u = 10^{-16}$ 안에서 살아남긴 하지만, 여기에 몇 번의 상쇄가 겹치면 사라진다. SPICE가 전류 수렴 판정에 상대 허용치(RELTOL, 기본 $10^{-3}$)와 절대 허용치(ABSTOL, 기본 $10^{-12}\ \mathrm{A}$)를 따로 두는 이유가 이것이다 — 상대만 쓰면 0에 가까운 전류에서 판정이 불가능하고, 절대만 쓰면 큰 전류에서 판정이 무의미하다.

절단 오차와 반올림 오차의 줄다리기

가장 단순한 예로 둘의 관계를 본다. 도함수를 전진 차분으로 근사하면

$$f'(x) \approx \frac{f(x+h) - f(x)}{h}$$

테일러 전개 $f(x+h) = f(x) + hf'(x) + \tfrac12 h^2 f''(\xi)$를 넣으면 절단 오차는 $\tfrac12 h f''(\xi)$, 즉 $h$에 비례한다. 한편 분자의 두 함수값에는 각각 $\lvert f \rvert u$ 급의 반올림 오차가 있고, 그것이 $h$로 나뉘므로 반올림 오차는 $u/h$에 비례한다.

$$\lvert \varepsilon_{\text{total}} \rvert \lesssim \underbrace{\tfrac12 h \lvert f'' \rvert}_{\text{절단}} + \underbrace{\frac{2u\lvert f \rvert}{h}}_{\text{반올림}}$$

첫 항은 $h$를 줄이면 줄고, 둘째 항은 $h$를 줄이면 는다. 합이 최소가 되는 곳은 $h_{\text{opt}} \sim \sqrt{u} \approx 10^{-8}$이고, 그때 오차도 $\sqrt{u} \approx 10^{-8}$ 수준이다. 배정도로 계산해도 도함수는 여덟 자리밖에 못 얻는다.

그림 1.12-1. 전진 차분으로 f′(x)를 구할 때의 오차를 걸음 폭 h에 대해 양대수 축에 그린 것. 오른쪽에서 왼쪽으로 h를 줄이면 절단 오차(기울기 +1인 직선)는 줄지만 반올림 오차(기울기 −1인 직선)가 늘어, 총오차는 두 직선이 만나는 h ≈ √u ≈ 10⁻⁸ 근처에서 바닥을 찍고 다시 올라간다. 바닥의 높이가 곧 그 방법으로 얻을 수 있는 최선의 정밀도다. 중심 차분처럼 절단 오차의 차수가 높으면 왼쪽 직선이 가팔라져 바닥이 더 낮고 더 오른쪽으로 온다.직접 작도

직관

이 그림이 이 장 전체의 그림이다. 모든 수치 방법은 이 V자 곡선 위의 어딘가에 있다. 방법을 개선한다는 것은 대개 절단 오차의 차수를 올려 왼쪽 직선을 가파르게 만드는 일이고(4절의 중심 차분·심프슨, 10절의 룽게–쿠타), 그러면 바닥이 낮아지고 오른쪽으로 옮겨 가 더 큰 $h$로도 같은 정밀도를 얻는다. 반올림 직선을 낮추는 방법은 하나뿐이다 — 더 긴 가수(4배 정밀도)를 쓰거나, 상쇄를 피하도록 식을 고쳐 쓰는 것.

그래프의 축이 둘 다 대수라는 점도 중요하다. 절단 오차가 $h^p$에 비례하면 왼쪽 직선의 기울기가 $p$다. 그래서 $h$를 절반으로 줄여 오차가 $1/2^p$이 되는지 보면 방법의 차수를 실험으로 잴 수 있다. 시뮬레이터의 격자 수렴 검사가 정확히 이것이다.

안정성 — 알고리즘의 성질과 문제의 성질

같은 문제를 푸는 두 알고리즘이 반올림 오차에 다르게 반응할 수 있다. 오차가 단계마다 일정 비율로 줄거나 그대로면 안정, 단계마다 불어나면 불안정하다. 앞의 이차방정식에서 근의 공식 그대로가 불안정한 알고리즘이었고, $x_2 = c/x_1$이 안정한 알고리즘이었다.

그런데 알고리즘을 아무리 잘 짜도 안 되는 문제가 있다. 입력이 $10^{-16}$만 흔들려도 답이 $10^{-6}$만큼 바뀌는 문제라면, 어떤 알고리즘도 여섯 자리 이상은 못 준다. 이것은 알고리즘이 아니라 문제 자체의 성질이고, 7절에서 조건수라는 이름으로 정량화한다. 두 개념을 나눠 두자 —

• 조건(conditioning)은 문제의 성질. 입력의 작은 변화에 답이 얼마나 민감한가.
• 안정성(stability)은 알고리즘의 성질. 계산 중에 생긴 오차를 얼마나 키우는가.

안정한 알고리즘으로 조건이 좋은 문제를 풀 때만 정확한 답이 나온다. 하나라도 빠지면 답의 자릿수가 깎인다.

의문

$u \approx 10^{-16}$인데 왜 시뮬레이션 결과는 $10^{-3}$ 수준에서 흔들리는가?

오차의 근원은 반올림이 아니라 절단인 경우가 대부분이다. 격자 간격 $h$가 만드는 이산화 오차는 $h^2$ 수준이고, 소자 시뮬레이션의 격자는 계산량 때문에 결코 촘촘하지 않다. 그 위에 조건수가 $10^{6}$쯤 되는 선형계를 풀면(7절) 반올림 $10^{-16}$이 $10^{-10}$으로 불어나고, 뉴턴 반복을 $10^{-5}$에서 끊으면(2절) 그 정밀도로 끝난다. 세 층이 겹쳐서 결국 세 자리다.

이 말은 거꾸로 정밀도를 올리고 싶으면 어디를 손대야 하는지도 말해 준다. 4배 정밀도를 켜도 소용없다. 격자를 촘촘히 하거나(절단 오차), 방정식의 스케일을 맞추거나(조건수), 수렴 판정을 조이는(반복 오차) 것이 순서다.

반복법과 뉴턴-랩슨

왜 반복인가

$f(x) = 0$을 푸는 일은 수치해석의 첫 문제이면서 가장 자주 만나는 문제다. 대수 방정식은 5차부터 근의 공식이 없고, $x = \cos x$ 같은 초월 방정식은 처음부터 닫힌 해가 없다. 그래서 답을 한 번에 구하는 대신, 대충 찍은 뒤 고쳐 나간다. 이것이 반복법(iteration)이고, 이 장의 6절·9절·10절·14절이 전부 같은 뼈대다.

반도체에서는 페르미 준위를 구하는 전하 중립 조건 $n(E_F) - p(E_F) + N_A^- - N_D^+ = 0$이 초월 방정식이고, 푸아송 방정식의 우변 $n_i(e^{\psi/V_T} - e^{-\psi/V_T})$가 미지수의 지수함수라 방정식 전체가 비선형이다. 문턱전압 아래의 반전층 전하, 이종접합의 밴드 정렬, 쇼트키 접합의 영상 전하 — 어디를 봐도 $f(x) = 0$이 있다.

고정점 반복

방정식을 $x = g(x)$ 꼴로 바꿔 놓고 반복한다.

$$x_{n+1} = g(x_n), \qquad n = 0,\ 1,\ 2,\ \ldots$$

수렴하면 극한 $s$는 $s = g(s)$를 만족하므로 원래 방정식의 해다. 바꾸는 방법은 여럿이다 — $x^2 - 3x + 1 = 0$은 $x = (x^2 + 1)/3$으로도, $x = 3 - 1/x$로도 쓸 수 있고 둘의 수렴은 다르다. 무엇이 갈라놓는가.

정리. $g$가 구간 $J$를 $J$ 안으로 보내고, $J$ 위에서 $\lvert g'(x) \rvert \le K < 1$이면, $J$ 안의 어느 초기값에서 시작하든 반복은 $J$ 안의 유일한 고정점으로 수렴한다.

증명은 평균값 정리 한 줄이다. $s = g(s)$를 빼면

$$x_{n+1} - s = g(x_n) - g(s) = g'(\xi_n)(x_n - s) \;\Longrightarrow\; \lvert x_{n+1} - s \rvert \le K \lvert x_n - s \rvert \le K^{n+1} \lvert x_0 - s \rvert$$

$K < 1$이면 오른쪽이 0으로 간다. 한 번 반복할 때마다 오차가 $K$배 — 이것을 선형 수렴이라 한다. $K = 0.9$이면 자릿수 하나 얻는 데 22번, $K = 0.1$이면 1번이다.

직관

$x = g(x)$의 해는 곡선 $y = g(x)$와 직선 $y = x$의 교점이다. 반복은 교점 근처에서 곡선의 기울기가 1보다 완만하면 안으로 말려 들어가고, 가파르면 튕겨 나간다. 기울기가 양수면 한쪽에서 계단처럼 다가가고, 음수면 좌우로 번갈아 뛰면서 다가간다(거미줄 그림). $g'(s) = 0$이면 $K$가 사실상 0이라 한 번에 여러 자리를 얻는다 — 다음의 뉴턴법은 이 조건을 일부러 만들어 낸 고정점 반복이다.

뉴턴법 — 접선을 따라 내려간다

$f(x_n)$과 $f'(x_n)$을 알면, $x_n$에서 곡선을 접선으로 바꿔 놓고 접선이 축과 만나는 곳을 다음 근사로 삼는다.

$$\boxed{\;x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}\;}$$

같은 식을 테일러 전개로도 얻는다. $0 = f(s) = f(x_n) + (s - x_n)f'(x_n) + \cdots$에서 2차 이상을 버리고 $s$에 대해 풀면 위 식이다. 접선은 1차 테일러 다항식이니 같은 말이다.

그림 1.12-2. 뉴턴법의 한 걸음. 현재 위치 xₙ에서 곡선 f(x)에 접선을 긋고, 접선이 가로축을 지나는 자리를 xₙ₊₁로 삼는다. 접선의 기울기가 f′(xₙ)이므로 가로축까지의 거리는 f(xₙ)/f′(xₙ)이다. 근 근처에서 곡선이 접선과 거의 같아지므로 두세 걸음이면 붙지만, 오른쪽처럼 f′ ≈ 0인 평평한 자리에서 출발하면 접선이 축과 아주 먼 곳에서 만나 엉뚱한 곳으로 튀어 나간다.직접 작도

예 — 제곱근. $f(x) = x^2 - a$이면

$$x_{n+1} = x_n - \frac{x_n^2 - a}{2x_n} = \frac12\left(x_n + \frac{a}{x_n}\right)$$

바빌로니아 사람들이 알던 공식이다. $a = 2$, $x_0 = 1$에서 시작하면 $1.5,\ 1.41667,\ 1.41422,\ 1.41421356\ldots$ — 네 번째에서 이미 열 자리가 맞는다.

이것이 뉴턴법의 특징인 2차 수렴이다. $g(x) = x - f(x)/f'(x)$로 놓고 미분하면 $g'(x) = f f'' / f'^2$이므로 근 $s$에서 $g'(s) = 0$이다. 앞의 정리에서 $K$가 0으로 내려앉는 경우이고, 2차 항까지 전개하면

$$\varepsilon_{n+1} \approx -\frac{f''(s)}{2f'(s)}\,\varepsilon_n^2$$

오차가 제곱으로 준다. 오차 $10^{-2}$가 다음에 $10^{-4}$, 그다음 $10^{-8}$, $10^{-16}$ — 반복마다 맞는 자릿수가 두 배가 된다. 배정도의 16자리는 네다섯 번이면 끝난다.

직관

2차 수렴이 왜 접선에서 나오는가. 접선은 곡선을 그 점에서 값과 기울기까지 맞춘 것이라, 둘의 차이는 $\tfrac12 f''(\xi)(x - x_n)^2$ — 거리의 제곱이다. 근이 $x_n$에서 $\varepsilon_n$만큼 떨어져 있으면 접선이 축과 만나는 곳과 진짜 근의 차이도 $\varepsilon_n^2$ 급이 된다. 곡선을 직선으로 갈음한 오차가 다음 걸음의 오차로 그대로 옮겨진 것이다.

파라미터로 읽으면 앞의 계수 $f''/2f'$가 수렴의 품질을 정한다. 곡률이 크고(큰 $f''$) 기울기가 완만할수록(작은 $f'$) 제곱은 되어도 계수가 커서 처음 몇 걸음이 느리다. 그리고 $f'(s) = 0$이면 계수가 무한대 — 2차 수렴이 깨진다. 아래 ‘실패하는 경우’ 참조.

실패하는 경우

뉴턴법은 근 가까이에서만 2차 수렴한다. 멀리서 시작하면 다음 셋 중 하나가 흔히 일어난다.

• 평평한 곳에서 출발. $f'(x_n) \approx 0$이면 나눗셈이 큰 값을 만들어 멀리 튄다. 그림의 오른쪽 경우다. $\tan^{-1} x = 0$을 $x_0 = 1.5$에서 시작하면 발산한다.
• 순환. $x_0 \to x_1 \to x_0$처럼 두 점 사이를 오간다. $x^3 - 2x + 2 = 0$을 $x_0 = 0$에서 시작하면 0과 1을 영원히 오간다.
• 중근. $f(s) = 0$이면서 $f'(s) = 0$인 경우(예: $(x-1)^2$). 수렴은 하지만 $\varepsilon_{n+1} \approx \tfrac12\varepsilon_n$ — 선형 수렴으로 떨어진다. 중복도가 $m$이면 $x_{n+1} = x_n - m f/f'$로 고쳐 2차를 되찾을 수 있다.

실용적인 처방은 두 가지다. 첫째, 이분법(bisection)으로 먼저 근을 가둔다. $f(a)f(b) < 0$인 구간의 중점을 취해 부호가 바뀌는 쪽 절반을 남기는 방법으로, 한 번에 자릿수 하나 얻기도 어려울 만큼 느리지만(오차가 매번 정확히 절반) 절대 실패하지 않는다. 둘째, 감쇠 (damping). 뉴턴 걸음 $\Delta x = -f/f'$를 그대로 밟지 않고 $\lvert f \rvert$가 줄어들 때까지 $\Delta x/2,\ \Delta x/4,\ \ldots$로 줄여 밟는다.

주의

소자 시뮬레이터의 뉴턴 감쇠는 물리에서 온다. 푸아송 방정식의 미지수인 퍼텐셜 $\psi$가 한 걸음에 열역학 전압 $V_T = kT/q \approx 26\ \mathrm{mV}$의 몇 배 이상 움직이면, 우변의 $e^{\psi/V_T}$가 $e^{10}$ 이상 뛰어 접선 근사가 무의미해진다. 그래서 걸음을 $V_T$의 몇 배로 잘라 내는(clamping) 것이 표준이고, 그 뒤에 남는 것이 순수한 2차 수렴이다. 시뮬레이터 로그에서 처음 몇 번은 잔차가 천천히 줄다가(감쇠 구간) 어느 순간 $10^{-3},\ 10^{-6},\ 10^{-12}$로 곤두박질치는 것이 이 두 단계다.

할선법 — 도함수 없이

$f'$을 구하기 어렵거나 비싸면 접선 대신 최근 두 점을 잇는 할선(secant)을 쓴다.

$$x_{n+1} = x_n - f(x_n)\,\frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})}$$

뉴턴 공식의 $f'(x_n)$을 차분으로 바꾼 것이다. 수렴 차수는 $\tfrac12(1 + \sqrt5) \approx 1.618$ — 2차보다는 느리지만 한 걸음에 함수값 하나만 새로 계산하므로, 함수값 하나가 비쌀 때는 뉴턴법보다 시간당 빠른 경우가 많다.

$x_n - x_{n-1}$와 $f(x_n) - f(x_{n-1})$이 근 근처에서 둘 다 작아지므로 1절의 상쇄가 난다. 그래서 근에 충분히 가까워지면 더 이상 정밀도가 오르지 않는데, 이때는 이미 필요한 자릿수를 얻은 뒤라 실용상 문제가 되지 않는다.

연립 방정식으로 — 야코비 행렬

미지수가 여럿이면 $\mathbf{F}(\mathbf{x}) = \mathbf{0}$, 즉 $F_1(x_1, \ldots, x_N) = 0,\ \ldots,\ F_N = 0$이다. 테일러 전개를 벡터로 옮기면 접선이 접평면이 되고, 도함수가 야코비 행렬 $J_{ij} = \partial F_i / \partial x_j$가 된다(1.1.5이중·삼중 적분).

$$\mathbf{F}(\mathbf{x}^{(m)}) + J(\mathbf{x}^{(m)})\,\Delta\mathbf{x} = \mathbf{0} \;\Longrightarrow\; J(\mathbf{x}^{(m)})\,\Delta\mathbf{x} = -\mathbf{F}(\mathbf{x}^{(m)}), \qquad \mathbf{x}^{(m+1)} = \mathbf{x}^{(m)} + \Delta\mathbf{x}$$

나눗셈 $f/f'$가 선형계 풀이 $J^{-1}\mathbf{F}$로 바뀐다. 뉴턴 한 걸음마다 $N \times N$ 선형계를 하나씩 풀어야 하고, 그 선형계를 어떻게 푸는가가 5·6절이다. 수렴 성질은 스칼라와 같다 — 해 근처에서 2차, 멀리서는 감쇠 필요.

직관

소자 시뮬레이션의 뼈대가 정확히 이 식이다. 격자점이 $N$개이고 미지수가 점마다 셋($\psi,\ n,\ p$)이면 $\mathbf{x}$는 $3N$차원이고, $\mathbf{F}$는 이산화한 푸아송·전자 연속·정공 연속 방정식이다. $J$는 $3N \times 3N$이지만 격자점은 이웃과만 결합하므로 한 행에 0이 아닌 성분이 열 개 남짓인 희소 행렬이고, 그래서 $N = 10^{5}$이라도 풀 수 있다.

파라미터를 읽으면, $J$의 대각 블록은 각 방정식이 자기 미지수에 얼마나 민감한지, 비대각 블록은 방정식 사이의 결합이다. 푸아송의 우변이 $n$에, 연속 방정식의 드리프트 항이 $\psi$에 걸리는 것이 그 결합이고, 결합이 강할수록(높은 주입, 큰 바이어스) 한꺼번에 풀지 않으면 수렴하지 않는다.

세 방정식을 한꺼번에 푸는 대신 번갈아 풀 수도 있다. 푸아송을 $n,\ p$ 고정으로 풀어 $\psi$를 얻고, 그 $\psi$로 두 연속 방정식을 풀어 $n,\ p$를 갱신하고, 다시 푸아송으로 — 이것이 검멜 반복(Gummel iteration)이다. 야코비 행렬의 비대각 블록을 버린 셈이라 한 걸음이 싸지만 수렴은 선형이고, 결합이 강하면 발산한다. 시뮬레이터가 낮은 바이어스에서 검멜로 시작해 결합 뉴턴으로 넘어가는 것이 이 때문이다. 6절에서 이 “번갈아 풀기”를 선형계 위에서 다시 만난다.

언제 멈추는가

참값을 모르니 $\lvert x_n - s \rvert$로는 멈출 수 없다. 대신 쓰는 것이 둘이다.

• 갱신량 $\lvert x_{n+1} - x_n \rvert$. 2차 수렴 중이면 이것이 오차의 좋은 추정이다 — 다음 오차가 이번 오차의 제곱이므로, 이번 걸음의 크기가 곧 남은 오차의 크기다.
• 잔차 $\lvert f(x_n) \rvert$. 방정식이 얼마나 안 맞는가. 그러나 $f'$이 작으면 잔차가 작아도 오차는 클 수 있다(7절의 조건수 문제가 여기서도 같은 얼굴로 나타난다).

둘 다 절대값과 상대값을 함께 봐야 하는 이유는 1절의 RELTOL·ABSTOL 이야기와 같다. 미지수가 $10^{20}$과 $10^{2}$를 함께 담고 있으면 하나의 허용치로는 판정이 안 된다.

보간과 스플라인

점 몇 개로 함수를 대신한다

측정값은 언제나 점이다. SIMS로 잰 도핑 프로파일은 깊이 200군데의 농도이고, 이동도 모델은 온도·도핑 몇 조합에서의 표 값이며, 소자 모델의 표 기반(table-based) 모델은 $I$–$V$ 격자점의 전류다. 그 사이의 값이 필요할 때 점들을 잇는 함수를 만드는 일이 보간(interpolation)이다. 점 사이를 채우면 보간, 점 범위 밖으로 나가면 외삽(extrapolation)이고, 뒤의 것은 훨씬 위험하다.

$n + 1$개의 점 $(x_0, f_0),\ \ldots,\ (x_n, f_n)$이 주어졌다고 하자($x_j$는 서로 다르다). 차수가 $n$ 이하인 다항식 중 이 점을 전부 지나는 것은 정확히 하나 있다. 있다는 것은 아래 라그랑주 공식이 직접 보여 주고, 하나뿐이라는 것은 두 개가 있다면 그 차가 $n$차 이하이면서 근을 $n + 1$개 가지므로 영다항식이어야 한다는 데서 나온다.

라그랑주 형식

$k$번째 마디에서 1이고 나머지 마디에서 0인 $n$차 다항식을 만들면 끝이다.

$$L_k(x) = \prod_{\substack{j=0 \\ j \neq k}}^{n} \frac{x - x_j}{x_k - x_j}, \qquad L_k(x_j) = \begin{cases} 1 & j = k \\ 0 & j \neq k \end{cases}$$

분자는 $x_k$를 뺀 모든 마디에서 0이 되게 하고, 분모는 $x = x_k$에서 값을 1로 맞춘다. 이것을 함수값으로 가중해 더한다.

$$p_n(x) = \sum_{k=0}^{n} f_k\, L_k(x)$$

각 마디에서 자기 항만 살아남으므로 $p_n(x_k) = f_k$다. 식이 명시적이라 이론에 편하고, 4절의 적분 공식이 전부 여기서 나온다.

직관

$L_k$는 기저다. 1.2의 언어로 말하면 $n$차 이하 다항식의 공간은 $n + 1$차원이고, 단항식 $1, x, \ldots, x^n$ 대신 $L_0, \ldots, L_n$을 기저로 잡은 것이다. 이 기저의 좌표가 그냥 함수값 $f_k$ 자체라는 점이 라그랑주 형식의 요점이다 — 계수를 구하려고 연립방정식을 풀 필요가 없다. 단항식 기저에서 같은 일을 하려면 방데르몽드 행렬을 풀어야 하고, 그 행렬은 조건이 매우 나쁘다(7절).

오차 공식 — 어디서 얼마나 틀리는가

$f$가 $n + 1$번 미분 가능하면 마디가 아닌 $x$에서

$$f(x) - p_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}\,(x - x_0)(x - x_1)\cdots(x - x_n)$$

여기서 $\xi$는 $x$와 마디들을 포함하는 구간 안의 어딘가다. 정확한 위치는 모르고, 그래서 쓸 때는 $\lvert f^{(n+1)} \rvert$의 최댓값으로 갈음한다. 이 식은 테일러 정리의 나머지항과 꼭 닮았다 — 실제로 마디를 전부 한 점으로 모으면 테일러 나머지항이 된다(1.5.3함수의 급수 전개와 근사의 정확도).

두 인자를 따로 읽는다.

• $f^{(n+1)}(\xi)$ — 함수의 성질. 함수가 매끄러울수록 작다. $e^x$는 모든 도함수가 $e^x$라 얌전하고, $1/(1 + 25x^2)$는 고계 도함수가 $n!$처럼 자라 사납다.
• $\prod (x - x_j)$ — 마디의 배치. 마디 근처에서는 0에 가깝고, 마디 사이 한가운데에서 크며, 양 끝 구간에서 가장 크다. 마디 밖(외삽)으로 나가면 다항식처럼 빠르게 커진다.

균등 마디에서 두 번째 인자의 크기를 $h$로 적으면 $\lvert \prod(x - x_j) \rvert \le$ 대략 $n!\,h^{n+1}/4$ 수준이므로, 오차는 $h^{n+1}$에 비례한다. 선형 보간($n = 1$)이면 $h^2$ — 마디 간격을 절반으로 줄이면 오차가 4분의 1이다.

뉴턴의 분할 차분 — 마디를 하나씩 더할 수 있게

라그랑주 형식은 마디를 하나 추가하면 모든 $L_k$를 다시 만들어야 한다. 뉴턴 형식은 그 불편을 없앤다. 먼저 분할 차분을 정의한다.

$$f[x_0, x_1] = \frac{f_1 - f_0}{x_1 - x_0}, \qquad f[x_0, x_1, x_2] = \frac{f[x_1, x_2] - f[x_0, x_1]}{x_2 - x_0}, \qquad \ldots$$

1차 분할 차분은 두 점을 잇는 직선의 기울기, 2차는 기울기의 변화율이다 — 차례로 도함수 $f',\ f''/2!,\ f'''/3!,\ \ldots$에 대응한다. 그러면

$$\begin{aligned} p_n(x) = f_0 &+ (x - x_0)\,f[x_0, x_1] + (x - x_0)(x - x_1)\,f[x_0, x_1, x_2] \\ &+ \cdots + (x - x_0)\cdots(x - x_{n-1})\,f[x_0, \ldots, x_n] \end{aligned}$$

테일러 급수의 이산판이다. $(x - x_0)^k$ 자리에 $(x - x_0)\cdots(x - x_{k-1})$이, $f^{(k)}/k!$ 자리에 $k$차 분할 차분이 들어갔다. 마디를 하나 더 넣으면 항 하나만 뒤에 붙는다. 계산은 표로 한다 — 첫 열에 $f_j$, 둘째 열에 인접 1차 차분, 셋째 열에 인접 2차 차분 … 표의 맨 윗줄이 그대로 계수다.

마디가 균등($x_j = x_0 + jh$)하면 분할 차분의 분모가 상수가 되어 전진 차분 $\Delta f_j = f_{j+1} - f_j$, $\Delta^2 f_j = \Delta f_{j+1} - \Delta f_j$로 간단해진다. $r = (x - x_0)/h$로 놓으면 Newton 전진 공식은

$$p_n(x) = f_0 + r\,\Delta f_0 + \frac{r(r-1)}{2!}\,\Delta^2 f_0 + \cdots + \frac{r(r-1)\cdots(r-n+1)}{n!}\,\Delta^n f_0 = \sum_{s=0}^{n} \binom{r}{s}\,\Delta^s f_0$$

이항계수 $\binom{r}{s}$에 정수가 아닌 $r$이 들어가는 것이 포인트다. 표의 끝에서 거꾸로 세는 후진 공식도 있는데, Kreyszig가 그것을 $\nabla$로 적기 때문에 이 백과는 식을 옮기지 않는다(장 앞 주의 ④). 필요하면 마디 번호를 뒤집어 전진 공식으로 쓰면 같은 결과다.

룽게 현상 — 차수를 올리면 나빠진다

오차 공식대로면 마디를 늘릴수록 $h^{n+1}$이 작아져 좋아질 것 같다. 그렇지 않다. $f(x) = 1/(1 + 25x^2)$를 $[-1, 1]$의 균등 마디로 보간하면 $n$을 올릴수록 양 끝에서 진동이 커져 발산한다. 오차 공식의 첫 인자 $f^{(n+1)}$이 $n!$보다 빨리 자라기 때문이다.

그림 1.12-3. f(x) = 1/(1 + 25x²)를 [−1, 1]에서 보간한 것. 균등한 11개 마디로 10차 다항식을 지나게 하면 가운데는 잘 맞지만 양 끝 구간에서 크게 출렁인다(룽게 현상). 같은 함수를 같은 마디로 3차 스플라인으로 이으면 출렁임이 없다. 문제는 점의 수가 아니라 '한 다항식으로 전부를 맞추려는' 데 있다 — 한 곳의 정보가 다항식을 통해 멀리까지 전달되기 때문이다.직접 작도

직관

왜 양 끝인가. 인자 $\prod(x - x_j)$는 마디 사이에서 번갈아 부호가 바뀌는 $n + 1$차 다항식인데, 그 크기가 구간 가운데보다 끝 구간에서 훨씬 크다. 가운데의 $x$는 양쪽 마디에 둘러싸여 곱이 서로 상쇄되지만, 끝의 $x$는 한쪽으로만 마디가 있어 상쇄가 없다. 마디를 양 끝에 몰아 놓으면 — 정확히는 $x_j = \cos\bigl(\tfrac{2j+1}{2n+2}\pi\bigr)$의 체비쇼프 마디로 잡으면 — 이 인자의 최댓값이 균등 배치보다 지수적으로 작아지고, 룽게 현상이 사라진다. 그래서 고차 다항식 보간은 마디를 고를 수 있을 때(적분·스펙트럼법)만 쓰고, 마디가 측정으로 정해져 있을 때는 아래의 스플라인으로 간다.

스플라인 — 낮은 차수를 이어 붙인다

한 다항식으로 전부를 맞추는 대신, 마디 사이마다 3차 다항식을 두고 마디에서 매끄럽게 잇는다. 제도사가 가는 자(spline)를 못에 걸어 휘어 곡선을 그리던 도구에서 온 이름이다.

정의. 마디 $x_0 < x_1 < \cdots < x_n$에 대해 함수 $g(x)$가 3차 스플라인이라 함은 ① 각 구간 $[x_{j}, x_{j+1}]$에서 3차 다항식이고, ② 전체에서 $g,\ g',\ g''$이 연속이며, ③ $g(x_j) = f_j$인 것이다.

미지수를 세어 보자. 구간 $n$개 × 계수 4 = $4n$. 조건은 마디값 $2n$(각 구간 양 끝), 안쪽 마디 $n - 1$개에서 $g'$ 연속과 $g''$ 연속 $2(n-1)$. 합이 $4n - 2$이니 둘이 모자란다. 양 끝에서 조건 두 개를 더 주어야 하고, 그 선택이 스플라인의 종류다.

• 고정(clamped) 스플라인 — 양 끝의 기울기 $g'(x_0) = k_0$, $g'(x_n) = k_n$을 지정. Kreyszig의 기본 형태다.
• 자연(natural) 스플라인 — 양 끝에서 $g'' = 0$. 끝 너머로 직선처럼 나간다. 끝 기울기를 모를 때의 표준.

Kreyszig는 각 마디의 기울기 $k_j = g'(x_j)$를 미지수로 잡는다. 구간마다 양 끝의 값과 기울기로 3차식이 하나로 정해지므로(에르미트 보간), 남는 조건은 $g''$의 연속뿐이고, 마디 간격이 균등($h$)하면 그것이

$$k_{j-1} + 4k_j + k_{j+1} = \frac{3}{h}\,(f_{j+1} - f_{j-1}), \qquad j = 1, \ldots, n-1$$

이다. 미지수 $k_1, \ldots, k_{n-1}$에 대한 삼중대각 선형계다 — 한 줄에 미지수가 셋뿐이라 5절의 토머스 알고리즘으로 $O(n)$에 풀린다. 우변은 중심 차분 $(f_{j+1} - f_{j-1})/2h$의 6배, 즉 “기울기의 가중 평균이 중심 차분과 같다”는 조건이다. 대각 성분 4가 비대각 1 + 1보다 크므로 6절의 대각 우세 조건이 자동으로 성립해 언제나 풀린다.

직관

스플라인이 출렁이지 않는 이유는 물리적이다. 3차 스플라인은 마디를 지나는 모든 $C^2$ 곡선 중 굽힘 에너지 $\int (g'')^2\,dx$가 가장 작은 곡선이다 — 가는 자가 못 사이에서 취하는 모양이 바로 이것이다(작은 변형의 탄성 에너지가 곡률의 제곱에 비례하므로). 그래서 스플라인은 한 곳의 정보를 이웃 구간 너머로 거의 전달하지 않는다. 마디 하나를 옮겨도 그 영향은 몇 구간 안에서 지수적으로 사그라진다. 고차 다항식은 정반대다 — 마디 하나가 전체 계수를 바꾼다.

파라미터로 보면, 삼중대각 계의 비대각/대각 비율 $1/4$이 그 감쇠율이다. 한 마디의 교란은 이웃 마디로 갈 때마다 대략 $1/(2 + \sqrt3) \approx 0.27$배로 준다.

주의

측정값 보간에서 흔한 실수.

① 대수 축에서 보간할 것. 도핑 농도나 전류처럼 열 자릿수를 오가는 양을 선형 축에서 보간하면 큰 값 사이의 작은 값이 전부 뭉개진다. $\log N$을 보간하고 지수를 취한다. 스플라인이 음수를 뱉어 $\log$가 깨지는 사고도 이것으로 막힌다.

② 스플라인은 잡음도 정확히 지난다. 조건 ③이 “모든 점을 지난다”이므로 측정 잡음까지 충실히 따라가 마디 사이에서 물결을 만든다. 잡음이 있는 데이터는 보간이 아니라 8절의 최소제곱이 맞는 도구다. 점을 지나야 하면 보간, 점을 요약해야 하면 회귀다.

③ 외삽은 하지 않는다. 스플라인 밖으로 나가면 마지막 구간의 3차식이 그대로 이어져 순식간에 발산한다. 표 기반 소자 모델이 격자 밖 바이어스에서 이상한 전류를 내는 원인 중 하나다.

수치 적분과 미분

적분 — 보간 다항식을 대신 적분한다

$\int_a^b f(x)\,dx$를 닫힌 꼴로 못 구하면, $f$를 마디 몇 개에서 값을 재고 그 점을 지나는 보간 다항식을 대신 적분한다. 3절의 라그랑주 형식을 넣으면

$$\int_a^b f(x)\,dx \approx \int_a^b p_n(x)\,dx = \sum_{k=0}^{n} f_k \int_a^b L_k(x)\,dx = \sum_{k=0}^{n} w_k\, f_k$$

가중치 $w_k$는 $f$와 무관하게 마디 배치만으로 정해진다. 그래서 모든 수치 적분 공식은 “함수값의 가중합”이고, 공식마다 다른 것은 마디와 가중치뿐이다.

사다리꼴과 심프슨

구간을 $n$등분($h = (b-a)/n$)하고 각 조각에서 $f$를 직선으로 잇는 것이 사다리꼴 공식이다.

$$J \approx h\left[\tfrac12 f_0 + f_1 + f_2 + \cdots + f_{n-1} + \tfrac12 f_n\right]$$

각 조각의 오차는 3절의 보간 오차 공식을 적분해 $-\tfrac{1}{12}h^3 f''(\xi_j)$이고, $n$개를 모으면 $n h^3 = (b-a)h^2$이므로

$$\varepsilon_{\text{trap}} = -\frac{b-a}{12}\,h^2\,f''(\xi)$$

$h^2$ 차수. 기울기 부호는 위로 볼록한 함수($f'' < 0$)에서 사다리꼴이 곡선 아래를 과소평가 한다는 그림과 맞는다.

각 조각 두 개(마디 셋)를 포물선으로 잇는 것이 심프슨 공식이다. $n$은 짝수여야 한다.

$$J \approx \frac{h}{3}\left[f_0 + 4f_1 + 2f_2 + 4f_3 + \cdots + 2f_{n-2} + 4f_{n-1} + f_n\right]$$

$1,\ 4,\ 2,\ 4,\ \ldots,\ 4,\ 1$의 가중치는 포물선 조각 하나의 적분 $\tfrac{h}{3}(f_0 + 4f_1 + f_2)$을 이어 붙인 결과다. 오차는

$$\varepsilon_{\text{Simp}} = -\frac{b-a}{180}\,h^4\,f^{(4)}(\xi)$$

포물선(2차) 보간이니 $h^3$을 기대할 텐데 **$h^4$**다. 포물선 조각의 대칭 때문에 3차 항의 오차가 저절로 사라지기 때문이다 — 심프슨 공식은 3차 다항식까지 정확히 적분한다. 같은 $h$에서 사다리꼴보다 $h^2$배 좋고, 1절의 V자 그림에서 왼쪽 직선이 기울기 2에서 4로 가팔라진 것이다.

직관

$h^4$가 실무에서 뜻하는 것. 마디를 두 배로 늘리면 사다리꼴은 오차가 $1/4$, 심프슨은 $1/16$이 된다. 오차를 $10^{-8}$로 내리려면 사다리꼴은 마디 $10^{4}$개, 심프슨은 $10^{2}$개 급이다. 함수값 하나가 비싼 계산(예: 각 점마다 시뮬레이션 한 번)에서 이 차이는 100배의 시간이다.

단, 이 차수는 $f^{(4)}$가 유한할 때의 이야기다. 적분 구간 안에 꺾임이나 특이점이 있으면 (예: $\sqrt{x}$의 $x = 0$) 고계 도함수가 발산해 심프슨의 이점이 사라지고 사다리꼴과 같은 차수로 떨어진다. 그럴 때는 특이점에서 구간을 나누거나 변수를 바꿔 특이점을 없앤 뒤 적분한다.

리처드슨 외삽과 롬베르크

오차의 차수를 안다는 것은 오차를 빼낼 수 있다는 뜻이다. 사다리꼴 근사를 걸음 $h$와 $h/2$로 두 번 계산해 $J_h,\ J_{h/2}$를 얻었다면, 오차가 $Ch^2$ 꼴이므로

$$J - J_h \approx Ch^2, \qquad J - J_{h/2} \approx \tfrac14 Ch^2 \;\Longrightarrow\; J \approx J_{h/2} + \frac{J_{h/2} - J_h}{3}$$

두 근사를 조합해 $h^2$ 오차를 없앴다. 이것이 리처드슨 외삽이고, 결과는 정확히 심프슨 공식과 같다. 같은 요령을 되풀이해 $h^4,\ h^6,\ \ldots$ 오차를 차례로 없애는 것이 롬베르크 적분이다. 사다리꼴 값을 $h,\ h/2,\ h/4,\ \ldots$로 표의 첫 열에 적고, 이웃 열을 조합해 다음 열을 만든다.

$$T_{m,k} = T_{m,k-1} + \frac{T_{m,k-1} - T_{m-1,k-1}}{4^k - 1}$$

표의 대각선이 빠르게 수렴한다. “오차의 차수를 알면 그 차수를 없앨 수 있다”는 원리는 10절의 걸음 폭 조절(RKF)과 11절의 격자 수렴 검사에서 다시 쓴다.

가우스 구적법 — 마디까지 최적화한다

지금까지는 마디를 균등하게 미리 정해 놓고 가중치만 골랐다. 마디도 자유롭게 고르면 미지수가 두 배($n$개 마디 + $n$개 가중치 = $2n$)이므로, $2n - 1$차 다항식까지 정확히 적분하는 공식을 만들 수 있다. 그 최적 마디가 르장드르 다항식 $P_n$의 영점이고(1.5.6르장드르 방정식과 르장드르 다항식), 그 공식이 가우스 구적법(가우스–르장드르)이다.

$$\int_{-1}^{1} f(t)\,dt \approx \sum_{j=1}^{n} A_j\, f(t_j), \qquad P_n(t_j) = 0$$

$n = 2$이면 $t_j = \pm 1/\sqrt3$, $A_j = 1$ — 함수값 두 개로 3차 다항식까지 정확하다. 심프슨은 세 개로 3차이니, 같은 정확도를 값 하나 덜 쓰고 얻는다. $n = 5$면 값 다섯 개로 9차. 일반 구간 $[a, b]$는 $x = \tfrac12(a+b) + \tfrac12(b-a)t$로 옮긴다.

왜 르장드르인가. 오차는 $f$를 $2n - 1$차 다항식 $q$와 나머지로 쪼개서 보면 $\int f - \sum A_j f(t_j) = \int (f - q)$ 꼴이 되고, $f - q$가 마디에서 0이 되게 하면 그것은 $\prod(t - t_j) \times (\text{저차 다항식})$이다. $\prod(t - t_j)$가 $n$차 이하의 모든 다항식과 직교하면 이 적분이 사라진다 — 그 직교 다항식이 $P_n$이다. 1.6의 직교 함수계가 적분 공식으로 다시 나온 것이다.

직관

가우스 구적법은 마디를 양 끝으로 몰아 놓는다. $P_n$의 영점은 가운데보다 $\pm1$ 근처에 빽빽하고, 이것은 3절의 체비쇼프 마디와 같은 이유다 — 끝 구간에서 보간 오차 인자가 크니 거기에 마디를 더 둔다. 그래서 가우스 공식은 룽게 현상 없이 $n$을 올릴 수 있고, 매끄러운 함수에서는 $n$에 대해 지수적으로 수렴한다. 균등 마디 공식이 $h^p$로 대수적으로 수렴하는 것과 급이 다르다.

가중 함수를 바꾸면 다른 직교 다항식이 나온다. $\int_0^\infty e^{-x} f(x)\,dx$는 라게르 다항식(가우스–라게르), $\int_{-\infty}^{\infty} e^{-x^2} f(x)\,dx$는 에르미트 다항식이 마디를 준다(1.5.11에르미트·라게르 다항식). 페르미–디랙 적분 $F_{1/2}(\eta) = \int_0^\infty \dfrac{\sqrt{x}\,dx}{1 + e^{x - \eta}}$처럼 무한 구간에 지수 꼬리가 있는 적분이 정확히 이 부류다 — 축퇴 반도체의 캐리어 농도를 계산할 때마다 이 적분을 수치로 한다.

적응 적분

함수가 어떤 곳에서는 급하고 어떤 곳에서는 완만하면 균등 마디는 낭비다. 각 구간에서 오차를 추정해(예: 심프슨을 $h$와 $h/2$로 두 번) 허용치보다 크면 그 구간만 반으로 나누는 것이 적응 적분이다. 재귀적으로 진행되어 급한 곳에 마디가 저절로 몰린다. 라이브러리의 quad 류 함수가 이것이고, 오차 허용치를 인자로 받는 이유가 이 추정 때문이다.

수치 미분 — 적분의 반대, 성질도 반대

미분은 정의 그대로 차분으로 바꾼다. 테일러 전개에서 바로 나오는 세 공식이 기본이다.

$$\begin{aligned} f'(x) &\approx \frac{f(x+h) - f(x)}{h} &&\text{(전진, 오차 } \tfrac12 h f'')\\[4pt] f'(x) &\approx \frac{f(x+h) - f(x-h)}{2h} &&\text{(중심, 오차 } \tfrac16 h^2 f''')\\[4pt] f''(x) &\approx \frac{f(x+h) - 2f(x) + f(x-h)}{h^2} &&\text{(중심 2계, 오차 } \tfrac1{12} h^2 f^{(4)}) \end{aligned}$$

중심 차분이 전진보다 한 차수 좋은 것은 $f(x \pm h)$의 전개에서 $h^2$ 항이 대칭이라 빼면 사라지기 때문이다. 셋째 줄이 11절 차분법의 재료다 — 라플라시안의 5점 공식은 이 식을 $x$와 $y$로 한 번씩 쓴 것이다.

문제는 1절의 V자 그림이다. 적분은 함수값의 오차를 평균 내어 줄이지만, 미분은 함수값의 오차를 $h$로 나누어 키운다. 전진 차분은 최선이 $\sqrt{u} \approx 10^{-8}$, 중심 차분은 $u^{2/3} \approx 10^{-11}$, 2계 중심 차분은 $\sqrt{u} \approx 10^{-8}$ 급에서 멈춘다.

직관

스펙트럼으로 보면 이유가 한 줄이다. 미분은 진동수 $\omega$ 성분에 $i\omega$를 곱하고 (1.10.4푸리에 변환의 성질), 적분은 $1/i\omega$를 곱한다. 잡음은 고주파에 몰려 있으니 미분은 잡음을 증폭하고 적분은 잡음을 억제한다. 측정 데이터를 미분하는 일은 그래서 언제나 평활화(smoothing)와 함께 간다 — 스플라인이나 국소 다항식으로 먼저 맞춘 뒤 그 곡선을 미분하거나, 사비츠키–골레이 필터로 미분과 평활을 한 번에 한다.

주의

반도체 계측에서 이 함정에 직접 걸리는 예가 C–V 도핑 프로파일이다. 공핍 근사에서

$$N(x) = \frac{C^3}{q\,\varepsilon_s A^2}\left(\frac{dC}{dV}\right)^{-1}, \qquad x = \frac{\varepsilon_s A}{C}$$

이므로 측정한 $C(V)$를 미분해야 농도가 나온다. $C$에 $0.1\%$의 잡음이 있으면 $dC/dV$에는 그것이 $1/\Delta V$배로 들어가고, $N$은 그 역수라 잡음이 그대로 프로파일의 톱니로 나타난다. $1/C^2$ 대 $V$를 그려 기울기를 읽는 관행은 이 미분을 직선 맞춤(8절)으로 바꾸어 잡음을 누르는 요령이다. 문턱전압을 $g_m$의 최댓값으로 정하는 방법, 서브스레숄드 기울기를 $\log I_D$의 미분으로 읽는 방법도 전부 같은 문제를 안고 있다.

가우스 소거와 LU 분해

선형계는 어디에나 있다

이 장의 절반은 결국 $A\mathbf{x} = \mathbf{b}$로 끝난다. 2절의 뉴턴 걸음, 3절의 스플라인, 8절의 정규방정식, 11절의 차분 방정식, 13절의 유한요소 방정식이 전부 선형계이고, 그것을 얼마나 빨리 정확히 푸는가가 시뮬레이션의 속도를 정한다. 1.2에서 가우스 소거를 개념으로 봤다면(1.2.2가우스 소거와 선형계의 해), 여기서는 알고리즘으로 본다 — 몇 번 계산하고, 어디서 오차가 생기며, 어떻게 재활용하는가.

가우스 소거의 비용

앞선 소거로 $A$를 상삼각으로 만든 뒤 후진 대입한다. 연산 횟수를 세 보면, $k$번째 열을 소거할 때 아래 $(n - k)$개 행 각각에 $(n - k)$개 성분을 갱신하므로

$$\sum_{k=1}^{n-1} (n-k)^2 \approx \frac{n^3}{3} \quad(\text{곱셈·덧셈 각각}), \qquad \text{후진 대입} \approx \frac{n^2}{2}$$

소거가 $n^3$, 대입이 $n^2$. $n$이 두 배가 되면 시간이 여덟 배다. 크래머 공식이 $n!$이니 비교가 안 되고, 역행렬을 구해 곱하는 것은 소거의 세 배라 역행렬을 구하는 것은 언제나 손해다(1.2.4역행렬). $n = 10^{4}$이면 $10^{12}$번 — 초 단위, $n = 10^{6}$이면 $10^{18}$번 — 불가능. 소자 시뮬레이션의 $n$이 그 사이에 있고, 그래서 아래의 희소성이 없으면 아무것도 못 푼다.

피벗팅 — 작은 수로 나누지 않는다

수학적으로 가우스 소거는 피벗(대각 성분)이 0만 아니면 된다. 수치적으로는 작기만 해도 망한다. 유효숫자 4자리로

$$\begin{aligned} 0.0004\,x_1 + 1.402\,x_2 &= 1.406 \\ 0.4003\,x_1 - 1.502\,x_2 &= 2.501 \end{aligned}$$

을 풀어 보자(Kreyszig의 예). 첫 행을 피벗으로 쓰면 곱수가 $0.4003/0.0004 \approx 1001$이고, 둘째 행은 $-1.502 - 1001 \times 1.402 = -1405$가 되어 원래의 $-1.502$가 자릿수 밖으로 밀려난다. 결과 $x_2 = 1.002$, $x_1 = -2.500$ — 참값 $x_2 = 1$, $x_1 = 10$과 딴판이다. 두 행을 바꿔 $0.4003$을 피벗으로 쓰면 곱수가 $0.001$이 되어 정확히 $x_1 = 10.00$, $x_2 = 1.000$이 나온다.

부분 피벗팅은 매 단계에서 그 열의 남은 성분 중 절댓값이 가장 큰 것을 피벗으로 올리는 규칙이다. 곱수가 언제나 1 이하가 되어, 소거 과정에서 성분이 폭발적으로 자라는 일을 막는다. 행렬이 대각 우세이거나 대칭 양의 정부호이면 피벗팅 없이도 안전하다는 것이 알려져 있고, 그 경우가 아니면 언제나 켠다.

직관

피벗이 작다는 것은 그 행이 다른 행과 거의 평행하다는 뜻이다(행렬식이 작다). 그 행을 기준으로 다른 행들을 깎아 내면 큰 수를 곱해 빼게 되고, 1절의 상쇄가 일어난다. 피벗팅은 “가장 평행하지 않은 행부터 쓴다”는 규칙이고, 그렇게 해도 행렬 자체가 거의 특이하면 소용없다 — 그것은 알고리즘이 아니라 문제의 성질이라(7절) 어떤 순서로도 못 피한다.

LU 분해 — 소거를 저장한다

소거 과정에서 쓴 곱수를 버리지 않고 모으면, 소거 결과와 합쳐 원래 행렬의 곱 분해가 된다.

$$A = LU, \qquad L = \begin{bmatrix} 1 & & \\ l_{21} & 1 & \\ l_{31} & l_{32} & 1 \end{bmatrix}, \quad U = \begin{bmatrix} u_{11} & u_{12} & u_{13} \\ & u_{22} & u_{23} \\ & & u_{33} \end{bmatrix}$$

$L$은 곱수를 아래에 채운 단위 하삼각(대각이 1), $U$는 소거가 끝난 상삼각이다. 이렇게 대각을 $L$ 쪽에서 1로 고정한 것이 둘리틀(Doolittle) 분해, $U$ 쪽에서 1로 고정하면 크라우트 (Crout) 분해다. 둘은 대각행렬 하나 차이라 본질이 같다.

분해가 있으면 $A\mathbf{x} = \mathbf{b}$는 두 번의 삼각계로 쪼개진다.

$$L\mathbf{y} = \mathbf{b} \quad(\text{전진 대입}), \qquad U\mathbf{x} = \mathbf{y} \quad(\text{후진 대입})$$

각각 $n^2/2$이므로 **한 번 분해해 두면 우변마다 $n^2$**이다. 우변이 여러 개인 문제 — 같은 구조에 바이어스만 바꾸는 소자 시뮬레이션, 같은 회로에 입력만 바꾸는 과도 해석, 역행렬의 열을 하나씩 구하는 일 — 에서 이것이 $n^3$과 $n^2$의 차이다. 회로 시뮬레이터가 한 시간 걸음의 뉴턴 반복 안에서 야코비 행렬이 크게 변하지 않으면 분해를 재활용하는 것도 같은 이유다(수정 뉴턴법).

행렬식도 공짜다. $\det A = \det L \det U = u_{11}u_{22}\cdots u_{nn}$ — 피벗의 곱이고, 행 교환마다 부호가 바뀐다. 피벗 중 하나가 0에 가까우면 그 행렬이 특이에 가깝다는 경고다.

콜레스키 — 대칭 양의 정부호이면 절반

$A$가 대칭이고 양의 정부호(1.3.7이차 형식과 주축 변환)이면 $U = L^{\mathsf{T}}$로 잡을 수 있다.

$$A = LL^{\mathsf{T}}, \qquad l_{jj} = \sqrt{a_{jj} - \sum_{k<j} l_{jk}^2}, \qquad l_{ij} = \frac{1}{l_{jj}}\Bigl(a_{ij} - \sum_{k<j} l_{ik}l_{jk}\Bigr)$$

절반의 저장 공간과 절반의 연산($n^3/6$)으로 끝나고, 피벗팅이 필요 없다 — 제곱근 안이 언제나 양수라는 것이 양의 정부호의 정의이고, 그래서 성분이 자라지 않는다. 제곱근 안이 음수가 되면 그 행렬은 양의 정부호가 아니라는 판정이 되기도 한다. 13절의 유한요소 강성 행렬, 8절의 정규방정식 $A^{\mathsf{T}}A$, 11절의 (부호를 바꾼) 이산 라플라시안이 전부 이 부류라 콜레스키가 기본 도구다.

띠와 희소 — $n^3$을 피하는 유일한 길

이산화한 PDE의 행렬은 각 격자점이 이웃과만 결합하므로 대부분의 성분이 0이다. 1차원이면 삼중대각이고, 이때 소거는 각 행에서 한 성분만 없애면 되므로 전체가 $O(n)$이다 — 이것이 토머스 알고리즘이고, 3절의 스플라인과 12절의 크랭크–니콜슨이 여기에 올라탄다.

$$\begin{bmatrix} b_1 & c_1 & & \\ a_2 & b_2 & c_2 & \\ & \ddots & \ddots & \ddots \\ & & a_n & b_n \end{bmatrix} \quad\longrightarrow\quad c'_j = \frac{c_j}{b_j - a_j c'_{j-1}}, \qquad d'_j = \frac{d_j - a_j d'_{j-1}}{b_j - a_j c'_{j-1}}$$

한 번 전진해서 $c'_j,\ d'_j$를 만들고 한 번 후진해서 $x_j = d'_j - c'_j x_{j+1}$로 푼다.

일반적으로 대각에서 거리 $w$ 안에만 성분이 있으면(띠 폭 $w$) 소거 비용이 $n w^2$이다. 2차원 $N \times N$ 격자를 행 순서로 번호 매기면 $n = N^2$, $w = N$이므로 비용은 $N^4 = n^2$; 3차원은 $n = N^3$, $w = N^2$이라 $N^7 = n^{7/3}$이다. 차원이 오를수록 직접법의 상대적 비용이 급히 커지는데, 6절의 반복법이 필요한 자리가 정확히 여기다.

직관

희소 행렬의 소거에서 문제는 채워짐(fill-in)이다. 띠 안에 0이었던 자리가 소거 도중 0이 아닌 값으로 채워져, $L$과 $U$가 원래 $A$보다 훨씬 빽빽해진다. 채워짐의 양은 미지수의 번호 순서에 따라 달라지므로, 희소 직접 해법의 핵심은 소거 자체가 아니라 채워짐을 최소화하는 순서 재배열이다(최소 차수 순서, 중첩 분할 등). 회로 시뮬레이터의 KLU, 소자 시뮬레이터의 PARDISO·MUMPS가 이 일을 하는 라이브러리이고, “행렬 재배열” 단계가 로그에 따로 찍히는 이유다.

물리적으로 읽으면, 미지수 하나를 소거하는 것은 그 노드를 회로에서 지우고 그 노드에 연결됐던 모든 이웃을 서로 직접 연결하는 일이다(별–삼각 변환). 이웃이 많은 노드를 먼저 지우면 새 연결이 폭발한다 — 그래서 이웃이 적은 노드부터 지운다는 것이 최소 차수 순서의 뜻이다.

반복 해법과 수렴(야코비·가우스-자이델)

직접법이 안 될 때

5절의 직접법은 유한 번의 연산으로 (반올림을 빼면) 정확한 답을 준다. 그런데 3차원 소자를 $100 \times 100 \times 100$ 격자로 이산화하면 미지수가 $10^{6}$개이고, 띠 폭이 $10^{4}$라 직접법의 비용 $n w^2 = 10^{14}$은 감당이 안 된다. 반면 행렬–벡터 곱 $A\mathbf{x}$는 0이 아닌 성분 수만큼, 즉 $7 \times 10^{6}$번이면 된다. 행렬–벡터 곱만으로 답에 다가가는 방법이 반복법이고, 반복 횟수가 $10^{8}$번보다 훨씬 적으면 이긴다.

야코비와 가우스–자이델

$A\mathbf{x} = \mathbf{b}$의 $i$번째 방정식을 $x_i$에 대해 푼 꼴로 다시 쓴다.

$$x_i = \frac{1}{a_{ii}}\Bigl(b_i - \sum_{j \neq i} a_{ij}\,x_j\Bigr)$$

우변에 현재 추정 $\mathbf{x}^{(m)}$을 넣어 좌변의 새 추정 $\mathbf{x}^{(m+1)}$을 얻는 것이 야코비 반복이다. 한 바퀴(sweep)에 모든 $i$를 돈다. 가우스–자이델 반복은 한 가지만 다르다 — $x_i$를 갱신할 때 이미 갱신된 $x_1, \ldots, x_{i-1}$의 새 값을 즉시 쓴다.

$$x_i^{(m+1)} = \frac{1}{a_{ii}}\Bigl(b_i - \sum_{j < i} a_{ij}\,x_j^{(m+1)} - \sum_{j > i} a_{ij}\,x_j^{(m)}\Bigr)$$

야코비는 두 벡터를 들고 있어야 하고 순서가 상관없어 병렬화가 쉽다. 가우스–자이델은 한 벡터를 제자리에서 덮어쓰고, 대개 야코비의 두 배 빠르다.

Kreyszig는 대각을 1로 정규화한 뒤($a_{ii} = 1$이 되게 각 행을 나눈다) $A = I + L + U$로 쪼개 적는다 — 여기서 $L,\ U$는 5절의 LU 인자가 아니라 $A$의 순수한 아래·위 삼각 부분이다. 그러면

$$\begin{aligned} \text{야코비:}\quad &\mathbf{x}^{(m+1)} = \mathbf{b} - (L + U)\,\mathbf{x}^{(m)} \\ \text{가우스–자이델:}\quad &\mathbf{x}^{(m+1)} = \mathbf{b} - L\,\mathbf{x}^{(m+1)} - U\,\mathbf{x}^{(m)} \;\Longrightarrow\; \mathbf{x}^{(m+1)} = (I + L)^{-1}\bigl(\mathbf{b} - U\mathbf{x}^{(m)}\bigr) \end{aligned}$$

주의

이 절의 $L,\ U$는 5절의 $L,\ U$가 아니다. 5절에서는 $A = LU$(곱), 이 절에서는 $A = I + L + U$(합)다. Kreyszig가 두 절에서 같은 글자를 쓰므로 그대로 두되, 이 절 안에서는 언제나 합 분해다. 또 이 절의 반복 행렬을 $C$로 쓰는데, 3절의 스플라인 계수나 회로의 커패시턴스와 무관하다.

수렴 조건 — 스펙트럼 반지름

두 방법 모두 $\mathbf{x}^{(m+1)} = C\mathbf{x}^{(m)} + \mathbf{c}$ 꼴이다(야코비는 $C = -(L + U)$, 가우스–자이델은 $C = -(I + L)^{-1}U$). 참해 $\mathbf{x}$도 같은 식을 만족하므로 빼면 오차 $\mathbf{e}^{(m)} = \mathbf{x} - \mathbf{x}^{(m)}$에 대해

$$\mathbf{e}^{(m+1)} = C\,\mathbf{e}^{(m)} \;\Longrightarrow\; \mathbf{e}^{(m)} = C^m\,\mathbf{e}^{(0)}$$

$C^m \to O$일 필요충분조건은 $C$의 모든 고유값의 절댓값이 1보다 작은 것, 즉 스펙트럼 반지름 $\rho(C) < 1$이다(1.3.5닮음 변환과 대각화 — 대각화하면 $C^m = X D^m X^{-1}$이고 $D^m$의 성분이 $\lambda_j^m$). 수렴 속도는 가장 큰 고유값이 정한다.

$$\lVert \mathbf{e}^{(m)} \rVert \sim \rho(C)^m\, \lVert \mathbf{e}^{(0)} \rVert \;\Longrightarrow\; \text{자릿수 하나 얻는 데 } \frac{\ln 10}{-\ln \rho(C)} \text{ 번}$$

$\rho = 0.5$면 3.3번, $\rho = 0.99$면 230번, $\rho = 0.999$면 2300번이다. 반복법의 좋고 나쁨은 전부 $\rho$가 1에 얼마나 가까운가의 문제다.

고유값을 미리 구하기는 어려우므로 실제로는 충분조건을 쓴다. 어떤 행렬 노름(7절)에서든 $\lVert C \rVert < 1$이면 $\rho(C) \le \lVert C \rVert < 1$이다. 행합 노름으로 쓰면 가장 쉬운 판정이 나온다 —

$$\lvert a_{ii} \rvert > \sum_{j \neq i} \lvert a_{ij} \rvert \quad(\text{모든 } i)$$

즉 강한 대각 우세이면 야코비도 가우스–자이델도 수렴한다. 이산 라플라시안은 대각이 $-4$, 비대각 합이 4라 등호가 성립하는 경계 상황인데, 경계 조건이 있는 행 몇 개에서 부등호가 성립하면 그것으로 충분하다는 것이 알려져 있다(약한 대각 우세 + 기약).

얼마나 느린가 — 라플라스 방정식의 경우

$N \times N$ 격자의 5점 라플라시안(11절)에 대해 야코비 반복 행렬의 고유값은 정확히 알려져 있다.

$$\lambda_{pq} = \tfrac12\Bigl(\cos\frac{p\pi}{N+1} + \cos\frac{q\pi}{N+1}\Bigr), \qquad \rho_J = \cos\frac{\pi}{N+1} \approx 1 - \frac{\pi^2}{2N^2}$$

$N = 100$이면 $\rho_J \approx 0.9995$ — 자릿수 하나에 4700번, 여섯 자리에 3만 번이다. 가우스–자이델은 $\rho_{GS} = \rho_J^2$이라 절반, 그래도 만 번이 넘는다. 게다가 $N$을 두 배로 하면 $\rho$가 1에 네 배 가까워져 반복 횟수가 네 배가 되고, 한 바퀴 비용도 네 배이니 전체가 열여섯 배다.

직관

왜 이렇게 느린가. 가우스–자이델 한 바퀴는 각 점을 이웃의 평균으로 바꾸는 일이다 — 그것은 확산 방정식을 한 걸음 진행시키는 것과 같다(12절의 명시적 차분과 비교해 보라). 한 점에 잘못 놓인 오차는 한 바퀴에 이웃 한 칸씩만 번져 나가고, 확산이니 거리 $N$을 건너는 데 $N^2$ 바퀴가 든다. 반복 횟수 $\propto N^2$은 이 확산 시간이다.

고유값으로 읽으면, 오차의 고주파 성분(이웃끼리 부호가 다른 톱니)은 $\lambda_{pq} \approx 0$이라 한두 바퀴에 사라지고, 저주파 성분(격자 전체에 걸친 완만한 언덕)은 $\lambda \approx 1$이라 거의 줄지 않는다. 이웃 평균은 국소 연산이라 전역 구조를 못 본다. 다중 격자법(multigrid)은 이 관찰에서 나온 방법이다 — 저주파 오차를 성긴 격자로 옮기면 거기서는 고주파가 되어 빨리 사라지고, 그것을 되돌려 보정한다. 격자 수준을 재귀적으로 오르내리면 전체 비용이 $O(n)$, 반복법의 이론적 최선이다.

그림 1.12-4. 같은 라플라스 문제에 대한 세 반복법의 오차를 반복 횟수에 대해 반대수 축에 그린 것. 직선의 기울기가 log ρ(C)이므로 기울기가 가파를수록 빠르다. 가우스–자이델은 야코비의 정확히 두 배 기울기(ρ_GS = ρ_J²)이고, 최적 과완화(SOR)는 격자가 커질수록 두 방법과의 격차가 벌어져 반복 횟수가 N²이 아니라 N에 비례한다. 오차가 곧장 직선을 타지 않고 처음 몇 번 빨리 떨어지는 것은 고주파 오차 성분이 먼저 사라지는 구간이다.직접 작도

과완화 — 한 걸음을 과장한다

가우스–자이델이 한 걸음에 $\Delta x_i$만큼 고친다면, 그보다 더 밀어 본다.

$$x_i^{(m+1)} = x_i^{(m)} + \omega\,\Delta x_i^{GS}, \qquad 1 < \omega < 2$$

$\omega = 1$이면 가우스–자이델이고, $\omega > 1$이 연속 과완화(SOR, successive over-relaxation)다. $\omega \ge 2$면 발산한다. 최적 $\omega$는 야코비의 스펙트럼 반지름으로 정해진다.

$$\omega_{\text{opt}} = \frac{2}{1 + \sqrt{1 - \rho_J^2}}, \qquad \rho_{\text{SOR}} = \omega_{\text{opt}} - 1 \approx 1 - \frac{2\pi}{N}$$

$1 - \pi^2/2N^2$이 $1 - 2\pi/N$으로 — $1$에서의 거리가 $N^2$ 급에서 $N$ 급으로 좋아졌다. $N = 100$이면 $\rho_{\text{SOR}} \approx 0.94$, 자릿수 하나에 37번이다. 가우스–자이델의 4700번에 비하면 백 배가 넘는다. 단 $\omega$가 최적에서 조금만 벗어나도 이득이 급히 줄고, 최적값은 문제마다 다르므로 실제로는 $\rho_J$를 추정해 가며 맞춘다.

직관

과완화가 왜 되는가. 가우스–자이델의 오차는 언제나 같은 부호로 단조롭게 준다 — 이웃 평균이 참값을 향해 조심스럽게 다가가기 때문이다. 방향이 맞는데 걸음이 짧다면 걸음을 늘리면 된다. $\omega$가 걸음 배율이고, 너무 늘리면(2 이상) 참값을 지나쳐 반대편에서 발산한다. 14절의 최급강하법에서 걸음 폭을 고르는 문제가 정확히 같은 구조이고, 거기서도 “곡률의 역수”로 정해진다.

크릴로프 부분공간법 — 현대의 표준

야코비·가우스–자이델·SOR은 고정 반복이다 — 매번 같은 $C$를 곱한다. 현대의 반복법은 다르게 간다. 지금까지 만든 벡터들 $\mathbf{r},\ A\mathbf{r},\ A^2\mathbf{r},\ \ldots$이 펼치는 공간(크릴로프 부분공간) 안에서 가장 좋은 근사를 매번 고른다.

• 켤레 기울기법(CG) — $A$가 대칭 양의 정부호일 때. 잔차가 서로 $A$-직교가 되게 방향을 고르며, 반복 횟수가 $\sqrt{\kappa(A)}$에 비례한다(7절). 라플라시안이면 $\kappa \sim N^2$이라 $N$번 — SOR과 같은 급이지만 $\omega$ 같은 조율 파라미터가 없다. 14절에서 최적화로 다시 나온다.
• GMRES · BiCGSTAB — 비대칭 행렬용. 소자 시뮬레이션의 야코비 행렬은 드리프트 항 때문에 비대칭이라 이쪽을 쓴다.

어느 방법이든 성패는 전처리(preconditioning)에 달려 있다. $A\mathbf{x} = \mathbf{b}$ 대신 $M^{-1}A\mathbf{x} = M^{-1}\mathbf{b}$를 풀되, $M$은 $A$와 비슷하면서 풀기 쉬운 행렬로 잡는다 (불완전 LU, 야코비의 대각, 다중 격자 한 바퀴). $M^{-1}A$의 고유값이 1 근처로 뭉치면 반복이 몇 십 번에 끝난다. 시뮬레이터 설정의 ILU, AMG 같은 항목이 이 $M$의 선택이다.

의문

2절의 검멜 반복은 이 절의 무엇인가?

블록 가우스–자이델이다. 미지수를 $\psi$ 블록, $n$ 블록, $p$ 블록으로 묶고, 한 블록을 풀 때 다른 블록은 최신 값으로 고정한다. 그래서 수렴 조건도 같다 — 블록 사이의 결합(야코비 행렬의 비대각 블록)이 대각 블록에 비해 약해야 한다. 낮은 바이어스에서는 $\psi$가 $n,\ p$에 약하게 걸려 잘 수렴하고, 높은 주입에서는 결합이 강해져 $\rho > 1$이 되어 발산한다. 그때 결합 뉴턴으로 넘어가는 것은 “블록 반복을 포기하고 전체 행렬을 직접 푼다”는 뜻이다.

노름·조건수와 불량 조건계

크기를 재는 자 — 노름

“오차가 작다”고 말하려면 벡터의 크기를 재는 자가 있어야 한다. 1.2에서 내적으로 정의한 유클리드 길이(1.2.5벡터공간과 내적공간)가 그 하나이고, 수치해석은 셋을 나란히 쓴다.

$$\lVert \mathbf{x} \rVert_1 = \sum_j \lvert x_j \rvert, \qquad \lVert \mathbf{x} \rVert_2 = \Bigl(\sum_j x_j^2\Bigr)^{1/2}, \qquad \lVert \mathbf{x} \rVert_\infty = \max_j \lvert x_j \rvert$$

$\lVert \cdot \rVert_1$은 성분의 합(택시 거리), $\lVert \cdot \rVert_2$는 유클리드 길이, $\lVert \cdot \rVert_\infty$는 가장 큰 성분이다. 어느 것을 써도 “0에 간다”는 결론은 같다 — 유한 차원에서는 상수배 안에서 서로 갈음되기 때문이다($\lVert \mathbf{x} \rVert_\infty \le \lVert \mathbf{x} \rVert_2 \le \sqrt{n}\, \lVert \mathbf{x} \rVert_\infty$). 다만 값과 편의가 다르다. $\infty$-노름은 계산이 가장 싸고 “최악의 성분”을 말해 주며, 2-노름은 회전에 불변이라 이론에 편하다.

행렬의 노름은 벡터를 얼마나 늘일 수 있는가로 정의한다.

$$\lVert A \rVert = \max_{\mathbf{x} \neq \mathbf{0}} \frac{\lVert A\mathbf{x} \rVert}{\lVert \mathbf{x} \rVert}$$

벡터 노름마다 하나씩 딸려 나오고(유도 노름), 셋은 닫힌 꼴이 있다.

노름	값	읽는 법
$\lVert A \rVert_1$	$\max_j \sum_i \lvert a_{ij} \rvert$	열합의 최대
$\lVert A \rVert_\infty$	$\max_i \sum_j \lvert a_{ij} \rvert$	행합의 최대
$\lVert A \rVert_2$	$\sqrt{\rho(A^{\mathsf{T}}A)} = \sigma_{\max}$	최대 특이값. $A$ 대칭이면 $\max \lvert \lambda \rvert$
$\lVert A \rVert_F$	$\bigl(\sum_{i,j} a_{ij}^2\bigr)^{1/2}$	프로베니우스. 유도 노름은 아니지만 계산이 쉽다

유도 노름은 정의에서 바로 $\lVert A\mathbf{x} \rVert \le \lVert A \rVert\, \lVert \mathbf{x} \rVert$와 $\lVert AB \rVert \le \lVert A \rVert\, \lVert B \rVert$를 만족한다. 6절의 수렴 판정 $\lVert C \rVert < 1$이 이 부등식이었다. 그리고 어느 노름에서든 $\rho(A) \le \lVert A \rVert$다 — 고유벡터를 넣으면 $\lVert A\mathbf{x} \rVert = \lvert \lambda \rvert\, \lVert \mathbf{x} \rVert$이기 때문이다.

조건수 — 문제가 얼마나 예민한가

$A\mathbf{x} = \mathbf{b}$에서 우변이 $\delta\mathbf{b}$만큼 흔들리면 해는 얼마나 흔들리는가. $A(\mathbf{x} + \delta\mathbf{x}) = \mathbf{b} + \delta\mathbf{b}$에서 $\delta\mathbf{x} = A^{-1}\delta\mathbf{b}$이므로

$$\lVert \delta\mathbf{x} \rVert \le \lVert A^{-1} \rVert\, \lVert \delta\mathbf{b} \rVert, \qquad \lVert \mathbf{b} \rVert \le \lVert A \rVert\, \lVert \mathbf{x} \rVert$$

두 부등식을 나누면 상대 변화 사이의 관계가 나온다.

$$\boxed{\; \frac{\lVert \delta\mathbf{x} \rVert}{\lVert \mathbf{x} \rVert} \le \underbrace{\lVert A \rVert\, \lVert A^{-1} \rVert}_{\kappa(A)}\; \frac{\lVert \delta\mathbf{b} \rVert}{\lVert \mathbf{b} \rVert} \;}$$

$\kappa(A) = \lVert A \rVert\, \lVert A^{-1} \rVert$가 조건수다. 입력의 상대 오차가 답의 상대 오차로 최대 $\kappa$배 증폭된다. $A$가 $\delta A$만큼 흔들릴 때도 1차 근사에서 같은 꼴이다. 언제나 $\kappa \ge 1$이고($\lVert I \rVert = 1 = \lVert AA^{-1} \rVert \le \kappa$), 직교행렬은 $\kappa_2 = 1$로 최선이다. $\kappa$가 크면 불량 조건(ill-conditioned)이라 한다.

직관

자릿수로 읽는다. 입력이 배정도로 정확하면 상대 오차 $10^{-16}$, 답은 최악의 경우 $\kappa \times 10^{-16}$. 그러므로 $\log_{10}\kappa$자리를 잃는다. $\kappa = 10^{6}$이면 답은 열 자리, $\kappa = 10^{12}$이면 네 자리, $\kappa = 10^{16}$이면 답에 맞는 자리가 하나도 없다 — 그 행렬은 배정도에서 특이행렬과 구별되지 않는다.

기하로 읽는다. $\kappa_2 = \sigma_{\max}/\sigma_{\min}$, 즉 $A$가 단위구를 타원체로 보낼 때 가장 긴 축과 가장 짧은 축의 비다. 어떤 방향은 크게 늘이고 어떤 방향은 거의 뭉개는 행렬이 조건이 나쁘다. 뭉개진 방향의 정보는 $\mathbf{b}$에 거의 남지 않으므로, $\mathbf{b}$의 작은 오차가 그 방향에서는 큰 $\delta\mathbf{x}$로 되살아난다. 대칭 행렬이면 특이값이 고유값의 절댓값이라 $\kappa_2 = \lvert \lambda \rvert_{\max} / \lvert \lambda \rvert_{\min}$이다.

그림으로 읽는다. 2차원에서 두 방정식은 두 직선이고 해는 교점이다. 두 직선이 거의 평행하면 한 직선을 조금만 옮겨도 교점이 멀리 달아난다 — 아래 그림.

그림 1.12-5. 2 × 2 선형계의 두 방정식을 직선으로 그린 것. 왼쪽은 두 직선이 거의 수직이라(조건이 좋다) 한 직선을 조금 옮겨도 교점이 조금만 움직인다. 오른쪽은 두 직선이 거의 평행이라(조건이 나쁘다) 같은 크기로 옮겼는데 교점이 멀리 달아난다. 두 경우 모두 행렬은 정칙이고 가우스 소거는 정상적으로 끝난다 — 조건수는 알고리즘이 아니라 문제의 성질이다.직접 작도

Kreyszig가 드는 표준 예가 힐베르트 행렬 $h_{ij} = 1/(i + j - 1)$이다. $3 \times 3$만 해도 $\kappa \approx 748$, $6 \times 6$은 $10^{7}$, $10 \times 10$은 $10^{13}$이다. 단항식 $1, x, x^2, \ldots$로 최소제곱 다항식 맞춤을 하면 정규방정식이 정확히 이 꼴이 되므로, 8절에서 말하듯 고차 다항식 맞춤은 단항식 기저로 하면 안 된다.

잔차는 오차가 아니다

계산한 $\tilde{\mathbf{x}}$가 얼마나 좋은지 보려고 잔차 $\mathbf{r} = \mathbf{b} - A\tilde{\mathbf{x}}$를 계산해 작으면 안심한다. 위험한 습관이다. $\mathbf{r} = A(\mathbf{x} - \tilde{\mathbf{x}})$이므로 위 정리에 $\delta\mathbf{b} = \mathbf{r}$을 넣으면

$$\frac{\lVert \mathbf{x} - \tilde{\mathbf{x}} \rVert}{\lVert \mathbf{x} \rVert} \le \kappa(A)\,\frac{\lVert \mathbf{r} \rVert}{\lVert \mathbf{b} \rVert}$$

잔차가 $10^{-10}$이어도 $\kappa = 10^{8}$이면 오차는 $10^{-2}$일 수 있다. 반대 방향의 부등식도 있어서, 잔차가 크면 오차도 확실히 크다. 즉 잔차는 오차의 하한은 주지만 상한은 $\kappa$를 알아야 준다. 2절의 뉴턴 수렴 판정에서 “잔차만으로는 안 된다”고 한 것이 이 식이다.

조건이 나쁜 문제를 만나면

조건수는 문제의 성질이라 알고리즘으로는 못 고친다. 할 수 있는 것은 문제를 바꾸는 것이다.

• 스케일링. 미지수와 방정식의 단위를 맞춘다. 행렬의 한 행이 $10^{20}$ 급이고 다른 행이 $10^{0}$ 급이면 그 자체로 $\kappa \ge 10^{20}$이다. 각 행·열을 적당한 대각행렬로 나누면 $\kappa$가 극적으로 준다.
• 기저 바꾸기. 단항식 대신 직교 다항식(르장드르·체비쇼프), 정규방정식 대신 QR(8절).
• 정칙화. 뭉개진 방향을 아예 버리거나($\sigma_{\min}$ 절단, 8절의 SVD) 작은 벌칙 항을 더해 그 방향을 고정한다($A^{\mathsf{T}}A + \mu I$, 14절의 레벤버그–마쿼트).

주의

소자 시뮬레이션이 스케일링 없이는 돌지 않는 이유. 드리프트–확산 방정식의 미지수는 $\psi \sim 1\ \mathrm{V}$, $n \sim 10^{20}\ \mathrm{cm^{-3}}$, 그리고 같은 격자에서 $p \sim 10^{0}\ \mathrm{cm^{-3}}$이다. 그대로 야코비 행렬을 세우면 조건수가 $10^{20}$을 넘어 배정도에서 특이행렬이다. 그래서 모든 시뮬레이터가 먼저 무차원화한다 — 길이는 디바이 길이, 퍼텐셜은 $V_T$, 농도는 진성 농도나 최대 도핑으로 나눈다(De Mari 스케일링). 그 위에 $n = n_i e^{(\psi - \phi_n)/V_T}$처럼 농도를 지수 안의 준페르미 퍼텐셜로 바꿔 미지수의 크기를 $V_T$ 단위로 맞추는 것(슬롯붐 변수)이 두 번째 층이다. 로그에 찍히는 잔차가 무차원인 것은 이 때문이고, 그것을 물리 단위로 되돌리려면 스케일 인자를 곱해야 한다.

최소제곱

방정식이 미지수보다 많을 때

측정은 언제나 미지수보다 많다. 저항 하나를 재려고 전류 열 점에서 전압을 재면 방정식이 열 개, 미지수는 하나다. $A\mathbf{x} = \mathbf{b}$에서 $A$가 $m \times n$, $m > n$이면 일반적으로 해가 없다 — $\mathbf{b}$가 $A$의 열공간 밖에 있다. 1.2에서 본 대로(1.2.4역행렬) 그때 할 수 있는 최선은 잔차의 길이를 최소로 하는 것이다.

$$\min_{\mathbf{x}} \lVert \mathbf{b} - A\mathbf{x} \rVert_2^2 = \min_{\mathbf{x}} \sum_{i=1}^{m}\Bigl(b_i - \sum_j a_{ij}x_j\Bigr)^2$$

제곱합을 쓰는 이유는 셋이다. 미분이 되어 닫힌 해가 나오고, 2-노름이라 기하(정사영)가 있으며, 잡음이 정규분포이면 이것이 최대가능도 추정이다(1.11.12회귀와 상관 (공정 데이터로의 적용은 19.9·19.12)). 각 $x_j$로 미분해 0으로 놓으면

$$\boxed{\;A^{\mathsf{T}}A\,\mathbf{x} = A^{\mathsf{T}}\mathbf{b}\;}$$

정규방정식이다. $n \times n$ 대칭 양의 정부호 계라 콜레스키(5절)로 풀린다. 기하로는 “잔차 $\mathbf{b} - A\mathbf{x}$가 $A$의 모든 열에 수직”이고, 그것이 정사영의 조건이다.

직선 맞춤

$y = a + bx$를 점 $(x_i, y_i)$ $m$개에 맞추면 $A$의 열이 $(1, \ldots, 1)$과 $(x_1, \ldots, x_m)$이므로 정규방정식은 $2 \times 2$다.

$$\begin{bmatrix} m & \sum x_i \\ \sum x_i & \sum x_i^2 \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} \sum y_i \\ \sum x_i y_i \end{bmatrix} \;\Longrightarrow\; b = \frac{m\sum x_i y_i - \sum x_i \sum y_i}{m \sum x_i^2 - (\sum x_i)^2}, \quad a = \bar y - b\bar x$$

$b$의 분모는 $m^2 \times (x\text{의 분산})$이고 분자는 $m^2 \times (x, y\text{의 공분산})$이다 — 기울기 = 공분산/분산. 다항식 $y = a_0 + a_1 x + \cdots + a_k x^k$이면 $A$의 열이 $x_i^0, \ldots, x_i^k$이고 정규방정식은 $(k+1) \times (k+1)$이다. 여기서 7절의 힐베르트 행렬이 나온다 — $x_i$가 $[0, 1]$에 고르게 있으면 $\sum x_i^{j+l} \approx m/(j + l + 1)$이라 정확히 힐베르트 꼴이다.

주의

반도체 추출에서 직선 맞춤이 숨어 있는 자리.

• TLM 접촉 저항. 간격 $d$가 다른 패드 사이의 저항 $R(d) = 2R_c + (R_{sh}/W)\,d$를 직선으로 맞춘다. 기울기가 면저항, 절편이 접촉 저항이다. 절편은 외삽이라 잡음에 예민하고, 그래서 $R_c$의 오차 막대가 $R_{sh}$보다 훨씬 크다.
• 아레니우스. $D = D_0 e^{-E_a/kT}$를 $\ln D$ 대 $1/T$로 그려 기울기에서 $E_a$를 읽는다. 로그를 취하는 순간 가중치가 바뀐다 — 작은 $D$의 상대 오차가 큰 $D$와 같은 무게로 들어간다. 원래 잡음이 절대량이면 로그 맞춤은 작은 값에 과도한 가중을 준다.
• $1/C^2$ 대 $V$. 4절에서 말한 대로 미분을 직선 맞춤으로 바꾼 것. 기울기에서 도핑, 절편에서 내부 전위.
• 문턱전압. $\sqrt{I_D}$(포화) 또는 $I_D/\sqrt{g_m}$을 $V_G$에 대해 직선으로 맞춰 $V_G$축 절편을 읽는다. 어느 구간을 맞추느냐로 값이 수십 mV 달라진다 — 모델이 직선인 구간이 좁기 때문이다.

정규방정식의 함정과 QR

정규방정식은 $\kappa(A^{\mathsf{T}}A) = \kappa(A)^2$이다(2-노름에서 $A^{\mathsf{T}}A$의 고유값이 $A$의 특이값의 제곱이므로). $\kappa(A) = 10^{6}$이면 정규방정식은 $10^{12}$ — 열 자리 얻을 문제를 네 자리로 만든다. $A^{\mathsf{T}}A$를 만드는 순간 자릿수를 절반 잃는다.

우회로는 $A$를 직접 직교 분해하는 것이다. 그람–슈미트(1.2.5벡터공간과 내적공간)나 하우스홀더 반사로 $A = QR$($Q$는 열이 정규직교인 $m \times n$, $R$은 $n \times n$ 상삼각)을 만들면

$$\lVert \mathbf{b} - A\mathbf{x} \rVert = \lVert \mathbf{b} - QR\mathbf{x} \rVert \;\Longrightarrow\; R\mathbf{x} = Q^{\mathsf{T}}\mathbf{b}$$

직교행렬은 길이를 보존하므로 조건수가 제곱되지 않고 $\kappa(A)$ 그대로다. 후진 대입 하나로 끝난다. 라이브러리의 lstsq가 기본으로 하는 것이 이것(또는 아래의 SVD)이고, 정규방정식은 교과서의 유도용이다.

특이값 분해 — 뭉개진 방향을 본다

1.2가 “1.12에서 다룬다”고 미뤄 둔 것이 이것이다. 임의의 $m \times n$ 행렬은

$$A = U\Sigma V^{\mathsf{T}}, \qquad \Sigma = \operatorname{diag}(\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_n \ge 0)$$

으로 쪼개진다. $U$($m \times m$)와 $V$($n \times n$)는 직교행렬, $\sigma_j$가 특이값이다. $A^{\mathsf{T}}A = V\Sigma^2 V^{\mathsf{T}}$이므로 $\sigma_j^2$은 $A^{\mathsf{T}}A$의 고유값, $V$의 열은 그 고유벡터다 — 1.3의 대칭 행렬 대각화를 정사각이 아닌 행렬로 확장한 것이다.

직관

SVD는 모든 행렬을 “회전 → 축별 늘이기 → 회전”으로 읽는다. $V^{\mathsf{T}}$가 입력을 돌려 특이 방향에 맞추고, $\Sigma$가 각 방향을 $\sigma_j$배 늘이고, $U$가 출력 공간에서 다시 돌린다. 단위구가 타원체로 가는 그림에서 반축 길이가 $\sigma_j$이고, 그래서 $\kappa_2 = \sigma_1/\sigma_n$이다. 1.3의 대각화가 “고유벡터 방향은 그대로”였다면, SVD는 “입력 방향과 출력 방향이 달라도 좋다”로 조건을 푼 대신 언제나 존재한다.

SVD가 있으면 최소제곱은 각 특이 방향에서 나눗셈이다.

$$\mathbf{x} = V\Sigma^{+}U^{\mathsf{T}}\mathbf{b} = \sum_{j} \frac{\mathbf{u}_j^{\mathsf{T}}\mathbf{b}}{\sigma_j}\,\mathbf{v}_j$$

$\Sigma^{+}$는 0이 아닌 $\sigma_j$를 역수로 바꾼 것이고, 이것이 1.2에서 이름만 나온 유사역행렬 $A^{+}$다. 식이 말해 주는 것이 두 가지다. 첫째, $\sigma_j$가 0이면 그 항을 버린다 — 계수가 모자란 경우에도 노름이 가장 작은 해를 준다. 둘째, $\sigma_j$가 0은 아닌데 아주 작으면 $1/\sigma_j$가 잡음을 폭발시킨다. 그때 그 항을 버리거나(절단 SVD) $\sigma_j/(\sigma_j^2 + \mu^2)$로 눌러 주는 것(티호노프 정칙화)이 7절 끝에서 말한 처방이다.

파라미터 추출에서 이것이 뜻하는 바 — 두 파라미터가 데이터에 같은 방향으로만 영향을 주면(예: 이동도와 채널 폭은 선형 영역 전류에 곱으로만 들어간다) 그 조합에 해당하는 $\sigma_j$가 0에 가깝고, 둘을 따로 추출하려는 시도는 잡음을 맞추는 일이 된다. SVD의 작은 특이값에 딸린 $\mathbf{v}_j$가 어느 파라미터 조합이 결정되지 않는지를 알려 준다.

가중 최소제곱

측정마다 잡음의 크기가 다르면 정확한 점에 더 무게를 준다.

$$\min \sum_i w_i\,(b_i - (A\mathbf{x})_i)^2, \qquad w_i = \frac{1}{\sigma_i^2}$$

잡음의 분산 $\sigma_i^2$의 역수가 통계적으로 맞는 가중치다. 앞의 아레니우스 예에서 로그를 취하면 자동으로 $w_i \propto D_i^2$이 되는 셈이고, 그것이 잡음의 성질과 맞는지가 로그 맞춤의 타당성을 정한다. 전류처럼 열 자릿수를 오가는 양을 맞출 때는 로그 잔차나 상대 잔차를 쓰지 않으면 큰 전류만 맞고 서브스레숄드 영역은 전혀 맞지 않는다.

비선형 최소제곱 — 가우스–뉴턴과 레벤버그–마쿼트

모델이 파라미터에 비선형이면($I = I_s(e^{V/nV_T} - 1)$에서 $n$, 타원계측의 굴절률 모델, SPICE 모델 파라미터 전부) 정규방정식이 없다. 2절의 뉴턴법을 끌어온다. 잔차 벡터 $\mathbf{r}(\mathbf{p}) = \mathbf{b} - \mathbf{f}(\mathbf{p})$를 현재 파라미터 근처에서 선형화하면 $\mathbf{r}(\mathbf{p} + \delta) \approx \mathbf{r} - J\delta$($J_{ij} = \partial f_i/\partial p_j$)이고, 이것의 최소제곱이

$$J^{\mathsf{T}}J\,\delta = J^{\mathsf{T}}\mathbf{r}$$

가우스–뉴턴법이다. 매 걸음이 선형 최소제곱이고, 해 근처에서 거의 2차로 수렴한다. 멀리서는 뉴턴법과 같은 이유로 튀므로 감쇠를 넣는다.

$$(J^{\mathsf{T}}J + \mu I)\,\delta = J^{\mathsf{T}}\mathbf{r}$$

$\mu$가 크면 $\delta \approx J^{\mathsf{T}}\mathbf{r}/\mu$ — 14절의 최급강하 방향으로 작은 걸음이고, $\mu \to 0$이면 가우스–뉴턴이다. 걸음이 잔차를 줄이면 $\mu$를 낮추고, 늘리면 $\mu$를 올린다. 이것이 레벤버그–마쿼트법이고, 곡선 맞춤 라이브러리의 기본 알고리즘이다. $\mu I$가 $J^{\mathsf{T}}J$의 작은 고유값을 밀어 올려 조건수를 낮추는 정칙화이기도 하다는 점이 7절과 잇는 고리다.

의문

맞춤이 “수렴했다”는 말과 “파라미터가 맞다”는 말은 다르다. 잔차가 작아도 파라미터가 틀릴 수 있는 경우가 둘이다.

① 파라미터가 서로 보상한다. 앞의 SVD 이야기다. $J^{\mathsf{T}}J$의 작은 고유값 방향으로는 파라미터를 얼마든지 움직여도 잔차가 거의 안 변한다. 그 방향의 파라미터 값은 초기값이 정한 것이지 데이터가 정한 것이 아니다. 신뢰구간(1.11.9신뢰구간)을 함께 보지 않으면 모른다. ② 모델이 틀렸다. 잔차가 무작위가 아니라 구조를 가지면(한 구간에서 계속 양수, 다른 구간에서 계속 음수) 파라미터가 아니라 모델을 바꿔야 한다. 잔차 그림을 보는 습관이 이래서 필요하다.

고유값의 수치 계산(멱승법·QR)

특성 다항식으로는 구하지 않는다

1.3에서 고유값을 $\det(A - \lambda I) = 0$의 근으로 정의했고, 손 계산은 그렇게 했다. 컴퓨터는 절대 그렇게 하지 않는다. 이유는 1.3의 1절에서 이미 봤다(1.3.1고유값 방정식과 특성 다항식 — 윌킨슨의 예) — 다항식의 근은 계수에 극도로 예민해서, 계수를 $10^{-10}$만큼 흔들면 근이 $10^{-1}$만큼 움직인다. 특성 다항식의 계수를 계산하는 과정에서 반올림이 반드시 들어가므로, 얻은 근은 원래 행렬의 고유값과 무관해진다.

그래서 수치 고유값 계산은 행렬에서 다항식을 거치지 않고 행렬을 직접 변형해 간다. 크게 두 부류다 — 고유값 하나(또는 몇 개)를 반복으로 얻는 멱승법 계열과, 전부를 닮음 변환으로 얻는 QR 계열.

게르슈고린 — 계산 없이 고유값이 있는 자리

먼저 계산 없이 아는 것이 있다. 모든 고유값은 복소평면의 원판

$$\lvert \lambda - a_{jj} \rvert \le \sum_{k \neq j} \lvert a_{jk} \rvert \qquad (j = 1, \ldots, n)$$

의 합집합 안에 있다(게르슈고린 정리). 증명은 한 줄이다 — 고유벡터의 가장 큰 성분을 $x_j$라 하고 $j$번째 방정식 $\sum_k a_{jk}x_k = \lambda x_j$를 $(\lambda - a_{jj})x_j = \sum_{k \neq j} a_{jk}x_k$로 정리한 뒤 $\lvert x_j \rvert$로 나누면 우변이 원판의 반지름 이하다.

대각 성분을 중심으로, 그 행의 비대각 절댓값 합을 반지름으로 하는 원판이다. 비대각이 작으면 고유값은 대각 성분 근처에 있고, 대각 우세면 원판이 원점을 포함하지 않아 행렬이 정칙이다 — 6절의 대각 우세 수렴 판정이 이 정리의 따름이다. 반복법의 초기값과 아래 QR의 이동량(shift)을 고를 때 쓴다.

멱승법 — 곱하고 또 곱한다

아무 벡터 $\mathbf{x}^{(0)}$에 $A$를 거듭 곱한다. 고유기저로 전개하면 (1.3.5닮음 변환과 대각화)

$$A^m\mathbf{x}^{(0)} = \sum_j c_j \lambda_j^m \mathbf{v}_j = \lambda_1^m\Bigl[c_1\mathbf{v}_1 + c_2\Bigl(\frac{\lambda_2}{\lambda_1}\Bigr)^m\mathbf{v}_2 + \cdots\Bigr]$$

$\lvert \lambda_1 \rvert > \lvert \lambda_2 \rvert \ge \cdots$이면 대괄호 안의 둘째 항부터 전부 0으로 가고, 절댓값이 가장 큰 고유값의 고유벡터만 남는다. 매번 노름으로 나눠 크기를 1로 유지하면 그것이 멱승법이다. 수렴 속도는 $\lvert \lambda_2/\lambda_1 \rvert^m$ — 두 고유값의 비가 수렴의 전부다. 비가 0.5면 빠르고, 0.99면 6절의 야코비 반복만큼 느리다.

고유값 자체는 레일리 몫으로 읽는다(1.3.7이차 형식과 주축 변환).

$$q = \frac{\mathbf{x}^{\mathsf{T}}A\mathbf{x}}{\mathbf{x}^{\mathsf{T}}\mathbf{x}}$$

$A$가 대칭이면 $q$는 고유값의 가중 평균이고, Kreyszig가 증명하는 오차 한계가 딸려 있다.

$$\lvert \lambda - q \rvert \le \delta, \qquad \delta^2 = \frac{\mathbf{x}^{\mathsf{T}}A^2\mathbf{x}}{\mathbf{x}^{\mathsf{T}}\mathbf{x}} - q^2 = \frac{\lVert A\mathbf{x} - q\mathbf{x} \rVert^2}{\lVert \mathbf{x} \rVert^2}$$

$\delta$는 “$\mathbf{x}$가 고유벡터에서 얼마나 벗어났는가”의 잔차이고, 어떤 고유값 하나는 반드시 $q$에서 $\delta$ 안에 있다. 고유벡터의 오차가 $\varepsilon$이면 레일리 몫의 오차는 $\varepsilon^2$이다 — 1.3의 변분 원리가 수치적으로 하는 말이다.

역반복 — 원하는 고유값을 골라 뽑는다

$A$ 대신 $(A - \sigma I)^{-1}$에 멱승법을 쓴다. 이 행렬의 고유값은 $1/(\lambda_j - \sigma)$이므로 $\sigma$에 가장 가까운 $\lambda_j$가 지배적이 되고, 수렴비는 $\lvert \lambda_j - \sigma \rvert / \lvert \lambda_k - \sigma \rvert$($\lambda_k$는 둘째로 가까운 것)다. $\sigma$를 목표 고유값 가까이 잡으면 이 비가 0에 가까워 두세 번에 끝난다. 역행렬은 구하지 않고, 매 반복 $(A - \sigma I)\mathbf{y} = \mathbf{x}$를 5절의 LU로 푼다 — 한 번 분해해 두고 우변만 바꾸는 바로 그 경우다.

$\sigma$를 매번 최신 레일리 몫으로 갱신하면 레일리 몫 반복이고, 대칭 행렬에서 수렴이 3차다 — 자릿수가 매번 세 배. 1.3의 7절이 “1.12에서 다룬다”고 미룬 것이 이것이다.

직관

$\sigma$가 고유값에 아주 가까우면 $A - \sigma I$는 거의 특이행렬이라 7절의 논리로 풀면 안 될 것 같다. 그런데 된다 — 오히려 그래서 된다. 조건이 나쁜 계는 오차를 “뭉개진 방향”, 즉 $\sigma$에 해당하는 고유벡터 방향으로 증폭시킨다. 우리가 원하는 것이 정확히 그 방향이다. 계산된 $\mathbf{y}$의 크기는 엉망이어도 방향은 정확하고, 정규화하면 크기는 버려진다. 불량 조건이 해가 되는 것은 답의 크기를 원할 때이고, 방향만 원하면 약이 된다.

QR 반복 — 전부 한꺼번에

고유값 전부가 필요하면 멱승법을 $n$번 돌리지 않는다. 행렬을 닮음 변환으로 삼각으로 만든다 (1.3.5닮음 변환과 대각화 — 닮음은 고유값을 보존한다). 절차는 두 단계다.

① 삼중대각(또는 헤센베르크)으로 줄인다. 하우스홀더 반사(직교행렬 $P = I - 2\mathbf{v}\mathbf{v}^{\mathsf{T}}$)를 $n - 2$번 적용해 대칭 행렬을 삼중대각 $T = P^{\mathsf{T}}AP$로 만든다. 비용 $O(n^3)$이 한 번. 비대칭이면 위 헤센베르크(대각 아래 한 줄까지만 0이 아님)까지 간다.

② QR 반복. $T$를 $T = QR$로 분해하고 순서를 바꿔 곱한다.

$$T_k = Q_k R_k, \qquad T_{k+1} = R_k Q_k = Q_k^{\mathsf{T}}\,T_k\,Q_k$$

$T_{k+1}$은 $T_k$와 닮음이므로 고유값이 같고, $k \to \infty$에서 **상삼각(대칭이면 대각)**으로 수렴한다. 대각선에 고유값이 절댓값 크기 순으로 늘어선다. 삼중대각이라 한 번의 QR 분해가 $O(n)$이고, 그래서 전체가 $O(n^2)$에 끝난다.

직관

QR 반복은 멱승법을 모든 방향에서 동시에 돌리는 것이다. $Q_1 Q_2 \cdots Q_k$의 첫 열은 멱승법의 $A^k\mathbf{e}_1$을 정규화한 것이고, 둘째 열은 첫 열에 직교한 채로 같은 일을 한 것, … $n$번째 열까지. 그래서 $(i+1, i)$ 성분의 수렴비가 $\lvert \lambda_{i+1}/\lambda_i \rvert$ — 멱승법의 수렴비가 인접한 고유값 쌍마다 하나씩이다.

인접한 고유값이 가까우면 느리다. 그래서 이동(shift)을 넣는다 — $T_k - \sigma I$를 분해하고 $R_kQ_k + \sigma I$로 되돌리면 수렴비가 $\lvert (\lambda_{i+1} - \sigma)/(\lambda_i - \sigma) \rvert$가 되고, $\sigma$를 마지막 대각 성분(윌킨슨 이동)으로 잡으면 역반복과 같은 이유로 맨 아래 고유값이 두세 번에 확정된다. 확정되면 그 행·열을 떼어 내고(수축, deflation) 크기를 하나 줄여 계속한다. 실제 라이브러리(LAPACK의 syev 류)가 이것이고, $n = 1000$이 1초 안에 끝난다.

큰 희소 행렬 — 몇 개만 필요할 때

밴드 계산이나 양자우물의 슈뢰딩거 방정식을 이산화하면 $n = 10^{4}$에서 $10^{7}$의 희소 행렬이 나오고, 필요한 것은 가장 낮은 몇 개의 고유값뿐이다(페르미 준위 근처의 밴드, 속박 준위 몇 개). 삼중대각화는 채워짐 때문에 못 쓰고, 행렬–벡터 곱만으로 가야 한다.

란초스법이 그 답이다. 멱승법이 마지막 벡터 하나만 남기고 버린 것을, 전부 보관해 정규직교화하면(6절의 크릴로프 부분공간) 그 부분공간에서 $A$의 표현이 작은 삼중대각 행렬이 되고, 그 고유값이 원래 행렬의 극단 고유값(가장 크고 가장 작은 것)에 빠르게 수렴한다. $m$번 곱하면 $m \times m$ 삼중대각 — 그것을 위의 QR로 푼다. 비대칭 버전이 아놀디법이고, 6절의 GMRES와 같은 뼈대다. 평면파 밀도범함수 계산이 매 자기일관 반복마다 하는 일이 이것이다 (정확히는 그 변형인 다비드슨법·켤레기울기 최소화).

주의

어느 고유값이 나왔는지 확인할 것. 멱승법은 절댓값 최대를 주고, 역반복은 $\sigma$에 가장 가까운 것을 주고, 란초스는 양 끝을 준다. 원하는 것이 “가장 작은 양의 고유값”이면 셋 중 어느 것도 곧바로 주지 않는다 — $\sigma = 0$의 역반복이 절댓값 최소를 주는데, 음의 고유값이 더 가까이 있으면 그것이 나온다. 양자우물의 바닥 상태를 구할 때 해밀토니안을 양의 정부호로 이동시켜 두는 습관이 이것 때문이다.

축퇴가 있으면 고유벡터는 정해지지 않는다(1.3.6축퇴와 고유기저). 멱승법·QR이 돌려주는 축퇴 고유벡터는 고유공간 안의 아무 기저이고, 실행마다 다를 수 있다. 대칭으로 분류하려면 1.3처럼 다른 연산자로 기저를 고정해야 한다.

ODE의 수치 해법(룽게-쿠타·다단계법)

초기값 문제를 걸어간다

$y' = f(x, y)$, $y(x_0) = y_0$. 1.4에서 해의 존재와 유일성을 보았고(1.4.3방향장과 해의 존재·유일성), 닫힌 해가 있는 경우도 몇 가지 보았다. 그 밖의 경우 — 이동도가 전기장에 비선형인 캐리어 운동, 비선형 소자가 든 회로의 과도 응답, 산화막 성장 속도식 — 는 시작점에서 출발해 한 걸음씩 전진한다. $x_n = x_0 + nh$에서의 근사 $y_n$을 앞의 값들로부터 만든다.

이 절의 이름표를 정리해 둔다. 일단계법(one-step)은 $y_{n+1}$을 $y_n$만으로 만들고 (오일러·룽게–쿠타), 다단계법(multistep)은 $y_n, y_{n-1}, \ldots$를 함께 쓴다(아담스). 명시적(explicit) 방법은 $y_{n+1}$이 우변에 나오지 않아 대입으로 끝나고, 암시적 (implicit) 방법은 $y_{n+1}$이 양변에 있어 매 걸음 방정식(비선형이면 2절의 뉴턴법)을 풀어야 한다. 암시적이 비싼데도 쓰는 이유가 이 절의 뒷부분이다.

오일러법 — 접선을 따라 한 걸음

가장 단순한 방법이다. 현재 점에서의 기울기 $f(x_n, y_n)$으로 $h$만큼 직진한다.

$$y_{n+1} = y_n + h\,f(x_n, y_n)$$

테일러 전개 $y(x_n + h) = y_n + hy' + \tfrac12 h^2 y'' + \cdots$의 첫 두 항이다. 버린 첫 항이 $\tfrac12 h^2 y''$이므로 한 걸음의 오차(국소 절단 오차)는 $O(h^2)$다. 그런데 구간 $[x_0, X]$를 건너려면 $(X - x_0)/h$ 걸음이 필요하고, 오차가 쌓이므로 끝점의 오차 (전역 오차)는 $O(h)$다. 이것을 1차 방법이라 한다 — 차수는 언제나 전역 오차로 센다.

1.4의 방향장(1.4.3방향장과 해의 존재·유일성) 위에서 오일러법은 화살표를 따라 짧은 직선을 잇는 일이다. 곡선이 휘는 쪽으로 항상 바깥에서 따라가므로, $y' = y$처럼 위로 볼록한 해에서는 늘 아래로 처지고, 그 처짐이 걸음마다 쌓인다.

개선된 오일러와 룽게–쿠타

처짐을 줄이려면 구간 안의 기울기를 여러 곳에서 재어 평균한다. 오일러법으로 끝점을 예측한 뒤 그곳의 기울기와 시작점의 기울기를 평균하는 것이 개선된 오일러법(호인법)이다.

$$\begin{aligned} k_1 &= h\,f(x_n, y_n) \\ k_2 &= h\,f(x_n + h,\ y_n + k_1) \\ y_{n+1} &= y_n + \tfrac12(k_1 + k_2) \end{aligned}$$

전역 오차 $O(h^2)$ — 4절의 사다리꼴 공식에 해당한다(실제로 $f$가 $x$만의 함수이면 정확히 사다리꼴이 된다). 기울기를 네 곳에서 재어 심프슨 가중치로 섞으면 고전적 4차 룽게–쿠타법이다.

$$\begin{aligned} k_1 &= h\,f(x_n,\ y_n) \\ k_2 &= h\,f(x_n + \tfrac12 h,\ y_n + \tfrac12 k_1) \\ k_3 &= h\,f(x_n + \tfrac12 h,\ y_n + \tfrac12 k_2) \\ k_4 &= h\,f(x_n + h,\ y_n + k_3) \\[4pt] y_{n+1} &= y_n + \tfrac16\,(k_1 + 2k_2 + 2k_3 + k_4) \end{aligned}$$

$k_1$은 시작점, $k_2$와 $k_3$은 중점(두 번째는 첫 번째 중점 기울기로 다시 예측한 중점), $k_4$는 끝점의 기울기이고, $1 : 2 : 2 : 1$의 가중치가 심프슨의 $1 : 4 : 1$이다. 전역 오차 $O(h^4)$. 한 걸음에 $f$를 네 번 계산하지만, 오일러법으로 같은 정확도를 얻으려면 걸음을 수천 배 잘게 해야 하므로 비교가 안 된다. RK4는 “특별한 사정이 없으면 이것”인 기본값이다.

주의

이 절의 $k_1, \ldots, k_4$는 기울기 × $h$, 즉 $y$의 증분이다. Kreyszig의 표기이고, 어떤 책은 $h$를 곱하지 않은 기울기 자체를 $k_i$로 쓰고 마지막 줄에 $h$를 곱한다. 둘 다 같은 방법이다. 11·12절에서 $k$는 시간 걸음 $\Delta t$로 뜻이 바뀐다.

걸음 폭 조절 — 오차를 재면서 걷는다

걸음 폭 $h$를 처음부터 고정하면 낭비이거나 위험하다. 해가 완만한 곳에서는 크게, 급한 곳에서는 작게 걸어야 한다. 그러려면 걸음마다 오차를 추정해야 하고, 4절의 리처드슨 요령이 그대로 쓰인다. 같은 걸음을 $h$로 한 번, $h/2$로 두 번 걸어 차이를 보면 오차가 추정된다. 더 싼 방법은 차수가 다른 두 공식이 같은 $k_i$를 공유하게 짜는 것이다 — 4차와 5차 공식이 여섯 개의 $k_i$를 나눠 쓰는 룽게–쿠타–펠베르크(RKF45)가 그것이고, 두 답의 차이가 곧 오차 추정이다. 추정이 허용치보다 크면 $h$를 줄여 다시 걷고, 훨씬 작으면 다음 걸음을 키운다.

라이브러리의 ode45, RK45가 이것이다. 사용자가 주는 것은 $h$가 아니라 허용 오차이고, $h$는 알고리즘이 정한다. SPICE의 시간 걸음도 마찬가지다 — .TRAN의 최대 걸음은 상한일 뿐, 실제 걸음은 국소 절단 오차 추정(LTE)이 정한다. “timestep too small” 오류는 오차 추정이 허용치 안으로 들어오지 않아 $h$를 줄이다 바닥에 닿았다는 뜻이다. 대개 회로 안의 불연속 (이상적 스위치, 꺾인 모델 방정식)이 원인이다.

다단계법 — 지난 걸음을 재활용한다

RK4는 한 걸음에 $f$를 네 번 계산하고 버린다. 다단계법은 지난 걸음에서 계산한 $f_n, f_{n-1}, f_{n-2}, \ldots$를 보관해 두고 그것으로 $f$를 다항식 보간(3절)한 뒤 적분(4절)한다. $y' = f$를 $x_n$에서 $x_{n+1}$까지 적분하면 $y_{n+1} - y_n = \int f\,dx$이고, 적분 안의 $f$를 지난 네 점의 3차 보간 다항식으로 바꾸면

$$y_{n+1} = y_n + \frac{h}{24}\,\bigl(55 f_n - 59 f_{n-1} + 37 f_{n-2} - 9 f_{n-3}\bigr)$$

4차 아담스–배시포스법이다. 계수 $55, -59, 37, -9$는 마디 밖(외삽)으로 적분하기 때문에 부호가 번갈아 크다 — 3절에서 경고한 외삽이다. 한 걸음에 $f$ 한 번으로 4차. 대신 처음 세 걸음은 RK4로 만들어야 하고(자기 시동이 안 된다), $h$를 바꾸려면 보관한 값을 다시 만들어야 한다.

외삽을 피하려면 $f_{n+1}$을 보간에 포함시킨다.

$$y_{n+1} = y_n + \frac{h}{24}\,\bigl(9 f_{n+1} + 19 f_n - 5 f_{n-1} + f_{n-2}\bigr)$$

아담스–몰턴법이고, 계수가 작고 얌전한 대신 $f_{n+1} = f(x_{n+1}, y_{n+1})$에 미지수가 들어 있어 암시적이다. 둘을 짝지어 쓴다 — 배시포스로 예측한 $y_{n+1}^{*}$를 몰턴의 우변에 넣어 수정한다(예측자–수정자). 두 값의 차이 $\lvert y_{n+1} - y_{n+1}^{*} \rvert$가 오차 추정이고, Kreyszig가 보이듯 실제 오차의 약 $1/15$가 그 차이다.

연립과 고계 — 벡터로 쓰면 같다

미지수가 여럿이거나 2계 이상이면 1.4에서 한 대로(1.4.9연립 ODE와 위상 평면) 1계 연립으로 바꾼다. $y'' = f(x, y, y')$는 $y_1 = y,\ y_2 = y'$로 놓아 $y_1' = y_2,\ y_2' = f(x, y_1, y_2)$이고, 위의 모든 공식은 $y,\ f,\ k_i$를 굵은 벡터로 바꾸면 그대로다. RK4의 $k_i$가 벡터가 되고 네 번의 $\mathbf{f}$ 계산이 각각 연립 전체를 평가한다.

강성 — 왜 암시적 방법이 필요한가

여기서부터가 이 절의 핵심이다. 시험 방정식 $y' = \lambda y$($\lambda < 0$, 해 $e^{\lambda x}$는 감쇠)에 오일러법을 쓰면

$$y_{n+1} = (1 + h\lambda)\,y_n = (1 + h\lambda)^{n+1}\,y_0$$

참해는 0으로 가는데 수치해가 0으로 가려면 $\lvert 1 + h\lambda \rvert < 1$, 즉 **$h < 2/\lvert \lambda \rvert$**여야 한다. 넘으면 부호를 바꾸며 폭발한다. 정확도가 아니라 안정성이 $h$의 상한을 정한 것이다.

연립에서는 $\lambda$가 야코비 행렬 $\partial\mathbf{f}/\partial\mathbf{y}$의 고유값들이고 (1.4.9연립 ODE와 위상 평면 — 선형화한 계의 모드), 가장 큰 $\lvert \lambda \rvert$가 $h$를 묶는다. 문제는 고유값들이 크게 다를 때 생긴다. 시정수 $1\ \mathrm{ps}$인 모드와 $1\ \mathrm{ms}$인 모드가 같은 계에 있으면, 빠른 모드는 몇 ps 만에 사라져 관심도 없는데 그것의 안정성 때문에 $h < 2\ \mathrm{ps}$로 묶여 $1\ \mathrm{ms}$까지 $10^{9}$ 걸음을 걸어야 한다. 이런 계를 강성(stiff)이라 하고, 강성비 $\lvert \lambda \rvert_{\max}/\lvert \lambda \rvert_{\min}$이 그 정도다.

암시적으로 바꾸면 사라진다. 후진 오일러법 $y_{n+1} = y_n + h f(x_{n+1}, y_{n+1})$을 시험 방정식에 쓰면

$$y_{n+1} = \frac{y_n}{1 - h\lambda} \;\Longrightarrow\; \Bigl\lvert \frac{1}{1 - h\lambda} \Bigr\rvert < 1 \quad\text{모든 } h > 0,\ \lambda < 0$$

어떤 $h$에서도 안정하다(A-안정). 빠른 모드는 $h$가 크면 한 걸음에 그냥 뭉개져 0이 되고 — 그것이 물리적으로 맞다, 어차피 사라질 모드였으니 — 느린 모드만 정확도에 맞는 $h$로 따라가면 된다. 걸음마다 $y_{n+1}$에 대한 방정식을 풀어야 하는 비용을 치르고 $10^{9}$ 걸음을 $10^{3}$ 걸음으로 바꾼다.

그림 1.12-6. 복소 hλ 평면에서 각 방법이 안정한 영역. 전진 오일러는 −1을 중심으로 한 반지름 1의 원판 안에서만, RK4는 그보다 조금 넓은 영역 안에서만 안정하다 — 둘 다 유계라서 |λ|가 큰 모드가 h를 묶는다. 후진 오일러는 +1을 중심으로 한 원판의 바깥 전부, 사다리꼴법은 왼쪽 반평면 전부가 안정 영역이라 감쇠하는 모드(Re λ < 0)에서는 h가 얼마든 안정하다. 강성 문제는 왼쪽 멀리(음의 실축 끝)에 고유값이 있는 문제이고, 그곳을 덮는 방법만 큰 h로 걸을 수 있다.직접 작도

직관

그림의 축이 $h\lambda$인 점을 보라. 방법의 안정 영역은 고정되어 있고, 문제의 고유값 $\lambda$에 $h$를 곱한 점이 그 안에 들어와야 한다. $h$를 줄이는 것은 모든 점을 원점 쪽으로 끌어당기는 일이라 언제나 결국 들어오지만, 가장 먼 점이 들어올 때까지 줄여야 한다. 안정 영역이 왼쪽 반평면 전체를 덮으면(A-안정) 감쇠 모드는 얼마나 멀리 있든 상관없다.

명시적 방법은 안정 영역이 언제나 유계다 — $y_{n+1}$이 $y_n$의 다항식이라 $\lvert h\lambda \rvert \to \infty$에서 반드시 발산하기 때문이다. 암시적 방법은 유리함수라 무한대에서 0이나 상수로 갈 수 있다. 강성은 명시적 방법으로는 원리적으로 못 넘는 벽이고, 이것이 SPICE와 소자 시뮬레이터가 전부 암시적인 이유다.

암시적 방법 중 실제로 쓰이는 것은 세 가지다.

방법	공식	차수	안정	비고
후진 오일러	$y_{n+1} = y_n + hf_{n+1}$	1	A-안정, 강하게 감쇠	빠른 모드를 과하게 죽인다(수치 감쇠)
사다리꼴	$y_{n+1} = y_n + \tfrac{h}{2}(f_n + f_{n+1})$	2	A-안정, 감쇠 없음	$h\lambda \to -\infty$에서 증폭률 $\to -1$: 부호가 번갈아 남는 링잉
기어(BDF)	$y_{n+1}$을 지난 값들의 후진 차분으로	1–6	2차까지 A-안정	SPICE의 METHOD=GEAR

사다리꼴법의 링잉이 SPICE 사용자가 만나는 “시간 걸음마다 부호가 뒤집히는 톱니”의 원인이다. 정확도는 좋지만 빠른 모드를 죽이지 않고 $\pm$로 튕기게 남긴다. 후진 오일러는 죽이지만 1차라 부정확하다. 그래서 실무는 사다리꼴로 걷다가 스위칭처럼 급한 순간에 후진 오일러 몇 걸음을 끼우거나, 2차 BDF(기어)를 쓴다.

의문

암시적 걸음마다 뉴턴법(2절)을 돌리는데, 그 안에서 선형계(5절)를 풀고, 그 선형계의 조건수 (7절)가 나쁘면 어떻게 되는가?

세 층이 그대로 겹친다. 회로 시뮬레이션의 한 시간 걸음은 ① 사다리꼴/BDF로 이산화 → ② 비선형 대수 방정식(KCL)을 뉴턴법으로 → ③ 걸음마다 희소 LU다. ②가 수렴하지 않으면 $h$를 줄인다 — $h$가 작을수록 $y_{n+1}$이 $y_n$에 가까워 뉴턴의 초기값이 좋아지기 때문이다. 그래서 “timestep too small”은 안정성·정확도만이 아니라 뉴턴이 수렴하지 않아서일 때가 많고, 그때 GMIN 스테핑이나 소스 스테핑(전원을 천천히 올려 연속적으로 해를 따라간다)이 처방이다. 연속법(continuation)이라 부르는 이 요령은 2절의 “좋은 초기값을 만드는” 기술이다.

타원형 PDE의 차분법과 경계 조건

격자 위에 방정식을 놓는다

1.7에서 2계 PDE를 타원·포물·쌍곡으로 나누었다(1.7.32계 PDE의 분류 — 타원·포물·쌍곡). 타원형은 경계값 문제다 — 시간이 없고, 영역 전체의 해가 경계 조건으로 한꺼번에 정해진다. 라플라스 $\nabla^2 u = 0$과 푸아송 $\nabla^2 u = f$가 대표이고(1.7.10라플라스·푸아송 방정식), 반도체에서는 정전 퍼텐셜이 곧 이것이다.

영역을 간격 $h$의 격자로 덮고 격자점 $(x_i, y_j) = (ih, jh)$에서의 값을 $u_{ij}$로 쓴다. 4절의 2계 중심 차분을 $x$와 $y$에 각각 쓰면

$$\nabla^2 u \approx \frac{u_{i+1,j} - 2u_{ij} + u_{i-1,j}}{h^2} + \frac{u_{i,j+1} - 2u_{ij} + u_{i,j-1}}{h^2}$$

이므로 푸아송 방정식은 격자점마다 하나씩

$$\boxed{\; u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4u_{ij} = h^2 f_{ij} \;}$$

5점 공식이다. 절단 오차는 $O(h^2)$. 라플라스($f = 0$)이면 $u_{ij} = \tfrac14(u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1})$ — 각 점은 네 이웃의 평균이다. 이것은 조화함수의 평균값 성질(1.1.9퍼텐셜 이론)의 이산판이고, Sadiku가 정전기장을 격자에서 손으로 풀 때 쓰는 규칙 그 자체다.

그림 1.12-7. 5점 공식의 스텐실과 세 종류의 경계. 안쪽 격자점(가운데)은 상하좌우 네 이웃과 자기 자신으로 식 하나를 만든다. 왼쪽 경계는 디리클레 조건 — 경계 위의 값이 주어져 있어 그 점은 미지수가 아니라 우변으로 넘어간다. 위쪽 경계는 노이만 조건 — 경계 바깥에 가상점을 하나 두고, 법선 도함수의 중심 차분이 주어진 값이 되게 가상점을 안쪽 점으로 표현해 소거한다. 오른쪽 아래처럼 경계가 격자와 어긋나면 그 점에서만 간격이 다른 비균등 공식을 쓴다.직접 작도

행렬로 모으면

$N \times N$개의 안쪽 격자점에 식이 $N^2$개 — $n = N^2$차 선형계 $A\mathbf{u} = \mathbf{b}$다. 격자점을 행 순서로 번호 매기면 $A$는 대각이 $-4$, 좌우 이웃이 $\pm1$ 자리, 상하 이웃이 $\pm N$ 자리에 1이 있는 블록 삼중대각 행렬이다.

$$A = \begin{bmatrix} T & I & & \\ I & T & I & \\ & \ddots & \ddots & \ddots \\ & & I & T \end{bmatrix}, \qquad T = \begin{bmatrix} -4 & 1 & & \\ 1 & -4 & 1 & \\ & \ddots & \ddots & \ddots \\ & & 1 & -4 \end{bmatrix}$$

성질 셋이 앞 절들을 불러온다. 대칭이고, $-A$가 양의 정부호이며(고유값이 6절에서 본 $-4 + 2\cos\frac{p\pi}{N+1} + 2\cos\frac{q\pi}{N+1}$, 전부 음수), 띠 폭이 $N$인 희소 행렬이다. 그래서 콜레스키(5절, 비용 $N^4$), 가우스–자이델(6절 — 타원형 문제에서는 리브만 반복이라 부른다), SOR, 그리고 켤레 기울기·다중 격자가 전부 후보다. 조건수는 $\kappa \approx 4N^2/\pi^2$ — 격자를 촘촘히 할수록 나빠진다(7절).

직관

$\kappa \propto N^2$을 물리로 읽으면, 이산 라플라시안의 가장 큰 고유값은 이웃끼리 부호가 반대인 체스판 모드($\approx -8/h^2$), 가장 작은 고유값은 영역 전체에 걸친 반파장 모드 ($\approx -2\pi^2/L^2$)다. 비는 $(L/h)^2 = N^2$. 즉 격자가 담을 수 있는 가장 짧은 파장과 가장 긴 파장의 비가 조건수다. 10절의 강성비와 같은 구조이고, 실제로 12절에서 이 행렬을 시간 방향으로 적분하면 그 고유값들이 그대로 강성을 만든다.

경계 조건 — 세 종류, 세 처리

디리클레(값이 주어짐, $u = g$). 경계 격자점의 값을 알고 있으므로 그 점은 미지수가 아니다. 안쪽 점의 식에 경계값이 나오면 우변으로 넘긴다. 가장 쉽다. 반도체에서는 옴 접촉 — 접촉 전극의 퍼텐셜은 인가 전압으로 고정된다.

노이만(법선 도함수가 주어짐, $\partial u/\partial n = g$). 경계 점의 값 자체가 미지수다. 경계 바깥에 가상점(ghost point) $u_{i,N+1}$을 두고, 법선 도함수의 중심 차분

$$\frac{u_{i,N+1} - u_{i,N-1}}{2h} = g_i \;\Longrightarrow\; u_{i,N+1} = u_{i,N-1} + 2hg_i$$

로 가상점을 안쪽 점으로 바꿔 5점 공식에 넣는다. 경계 점에도 식이 하나 생기고, 그 식은 $2u_{i,N-1}$ 항을 갖는다. 전진 차분 $(u_{i,N} - u_{i,N-1})/h = g_i$로 해도 되지만 $O(h)$라 전체의 $O(h^2)$를 망친다 — Kreyszig가 가상점을 쓰는 이유다. 반도체에서는 대칭면 (소자의 절반만 풀 때 $\partial\psi/\partial n = 0$)과 인공 경계(멀리 떨어진 기판 바닥, 전기장이 0이라 놓는다)가 노이만이다.

혼합·로빈($\alpha u + \beta\,\partial u/\partial n = g$). 위 둘의 조합으로, 같은 가상점 요령이 든다. 반도체에서는 쇼트키 접촉(전류가 표면 재결합 속도에 비례)이나 대류 열전달 경계가 이것이다.

주의

순수 노이만 문제는 해가 하나로 정해지지 않는다. 모든 경계에서 도함수만 주면 $u + C$도 해라서 $A$가 특이행렬이다(상수 벡터가 핵). 물리적으로는 “퍼텐셜의 영점을 어디에 둘지 말하지 않았다”는 뜻이고, 열 문제로는 “경계로 드나드는 열의 합이 0이 아니면 정상 상태가 없다”는 뜻이다. 한 점이라도 디리클레로 고정하거나, 특이한 계를 그대로 두고 반복법을 써서 해 하나를 고른다. 소자 시뮬레이션이 접촉 하나 없는 구조를 거부하는 이유가 이것이다.

노이만 경계에서 5점 공식은 대각 우세의 등호가 깨지지 않는다 — 가상점 소거로 비대각 계수가 2가 되어도 대각 4와 같으므로 6절의 수렴 판정이 유지된다. 그러나 대칭성은 깨지므로, 대칭을 원하면 경계 식을 $\tfrac12$로 나눠 맞춘다.

비균등 격자와 계면

경계가 곡선이면 격자와 어긋난다. 어긋난 점에서는 이웃까지의 거리가 $h$가 아니라 $ah,\ bh$ ($a, b < 1$)이고, Kreyszig의 비균등 공식

$$\frac{2}{h^2}\left[\frac{u_{i+1,j}}{a(1+a)} + \frac{u_{i-1,j}}{1+a} - \frac{u_{ij}}{a}\right] \quad(x\text{ 방향, 오른쪽 이웃까지 } ah)$$

를 그 점에서만 쓴다. 절단 오차가 그 점에서 $O(h)$로 떨어지지만 몇 점뿐이라 전체에는 크게 영향하지 않는다.

더 중요한 것은 계면이다. 소자에는 산화막/실리콘 경계가 있고, 그곳에서 유전율이 $3.9\varepsilon_0$에서 $11.7\varepsilon_0$으로 뛴다. 푸아송 방정식은 $\nabla\cdot(\varepsilon\nabla\psi) = -\rho$이고, 계면에서 성립해야 하는 것은 $\psi$의 연속과 전속 밀도 $\varepsilon\,\partial\psi/\partial n$의 연속이다. 5점 공식을 그냥 쓰면 두 번째가 깨진다. 표준 처방은 유한 체적(box) 이산화다 — 각 격자점 주위에 상자를 두고 발산 정리(1.1.7발산정리와 스토크스 정리)로 방정식을 “상자 표면을 지나는 전속의 합 = 상자 안의 전하”로 바꾼다. 상자 면마다 그 면이 놓인 재료의 $\varepsilon$을 쓰면 계면 조건이 저절로 들어가고, 전하 보존이 격자 수준에서 정확히 성립한다. 균일한 재료에서는 5점 공식과 같은 식이 나오므로 손해가 없다. 소자 시뮬레이터가 전부 이 방식이고, 6·7절에서 “이산 라플라시안”이라 부른 행렬은 정확히는 이 상자 행렬이다.

직관

격자 간격은 물리가 정한다. 퍼텐셜이 변하는 가장 짧은 길이 척도를 격자가 담아야 한다. 반도체에서 그것은 디바이 길이 $L_D = \sqrt{\varepsilon_s V_T / qN}$이고, 도핑 $10^{18}\ \mathrm{cm^{-3}}$에서 약 4 nm다. 공핍층 가장자리, 반전층, 접합 근처는 $h \lesssim L_D$로 잘게, 중성 벌크는 성기게 — 그래서 소자 시뮬레이션의 격자는 언제나 비균등이고, 격자 생성이 시뮬레이션의 절반이다. $h$를 $L_D$보다 크게 잡으면 5점 공식이 틀리는 것이 아니라, 지수적으로 변하는 전하 밀도를 격자가 보지 못해 답이 통째로 어긋난다. 12절의 안정성 조건이 시간 걸음을 묶듯, 디바이 길이가 공간 걸음을 묶는다.

정확도의 확인 — 격자를 절반으로

수치해가 맞는지는 참해를 모르니 직접 확인할 수 없다. 대신 격자를 절반으로 줄여 다시 푼다. 오차가 $Ch^2$이면 두 해의 차이가 오차의 $3/4$이고(4절의 리처드슨 논리), 차이가 $1/4$로 줄면 방법이 2차로 수렴하고 있다는 증거다. 차이가 줄지 않으면 격자가 아직 물리를 담지 못한 것이고(위의 디바이 길이), $1/2$로만 줄면 어딘가에 $O(h)$짜리 처리(전진 차분 경계, 비균등 점)가 지배하고 있는 것이다. 이 검사를 하지 않은 시뮬레이션 결과는 그림이지 답이 아니다.

포물형·쌍곡형 PDE와 안정성

시간이 들어오면 — 공간은 차분, 시간은 전진

열·확산 방정식 $u_t = c^2 u_{xx}$(1.7.6열·확산 방정식)와 파동 방정식 $u_{tt} = c^2 u_{xx}$(1.7.5파동 방정식과 달랑베르 해)는 초기값 문제다. 공간은 11절처럼 격자로 놓고 ($x_i = ih$), 시간은 10절처럼 걸어간다($t_j = jk$). 이 절에서 $k$는 시간 걸음이다. $u_{ij} = u(x_i, t_j)$이고, 첫 인덱스가 공간, 둘째가 시간이다.

공간을 먼저 이산화하면 격자점마다 ODE 하나가 생긴다 — $du_i/dt = c^2(u_{i+1} - 2u_i + u_{i-1})/h^2$. 이것은 $n$차 연립 ODE이고(선분법, method of lines), 그 야코비 행렬은 11절의 이산 라플라시안이다. 그러므로 10절의 모든 결론이 적용된다 — 고유값의 범위가 $h$에 따라 $-4c^2/h^2$까지 뻗으니 강성이고, 명시적 방법은 $k$가 묶인다. 그것이 아래의 안정성 조건이다.

명시적 방법과 $r \le 1/2$

시간에 전진 차분, 공간에 중심 차분(FTCS)을 쓰면

$$\frac{u_{i,j+1} - u_{ij}}{k} = c^2\,\frac{u_{i+1,j} - 2u_{ij} + u_{i-1,j}}{h^2}$$ $$\Longrightarrow\quad u_{i,j+1} = (1 - 2r)\,u_{ij} + r\,(u_{i+1,j} + u_{i-1,j}), \qquad r = \frac{c^2 k}{h^2}$$

새 값이 지난 시각의 세 값으로 대입 한 번에 나온다. 격자비 $r$이 방법의 전부다. $r = 1/2$이면 $u_{i,j+1} = \tfrac12(u_{i+1,j} + u_{i-1,j})$ — 자기 값은 사라지고 양쪽 이웃의 평균이 된다. $r > 1/2$이면 $u_{ij}$의 계수가 음수가 되고, 그때부터 문제가 시작된다.

폰 노이만 안정성 해석. 오차를 푸리에 모드로 쪼개(1.10.6이산 푸리에 변환과 FFT) $u_{ij} = G^j e^{\mathrm{i}\theta i}$($\theta$는 격자 파수, 여기서 $\mathrm{i}$는 허수단위)를 넣으면 한 시간 걸음의 증폭률이

$$G(\theta) = 1 - 2r + r\,(e^{\mathrm{i}\theta} + e^{-\mathrm{i}\theta}) = 1 - 4r\sin^2\frac{\theta}{2}$$

모든 $\theta$에서 $\lvert G \rvert \le 1$이어야 오차가 자라지 않는다. $G$의 최솟값은 $\theta = \pi$(이웃끼리 부호 반대인 체스판 모드)에서 $1 - 4r$이고, 이것이 $-1$ 이상이려면

$$\boxed{\;r = \frac{c^2 k}{h^2} \le \frac12\;}$$

넘으면 체스판 모드가 부호를 바꾸며 지수적으로 자란다 — 처음엔 안 보이다가 몇 십 걸음 뒤에 톱니가 폭발하는 전형적 불안정이다.

직관

물리로 읽으면 한 시간 걸음 동안 확산이 퍼지는 거리 $\sqrt{2c^2 k}$가 격자 간격 $h$를 넘으면 안 된다는 뜻이다. 명시적 공식은 한 걸음에 이웃 한 칸까지만 정보를 보내는데, 물리는 그보다 멀리 보내려 하니 격자가 따라가지 못한다.

대가가 크다. $k \le h^2/2c^2$이므로 공간을 두 배 촘촘히 하면 시간 걸음은 네 배 줄어야 하고, 계산량은 여덟 배다. 확산 시뮬레이션에서 $h$를 나노미터로 잡으면 $k$가 피코초급이 되어 초 단위의 열처리를 명시적으로 걷는 것은 불가능하다. 10절의 강성 이야기가 그대로다 — 야코비 행렬의 최대 고유값 $4c^2/h^2$에 오일러의 $k < 2/\lvert\lambda\rvert$를 넣으면 정확히 $r \le 1/2$이 나온다.

크랭크–니콜슨 — 암시적으로

공간 차분을 지난 시각과 새 시각의 평균으로 잡으면(10절의 사다리꼴법에 해당)

$$(2 + 2r)\,u_{i,j+1} - r\,(u_{i+1,j+1} + u_{i-1,j+1}) = (2 - 2r)\,u_{ij} + r\,(u_{i+1,j} + u_{i-1,j})$$

크랭크–니콜슨법이다. 좌변에 새 시각의 값 세 개가 얽혀 있으므로 매 걸음 삼중대각 선형계를 풀어야 한다 — 5절의 토머스 알고리즘으로 $O(n)$이니 명시적 방법과 비용이 거의 같다. 증폭률은

$$G(\theta) = \frac{1 - 2r\sin^2(\theta/2)}{1 + 2r\sin^2(\theta/2)} \;\Longrightarrow\; \lvert G \rvert \le 1 \quad\text{모든 } r$$

무조건 안정이고, 시간·공간 모두 2차 정확도 $O(k^2 + h^2)$다. Kreyszig는 $r = 1$ ($k = h^2/c^2$)을 권하는데, 그때 좌변이 $4u_{i,j+1} - u_{i+1,j+1} - u_{i-1,j+1}$로 특히 간단해진다.

주의

무조건 안정이 무조건 정확은 아니다. $r$이 크면 $\theta = \pi$에서 $G \to -1$이라 체스판 모드가 사라지지 않고 부호만 바꾸며 남는다. 초기 조건에 불연속(계단형 도핑 프로파일)이 있으면 그 자리에서 진동이 오래 남는 것이 이것이다. 10절의 사다리꼴 링잉과 같은 현상이고, 처방도 같다 — 처음 몇 걸음을 후진 오일러($G = 1/(1 + 4r\sin^2(\theta/2))$, 언제나 양수이고 큰 $r$에서 0으로 감)로 걸어 불연속을 누그러뜨린 뒤 크랭크–니콜슨으로 넘어간다.

파동 방정식과 CFL 조건

$u_{tt} = c^2 u_{xx}$는 시간에도 2계 중심 차분을 쓴다.

$$\frac{u_{i,j+1} - 2u_{ij} + u_{i,j-1}}{k^2} = c^2\,\frac{u_{i+1,j} - 2u_{ij} + u_{i-1,j}}{h^2}$$ $$\Longrightarrow\quad u_{i,j+1} = 2(1 - r^{*})\,u_{ij} + r^{*}\,(u_{i+1,j} + u_{i-1,j}) - u_{i,j-1}, \qquad r^{*} = \frac{c^2 k^2}{h^2}$$

새 값이 두 시각 전의 값까지 쓰므로 첫 걸음을 따로 만들어야 한다. 초기 속도 $u_t(x, 0) = g(x)$를 중심 차분 $(u_{i,1} - u_{i,-1})/2k = g_i$로 쓰고 가상 시각 $u_{i,-1}$을 소거하면

$$u_{i,1} = (1 - r^{*})\,u_{i,0} + \tfrac12 r^{*}\,(u_{i+1,0} + u_{i-1,0}) + k\,g_i$$

폰 노이만 해석은 $G$에 대한 2차 방정식이 되고, 결과는

$$\boxed{\;r^{*} = \frac{c^2 k^2}{h^2} \le 1 \quad\Longleftrightarrow\quad c\,k \le h\;}$$

쿠랑–프리드리히스–레비(CFL) 조건이다. 한 시간 걸음 동안 파동이 진행하는 거리 $ck$가 격자 간격 $h$를 넘으면 안 된다.

그림 1.12-8. CFL 조건의 기하. 파동 방정식에서 점 (xᵢ, tⱼ₊₁)의 참값은 아래로 뻗은 특성선 x ± ct 안쪽(물리적 의존 영역)의 초기값으로 정해진다. 차분 공식은 그 점을 아래 시각의 이웃 세 점으로만 만들므로, 수치적 의존 영역은 기울기 h/k의 삼각형이다. 왼쪽처럼 ck ≤ h이면 수치 삼각형이 물리 삼각형을 품어 필요한 정보가 전부 들어온다. 오른쪽처럼 ck > h이면 물리 삼각형이 더 넓어, 답에 영향을 줘야 할 초기값이 계산에 들어오지 않는다 — 어떤 공식으로도 맞을 수 없고, 실제로는 폭발한다.직접 작도

직관

CFL 조건은 안정성 계산이기 전에 인과율이다. 1.7의 달랑베르 해에서 $(x, t)$의 값은 $[x - ct, x + ct]$의 초기값으로 정해진다(의존 영역). 차분 공식이 $j$ 걸음 뒤에 볼 수 있는 초기값의 범위는 $x \pm jh$, 즉 $x \pm (h/k)t$다. 물리가 정보를 $c$의 속도로 보내는데 격자가 $h/k$의 속도로만 받으면, $h/k < c$일 때 격자는 필요한 정보를 받지 못한 채 답을 만든다. 그것은 정확도의 문제가 아니라 원리적으로 틀린 답이고, 안정성 해석은 그 틀림이 폭발로 나타난다는 것을 확인할 뿐이다. 이것이 필요조건인 이유이며, 어떤 명시적 방법도 이 조건을 피하지 못한다.

$r^{*} = 1$, 즉 $k = h/c$이면 차분 공식이 달랑베르 해를 정확히 재현한다 — 격자가 특성선 위에 정확히 놓이기 때문이다. Kreyszig가 예제에서 $k = h$($c = 1$)로 잡는 이유이고, 이때 공식은 $u_{i,j+1} = u_{i+1,j} + u_{i-1,j} - u_{i,j-1}$로 줄어든다. 파동은 확산과 달리 “정확히 격자 속도로” 걷는 것이 최선이고, 더 잘게 걷는다고 좋아지지 않는다(수치 분산이 생긴다).

세 방법을 한 줄로 — 라크스의 동치 정리

이 절이 세운 조건들을 하나로 묶는 정리가 있다. 차분 공식이 일관(consistent: $h, k \to 0$에서 원래 PDE로 돌아간다)하고 안정(오차가 유계)하면, 그리고 그때에만, 수치해가 참해로 수렴한다. 일관성은 테일러 전개로 확인하고(절단 오차가 0으로 가는가), 안정성은 폰 노이만 해석으로 확인한다. 둘 다 국소적 계산인데 그 결과가 전역적 수렴을 보장한다 — 4절에서 “오차의 차수를 알면 그 차수로 수렴한다”고 한 것의 PDE 버전이다.

이류항 — 풍상 차분과 셰르페터–검멜

반도체의 연속 방정식은 순수 확산이 아니다.

$$\frac{\partial n}{\partial t} = \frac{1}{q}\nabla\cdot J_n + G - R, \qquad J_n = q\mu_n n E + qD_n \nabla n$$

드리프트 항 $\mu_n nE$는 **1계 미분(이류)**이고, 이류–확산 방정식은 이 절의 어느 쪽도 아닌 중간이다. 이류항을 중심 차분으로 이산화하면 격자 페클레 수 $\mathrm{Pe} = \mu E h / 2D = Eh / 2V_T$가 1을 넘을 때 — 즉 한 격자 간격의 전압 강하가 $2V_T \approx 52\ \mathrm{mV}$를 넘을 때 — 농도가 이웃끼리 부호를 바꾸며 진동한다(음의 농도!). 공핍층의 전기장 $10^{5}\ \mathrm{V/cm}$에서 이것은 $h > 5\ \mathrm{nm}$면 벌써 깨진다는 뜻이다.

처방은 풍상 차분(upwind) — 이류 방향의 상류 값만 쓰는 것 — 인데 $O(h)$이고 수치 확산이 크다. 반도체는 더 좋은 답을 갖고 있다. 두 격자점 사이에서 전류 $J$가 일정하고 $E$가 상수라 가정하면 $J = q\mu nE + qD\,dn/dx$는 $n$에 대한 1계 선형 ODE라 정확히 풀리고, 그 해가

$$J_{i+1/2} = \frac{qD_n}{h}\Bigl[n_{i+1}\,B\Bigl(\frac{\psi_{i+1} - \psi_i}{V_T}\Bigr) - n_i\,B\Bigl(\frac{\psi_i - \psi_{i+1}}{V_T}\Bigr)\Bigr], \qquad B(x) = \frac{x}{e^x - 1}$$

셰르페터–검멜 이산화다. $B$는 베르누이 함수이고, $\lvert\psi_{i+1} - \psi_i\rvert \ll V_T$이면 $B \approx 1 - x/2$가 되어 중심 차분으로, $\gg V_T$이면 한쪽 $B$가 0으로 가서 풍상 차분으로 저절로 넘어간다. 1절에서 예고한 대로 드리프트와 확산을 하나의 지수 식으로 합쳐 상쇄를 피했고, 어떤 $h$에서도 농도가 양수로 유지된다. 소자 시뮬레이터가 예외 없이 이것을 쓰는 이유이며, 격자 페클레 수의 제약이 사라져 공핍층을 굳이 $V_T/E$보다 잘게 자르지 않아도 된다.

의문

안정하다는 것과 정확하다는 것과 물리적이라는 것은 서로 다른 세 조건인가?

그렇다. 크랭크–니콜슨은 안정하고 2차 정확하지만 큰 $r$에서 음의 농도를 만들 수 있다 (체스판 모드의 $G \to -1$). 후진 오일러는 안정하고 양수를 보존하지만 1차다. 명시적 FTCS는 $r \le 1/2$에서 셋 다 만족하지만 그 $r$이 감당이 안 된다. 셰르페터–검멜은 정확도를 조금 양보하고(균일한 $E$의 가정) 양수 보존을 얻은 것이다. 시뮬레이터가 방법을 고르는 기준은 대개 세 번째 — 음의 농도, 음의 전하, 보존되지 않는 전하는 어떤 정확도로도 용서되지 않는다.

유한요소법과 모멘트법 개요

차분법이 불편한 자리

11·12절의 차분법은 격자가 직교해야 한다. 곡선 경계, 비스듬한 계면, 국소적으로만 촘촘해야 하는 영역에서는 직교 격자가 낭비이거나 부정확하다. 그리고 미분을 차분으로 바꾸는 방식이라 해가 매끄럽지 않은 곳(재료 경계의 꺾임)에서 근거가 약하다. 두 방법이 이 불편을 다른 길로 피한다 — 유한요소법(FEM)은 미분 방정식을 적분 형태(변분 원리)로 바꾸고, 모멘트법 (MoM)은 아예 적분 방정식에서 출발한다. Sadiku의 14장이 정전기 문제로 둘을 소개하는 순서를 그대로 따른다.

유한요소법 — 에너지를 최소화하는 조각 함수

① 출발점은 미분 방정식이 아니라 에너지다. 라플라스 방정식 $\nabla^2 V = 0$의 해는 경계 조건을 만족하는 모든 함수 중 정전 에너지

$$W = \frac{\varepsilon}{2}\int_{\Omega} \lvert \nabla V \rvert^2\, dS$$

를 최소로 하는 함수다(1.6의 변분 원리 — 1.6.6변분법과 오일러 방정식. 오일러 방정식을 쓰면 $\nabla^2 V = 0$이 나온다). 그러니 방정식을 푸는 대신 $W$를 최소화하면 된다. 1.6의 레일리–리츠법(1.6.8레일리-리츠 변분 근사)이 정확히 이것인데, 거기서는 시험 함수를 영역 전체에 걸친 몇 개의 함수로 잡았다. 유한요소법의 새 발상은 시험 함수를 작은 조각마다 따로 잡는 것이다.

② 영역을 요소로 쪼갠다. 2차원이면 삼각형, 3차원이면 사면체. 크기가 제각각이어도 되고, 경계 모양을 따라가도 되며, 필요한 곳만 잘게 해도 된다. 꼭짓점을 절점(node)이라 하고, 미지수는 절점에서의 값 $V_1, \ldots, V_n$이다.

③ 요소 안에서는 1차식. 삼각형 요소 $e$ 안에서 $V_e(x, y) = a + bx + cy$로 잡는다. 세 꼭짓점의 값 $V_{e1}, V_{e2}, V_{e3}$이 $a, b, c$를 정하므로

$$V_e(x, y) = \sum_{i=1}^{3} \alpha_i(x, y)\, V_{ei}, \qquad \alpha_i(\text{꼭짓점 } j) = \delta_{ij}$$

$\alpha_i$가 형상 함수(shape function)다 — 3절의 라그랑주 기저 $L_k$의 2차원판으로, 자기 꼭짓점에서 1, 다른 꼭짓점에서 0인 평면이다. Sadiku의 표기로 $\alpha_1 = \frac{1}{2A}[(x_2 y_3 - x_3 y_2) + (y_2 - y_3)x + (x_3 - x_2)y]$이고($A$는 삼각형 넓이), 나머지는 첨자를 순환한 것이다.

그림 1.12-9. 1차원에서 본 유한요소 기저. 위: 절점마다 자기 절점에서 1이고 이웃 절점에서 0으로 떨어지는 모자 함수(hat function)가 하나씩 있다. 아래: 절점값 V_j에 모자 함수를 곱해 더하면 절점을 지나는 꺾은선이 되고, 그것이 유한요소 근사다. 절점값 하나를 바꾸면 그 절점에 붙은 두 요소만 바뀐다 — 강성 행렬이 희소한 이유이자 3절의 스플라인과 같은 국소성이다. 2차원의 형상 함수 αᵢ는 이 모자를 삼각형 위의 피라미드로 세운 것이다.직접 작도

④ 요소의 에너지를 절점값으로 쓴다. $\nabla V_e = \sum_i V_{ei}\nabla\alpha_i$이고 1차식이라 $\nabla\alpha_i$는 요소 안에서 상수다. 따라서

$$W_e = \frac{\varepsilon}{2}\int_e \lvert \nabla V_e \rvert^2\, dS = \frac{\varepsilon}{2}\sum_{i,j} V_{ei}\, C^{(e)}_{ij}\, V_{ej}, \qquad C^{(e)}_{ij} = \int_e \nabla\alpha_i \cdot \nabla\alpha_j\, dS$$

$C^{(e)}$가 $3 \times 3$ 요소 계수 행렬(강성 행렬)이고, 적분이 상수 × 넓이라 닫힌 꼴이 있다. Sadiku의 식으로

$$C^{(e)}_{12} = \frac{1}{4A}\bigl[(y_2 - y_3)(y_3 - y_1) + (x_3 - x_2)(x_1 - x_3)\bigr], \qquad C^{(e)}_{11} = \frac{1}{4A}\bigl[(y_2 - y_3)^2 + (x_3 - x_2)^2\bigr]$$

이며 나머지는 첨자 순환이다. 각 행의 합이 0이다($\sum_j C_{ij} = 0$ — 상수 퍼텐셜은 에너지가 0). 1.3의 이차 형식(1.3.7이차 형식과 주축 변환)이 여기서 에너지의 얼굴로 나타난다.

⑤ 조립한다. 전체 에너지는 요소 에너지의 합 $W = \sum_e W_e = \frac{\varepsilon}{2}\mathbf{V}^{\mathsf{T}}C\mathbf{V}$이고, 전역 행렬 $C$는 각 요소의 $3 \times 3$을 그 요소의 절점 번호 자리에 더해 넣어 만든다. 한 절점이 여러 요소에 속하면 그만큼 겹쳐 더해진다. $C_{ij}$는 절점 $i$와 $j$가 같은 요소에 속할 때만 0이 아니므로 희소·대칭·(반)양의 정부호 — 11절의 이산 라플라시안과 같은 부류다.

⑥ 최소화. $\partial W/\partial V_k = 0$ (모든 자유 절점 $k$)에서 $C\mathbf{V} = \mathbf{0}$. 경계에서 값이 지정된 절점(고정 절점)과 나머지(자유 절점)로 나누면

$$\begin{bmatrix} C_{ff} & C_{fp} \\ C_{pf} & C_{pp} \end{bmatrix} \begin{bmatrix} \mathbf{V}_f \\ \mathbf{V}_p \end{bmatrix} \;\Longrightarrow\; C_{ff}\,\mathbf{V}_f = -C_{fp}\,\mathbf{V}_p$$

우변이 알려진 경계값이고, 좌변이 5·6절로 푸는 선형계다. 푸아송 방정식이면 우변에 전하 항 $\int \rho\,\alpha_i\, dS$가 더해진다.

직관

유한요소법이 차분법과 결정적으로 다른 점 둘.

노이만 조건이 공짜다. 에너지 $W$를 최소화하는 과정에서 경계 항이 $\oint V\,\partial V/\partial n\, dl$ 꼴로 나오는데, 값을 지정하지 않은 경계에서는 이것이 0이 되는 것이 최소화의 결과로 따라온다. 즉 아무것도 하지 않은 경계는 저절로 $\partial V/\partial n = 0$이다(자연 경계 조건). 11절에서 가상점을 두고 소거하던 수고가 없다.

해가 어디에 사는지 정해져 있다. 근사해는 언제나 “조각별 1차식”이라는 정해진 공간 안에 있고, 그 공간에서 에너지 기준으로 가장 좋은 원소다. 격자를 잘게 하면 공간이 넓어져 참해에 가까워지고, 오차는 에너지 노름에서 $O(h)$, 값에서 $O(h^2)$다. 2차 형상 함수 (변마다 절점을 하나씩 더 둔 6절점 삼각형)로 올리면 한 차수씩 오른다.

에너지가 없는 문제(비대칭 연산자, 이류항)에도 같은 뼈대가 통한다. 방정식 $\mathcal{L}u = f$의 잔차 $\mathcal{L}\tilde u - f$가 형상 함수 각각과 직교하도록 요구하는 것 — $\int (\mathcal{L}\tilde u - f)\,\alpha_i\, dS = 0$ — 이 갈러킨법이고, 에너지가 있으면 레일리–리츠와 같은 방정식을 준다. 상용 유한요소 도구(COMSOL, ANSYS)의 내부가 이것이고, 반도체에서는 채널 응력의 탄성 해석, 패키지의 열–기계 해석, 웨이퍼 휨이 유한요소법으로 풀린다. 소자 시뮬레이터가 쓰는 유한 체적법(11절)은 형상 함수 대신 상자를 쓴 사촌이고, 둘 다 비정형 격자를 다룬다.

모멘트법 — 적분 방정식을 행렬로

정전기 문제를 거꾸로 세울 수 있다. 퍼텐셜을 구하는 대신, 도체 표면의 전하 분포를 미지수로 잡는다. 전하가 만드는 퍼텐셜은 1.1·1.8의 적분(1.8.42·3차원 문제)으로 쓰이고, 도체 표면에서 그 값이 인가 전압과 같아야 한다.

$$V(\mathbf{r}) = \frac{1}{4\pi\varepsilon_0}\int_{L} \frac{\rho_L(\mathbf{r}')\, dl'}{\lvert \mathbf{r} - \mathbf{r}' \rvert} = V_0 \qquad (\mathbf{r} \text{는 도체 위})$$

미지 함수 $\rho_L$이 적분 안에 있다 — 적분 방정식이다. Sadiku의 예는 길이 $L$, 반지름 $a$의 가는 도선에 $V_0 = 1\ \mathrm{V}$를 건 경우다.

① 도선을 $N$개의 조각으로 나누고, 각 조각에서 $\rho_L$이 상수 $\rho_n$이라 놓는다 (펄스 기저). 미지수 $N$개.

② 각 조각의 중점 $\mathbf{r}_m$에서 방정식이 성립하도록 요구한다(점 정합, point matching). 방정식 $N$개.

$$\sum_{n=1}^{N} A_{mn}\,\rho_n = V_0, \qquad A_{mn} = \frac{1}{4\pi\varepsilon_0}\int_{\text{조각 } n} \frac{dl'}{\lvert \mathbf{r}_m - \mathbf{r}' \rvert} \approx \frac{\Delta}{4\pi\varepsilon_0\,\lvert x_m - x_n \rvert} \quad (m \neq n)$$

$\Delta = L/N$이 조각 길이다. $m = n$이면 분모가 0이라 근사가 깨지는데, 조각 안에서 도선의 반지름 $a$를 살려 적분하면 자기 항 $A_{mm} = \frac{2}{4\pi\varepsilon_0}\ln\frac{\Delta}{a}$이 나온다. 이 항이 대각을 지배하므로 행렬은 대각 우세다.

③ $N \times N$ 선형계를 푼다. $A$는 꽉 찬(dense) 행렬이다 — 모든 조각이 모든 조각에 퍼텐셜을 만들기 때문이다. 그래서 $N$을 크게 못 잡지만($10^{3}$ 급), 미지수가 표면에만 있어 3차원 문제가 2차원 미지수로 줄어든다. 얻은 $\rho_n$을 더하면 총전하 $Q = \sum \rho_n\Delta$이고, $C = Q/V_0$이 도선의 커패시턴스다.

직관

언제 어느 쪽인가. 차분법·유한요소법은 영역을 채우고, 모멘트법은 경계만 채운다. 그래서 —

• 영역 안이 불균일하면(도핑·유전율·재료가 위치마다 다른 소자) 영역법. 모멘트법은 균일한 매질의 그린 함수가 있어야 성립한다.
• 영역이 무한히 열려 있으면(안테나, 패키지 밖의 전자기장, 배선 사이의 정전 결합) 모멘트법. 영역법은 무한 영역을 어딘가에서 인공적으로 잘라야 하고(11절의 인공 노이만 경계), 그 절단이 오차다.

배선 기생 커패시턴스 추출이 두 번째다 — 도체 표면에 전하를 놓고 모든 쌍의 상호 커패시턴스를 구한다. $N$이 커지면 꽉 찬 행렬을 반복법(6절)으로 풀되 행렬–벡터 곱을 다중극 전개로 $O(N)$에 근사하는 것(빠른 다중극법)이 그 도구의 심장이다. 밴드 계산의 평면파 기저와 원자 궤도 기저의 선택도 같은 대비다 — 전역 기저인가 국소 기저인가.

일반화. 펄스 기저와 점 정합은 가장 거친 선택이다. 미지 함수를 기저 $\{f_n\}$으로 전개하고 잔차를 가중 함수 $\{w_m\}$과 내적해 0으로 놓는 것이 일반적인 모멘트법이고, $w_m$이 델타 함수이면 점 정합, $w_m = f_n$이면 갈러킨이다. 즉 유한요소법과 모멘트법은 같은 가중 잔차법의 두 얼굴이고, 다른 것은 연산자가 미분(국소·희소)이냐 적분(전역·꽉 참)이냐뿐이다. “모멘트”라는 이름은 잔차와 가중 함수의 내적을 확률의 모멘트에 빗댄 것이다.

최급강하법과 최적화

방정식이 아니라 최솟값

지금까지는 “이 식을 만족하는 $\mathbf{x}$“를 찾았다. 마지막 절은 “이 값을 가장 작게 하는 $\mathbf{x}$“다. 8절의 최소제곱이 이미 그것이었고, 13절의 에너지 최소화도 그것이었다. 반도체에서 최적화는 어디에나 있다 — 소자 모델 파라미터 추출(측정과의 차이를 최소로), 공정 조건 최적화(수율을 최대로), 광 근접 보정(웨이퍼 패턴과 설계의 차이를 최소로), 그리고 신경망 학습(손실 함수를 최소로).

$\min f(\mathbf{x})$에서 $f$를 목적 함수라 한다. 미분 가능하면 최솟값에서 $\nabla f = \mathbf{0}$이다(필요조건). 이것은 $n$개의 미지수에 대한 $n$개의 방정식이니 2절의 뉴턴법으로 풀 수 있다 — 그러면 도함수 $J$가 $f$의 헤세 행렬 $H_{ij} = \partial^2 f/\partial x_i\partial x_j$가 된다. 그런데 헤세 행렬은 $n^2$개의 2계 도함수라 비싸고, 멀리서는 뉴턴법이 튄다. 그래서 1계 도함수만 쓰는 방법이 먼저 온다.

최급강하법 — 가장 가파른 내리막으로

현재 점에서 $f$가 가장 빨리 줄어드는 방향은 $-\nabla f$다(1.1.2기울기·발산·회전 — 기울기는 가장 빨리 커지는 방향). 그 방향으로 얼마나 갈지($t$)는 그 직선 위에서 $f$가 최소가 되는 곳으로 정한다.

$$\mathbf{x}_{k+1} = \mathbf{x}_k - t_k\,\nabla f(\mathbf{x}_k), \qquad t_k = \arg\min_{t > 0}\ g(t), \quad g(t) = f\bigl(\mathbf{x}_k - t\,\nabla f(\mathbf{x}_k)\bigr)$$

$g(t)$는 1변수 함수라 2절의 방법이나 해석적으로 최소화한다(선 탐색). Kreyszig의 예로 $f = x_1^2 + 3x_2^2$, 시작점 $(6, 3)$이면 $\nabla f = (2x_1, 6x_2)$이고 $g(t)$가 $t$의 2차식이라 $t_k$가 닫힌 꼴로 나온다.

$$(6,\ 3) \to (3.484,\ -0.774) \to (1.327,\ 0.664) \to (0.771,\ -0.171) \to \cdots$$

$x_2$의 부호가 매번 뒤집히며 지그재그로 원점에 다가간다. 수렴은 하되 느리다.

그림 1.12-10. f = x₁² + 3x₂²의 등고선(타원) 위에서 최급강하법이 걷는 길. 각 걸음은 등고선에 수직인 방향으로 출발해, 그 직선 위에서 f가 최소인 곳 — 직선이 어떤 등고선에 접하는 곳 — 에서 멈춘다. 다음 걸음은 그 접점의 등고선에 수직이므로 직전 걸음과 정확히 직각이다. 그래서 길이 지그재그가 되고, 타원이 길쭉할수록(조건수가 클수록) 골짜기 벽을 오가느라 바닥을 따라 내려가지 못한다. 켤레 기울기법은 두 번째 걸음부터 방향을 보정해 2차원 이차 함수를 정확히 두 걸음에 끝낸다.직접 작도

직관

왜 지그재그인가. 선 탐색이 끝난 점에서는 직선 방향으로 $f$가 더 안 줄어드니 $\nabla f \cdot (\text{직선 방향}) = 0$ — 새 기울기가 직전 걸음과 수직이다. 매 걸음이 직전과 직각이면 좁은 골짜기에서는 벽에서 벽으로 튕길 수밖에 없다.

얼마나 느린가. $f = \tfrac12\mathbf{x}^{\mathsf{T}}A\mathbf{x}$(대칭 양의 정부호 $A$, 등고선이 타원)에서 최급강하의 오차 감소율은

$$\frac{\lVert \mathbf{e}_{k+1} \rVert_A}{\lVert \mathbf{e}_k \rVert_A} \le \frac{\kappa - 1}{\kappa + 1}, \qquad \kappa = \frac{\lambda_{\max}}{\lambda_{\min}}$$

7절의 조건수가 다시 나온다. 타원의 긴 축과 짧은 축의 비가 $\sqrt{\kappa}$이고, $\kappa = 1$ (원)이면 한 걸음에 끝나며, $\kappa = 100$이면 걸음마다 오차가 $0.98$배 — 자릿수 하나에 백 번이다. 6절의 반복법과 정확히 같은 구조이고, 실제로 $A\mathbf{x} = \mathbf{b}$를 푸는 일은 $\tfrac12\mathbf{x}^{\mathsf{T}}A\mathbf{x} - \mathbf{b}^{\mathsf{T}}\mathbf{x}$를 최소화하는 일과 같다. 켤레 기울기법은 각 걸음의 방향을 직전 방향과 $A$-직교($\mathbf{d}_{k+1}^{\mathsf{T}}A\mathbf{d}_k = 0$)가 되게 고쳐 지그재그를 없앤 것이고, 감소율이 $(\sqrt\kappa - 1)/(\sqrt\kappa + 1)$로 좋아진다 — $\kappa$가 $\sqrt\kappa$로. 6절에서 “반복 횟수가 $\sqrt\kappa$에 비례한다”고 한 것이 이 식이다.

뉴턴과 준뉴턴

목적 함수를 2차까지 전개하면 $f(\mathbf{x} + \delta) \approx f + \nabla f^{\mathsf{T}}\delta + \tfrac12\delta^{\mathsf{T}}H\delta$이고, 이것을 $\delta$에 대해 최소화하면

$$H\,\delta = -\nabla f$$

뉴턴법이다. 이차 함수라면 한 걸음에 정확히 최솟값에 간다 — 최급강하가 골짜기 벽을 오갈 때 뉴턴은 골짜기의 모양($H$)을 알고 바닥으로 곧장 간다. 해 근처에서 2차 수렴. 그러나 $H$가 양의 정부호가 아닌 곳(안장점, 오목한 곳)에서는 그 방향이 오르막일 수 있고, $H$를 만들고 푸는 비용이 $n^2$과 $n^3$이다.

준뉴턴법(BFGS)은 $H$를 계산하지 않고, 걸음마다 관측된 기울기의 변화 $\nabla f_{k+1} - \nabla f_k$로 $H^{-1}$의 근사를 갱신해 간다 — 2절의 할선법을 벡터로 올린 것이다. 1계 도함수만으로 초선형 수렴을 얻고, 변수가 수천 개인 곡선 맞춤·모델 보정의 표준이다. 8절의 레벤버그–마쿼트는 최소제곱에 특화된 같은 부류로, $H \approx J^{\mathsf{T}}J$로 놓은 것이 다르다.

변수가 수백만 개(신경망, 픽셀 단위의 마스크 최적화)면 $H^{-1}$ 근사조차 못 들고 있으므로 최급강하로 돌아간다 — 다만 선 탐색 대신 고정 걸음(학습률)과 운동량 항으로 지그재그를 누르고, 기울기를 데이터 일부로만 추정한다(확률적 경사 하강). 걸음 폭이 6절의 $\omega$·10절의 $h$와 같은 역할이라 너무 크면 발산하고 너무 작으면 느리다.

구속 조건 — 라그랑주 승수와 선형 계획

$g(\mathbf{x}) = 0$을 만족하면서 $f$를 최소화하는 문제는 1.6의 라그랑주 승수(1.6.7구속 조건이 있는 변분과 라그랑주 승수)로 $\nabla f = \lambda\nabla g$를 풀거나, 구속을 벌칙 항 $\mu g^2$으로 목적 함수에 더해 무구속 문제로 바꾼다. 부등식 구속($g \le 0$)이면 승수의 부호 조건이 붙는다(KKT 조건).

$f$와 구속이 전부 선형이면 선형 계획법(Kreyszig 22.2)이다. 등고선이 평행한 직선이고 가능 영역이 다각형(다면체)이므로, 최솟값은 언제나 꼭짓점에 있다 — 안쪽에서는 기울기 방향으로 더 갈 수 있고, 변 위에서도 한쪽 꼭짓점으로 가면 더 좋아지기 때문이다. 심플렉스법은 꼭짓점에서 이웃 꼭짓점으로 목적 함수가 줄어드는 쪽으로만 옮겨 가는 방법이고, 웨이퍼 배정· 장비 스케줄링 같은 생산 계획이 이것으로 풀린다.

국소와 전역

이 절의 방법은 전부 내리막으로만 간다. 그래서 도착하는 곳은 시작점에서 가장 가까운 골짜기이지 가장 깊은 골짜기가 아니다. 목적 함수가 볼록(등고선이 한 겹의 타원처럼 생김)이면 골짜기가 하나라 문제없다 — 8절의 선형 최소제곱, 13절의 에너지가 그렇다. 비선형 모델 맞춤, 공정 최적화, 역 리소그래피는 볼록하지 않고 골짜기가 여럿이다.

처방은 세 종류다. 여러 시작점에서 출발해 가장 좋은 것을 취하거나(다중 시작), 오르막도 확률적으로 허용해 얕은 골짜기를 빠져나오거나(모의 담금질 — 온도를 서서히 낮추는 열처리의 비유가 이름이다), 목적 함수 자체를 몇 번의 평가로 대리 모델로 만들어 그 위에서 탐색한다 (베이즈 최적화 — 시뮬레이션 한 번이 몇 시간인 TCAD 보정에서 쓴다). 어느 쪽이든 찾은 최솟값이 전역이라는 보장은 없다는 것이 이 절의 마지막 경고다.

의문

파라미터 추출에서 “최적화가 수렴했다”는 말은 8절의 경고에 이 절의 경고를 더한 것이다. 잔차가 작아도 파라미터가 보상하고 있을 수 있고(SVD의 작은 특이값), 그 잔차가 국소 최솟값의 잔차일 수도 있다(다른 시작점에서 더 작은 잔차가 나온다). 그래서 실무의 추출은 파라미터를 한꺼번에 풀지 않는다 — 문턱전압은 선형 영역 낮은 $V_D$에서, 이동도는 그다음, 포화 파라미터는 그다음 식으로 한 번에 한두 개씩, 앞 단계의 값을 고정한 채 맞춘다. 골짜기가 하나인 작은 문제의 연쇄로 쪼개는 것이고, 2절의 검멜 반복·6절의 블록 반복과 같은 발상이다.

이 장이 반도체에서 다시 나오는 곳

• 소자 시뮬레이션(TCAD). 유한 체적 이산화(11절) → 셰르페터–검멜(12절) → 결합 뉴턴 또는 검멜 반복(2·6절) → 스케일링된 희소 야코비 행렬(7절)을 직접법·크릴로프법으로(5·6절) → 과도 해석은 후진 오일러·BDF(10절). 격자는 디바이 길이가 정한다(11절).
• 회로 시뮬레이션(SPICE). 수정 노드 해석의 비선형 대수–미분 방정식을 사다리꼴·기어법으로 이산화(10절)하고 걸음마다 뉴턴–LU(2·5절). RELTOL·ABSTOL은 1절의 상대·절대 오차, “timestep too small”은 10절의 오차 추정과 2절의 수렴 실패.
• 공정 시뮬레이션. 확산·산화는 비선형 포물형 PDE — 크랭크–니콜슨·후진 오일러(12절), 격자 이동. 이온 주입은 몬테카를로(1.11).
• 밴드 구조·양자 소자. 큰 에르미트 행렬의 가장 낮은 고유값 몇 개 — 란초스·다비드슨(9절). 슈뢰딩거–푸아송 자기일관 반복은 2절의 고정점 반복이고, 혼합(mixing) 파라미터가 6절의 $\omega$다.
• 계측과 파라미터 추출. 곡선 맞춤은 레벤버그–마쿼트(8절), C–V 미분은 4절의 경고, TLM·아레니우스는 직선 맞춤(8절), 신뢰구간과 파라미터 보상은 SVD(8절)와 국소 최솟값(14절).
• 배선·패키지. 기생 커패시턴스·인덕턴스 추출은 모멘트법 + 빠른 다중극(13절), 응력·열은 유한요소법(13절), 신호 무결성은 FDTD — 12절의 CFL 조건이 시간 걸음을 정한다.
• 리소그래피. 광 근접 보정과 역 리소그래피는 픽셀 단위 최급강하(14절), 공정 창 최적화는 구속 최적화.

전부 “연속을 이산으로 바꾸고 무엇을 잃는지 잰다”의 변주이고, 잃는 것은 언제나 절단·반올림· 반복·조건의 네 층 중 하나다.

이 장이 쓰이는 곳

• 1.11 확률과 통계
• 2.2 경계값 문제