편미분방정식의 수치 해법

격자와 유한 차분 근사

구간을 간격 $h = 1/N$로 나눈 격자(grid, mesh) $x_i = ih$ 위에서 $u_i \approx u(x_i)$를 미지수로 둡니다. 9장의 중앙 차분을 2계 도함수에 쓰면, 테일러 전개 $u(x \pm h) = u \pm hu' + \frac{h^2}{2}u'' \pm \frac{h^3}{6}u''' + \frac{h^4}{24}u'''' + \cdots$를 더해 홀수 항이 지워지고

이를 각 내부 점에 쓰면 $-u_{i-1} + 2u_i - u_{i+1} = h^2 f_i$, $i = 1, \dots, N - 1$이고, 경계값 $u_0 = u_N = 0$을 넣으면 삼중대각 연립방정식 $T\mathbf u = h^2\mathbf f$입니다. 행렬 $T = \operatorname{tridiag}(-1, 2, -1)$은 대칭 양정치이며, 12장의 토머스 알고리즘으로 $O(N)$ 연산에 풀립니다.

푸아송 방정식: 5점 스텐실과 희소 선형계

$x$와 $y$ 방향에 앞 절의 중앙 차분을 각각 쓰면, 간격 $h$의 정사각 격자 점 $(x_i, y_j)$에서

한 변에 내부 점이 $n = N - 1$개면 미지수는 $n^2$개, 행렬은 $n^2\times n^2$입니다. $n = 100$이면 미지수 9801개($N = 100$), 밀집 행렬로 저장하면 약 $10^8$개의 칸이지만 0이 아닌 값은 5만 개 미만입니다. 그래서 행렬을 만들지 않고 스텐실만으로 계산하는 희소 방법이 중요합니다. 밀집 가우스 소거는 $(n^2)^3/3 \approx 3\times10^{11}$번, 띠 촐레스키는 띠 너비 $n$을 이용해 약 $n^2\cdot n^2/2 = 5\times10^7$번입니다. 3차원이면 차이가 훨씬 커집니다.

큰 희소계를 푸는 방법은 13장의 주제였습니다. 이 문제에서 야코비 반복의 스펙트럼 반지름은 정확히 $\rho_J = \cos\pi h$, 가우스–자이델은 $\rho_J^2$, 최적 SOR($\omega = 2/(1 + \sin\pi h)$)은 $\omega - 1$입니다. 오차를 $10^{-6}$배 줄이는 데 필요한 반복 수는 다음과 같습니다.

야코비는 $N^2$에, SOR은 $N$에 비례해 늘어납니다. 격자 간격이 작아지면 행렬의 조건수가 $\kappa \approx 4/(\pi h)^2$로 커지기 때문입니다. 켤레 기울기법도 $\sqrt\kappa \propto N$번이 들고, 격자 크기와 무관한 반복 수로 푸는 것이 13장의 다중 격자법입니다. 직사각형 영역의 상수 계수 문제라면 8장의 FFT로 사인 변환을 해서 $O(n^2\log n)$에 직접 풀 수도 있습니다.

열 방정식: 명시적 FTCS, 음적 BTCS, 크랭크–니콜슨

공간을 먼저 차분화하면 열 방정식 $u_t = \alpha u_{xx}$는 내부 점 $n$개의 연립 ODE $\mathbf u' = -\frac{\alpha}{\Delta x^2}T\mathbf u$가 됩니다(선의 방법(method of lines)). $T$의 고유값은 $4\sin^2\frac{k\pi\Delta x}{2}$($k = 1, \dots, n$)로, 가장 큰 것이 약 4라 이 ODE의 고유값은 $0$에 가까운 것부터 $-4\alpha/\Delta x^2$까지 퍼져 있습니다. 강성비가 $\Delta x^{-2}$에 비례하는 전형적인 강성 계입니다. 이 ODE에 11장의 세 방법을 쓰면 세 가지 PDE 방법이 나옵니다.

크랭크–니콜슨 방법은 1947년 크랭크와 니콜슨(John Crank, Phyllis Nicolson)이 발표했습니다. 음적 방법의 대가는 매 걸음의 선형계인데, 1차원에서는 토머스 알고리즘(12장)으로 $O(n)$, 곧 FTCS 한 걸음과 같은 차수의 비용입니다. 다차원에서는 희소 선형계를 반복법이나 교대 방향 음적법(ADI: 방향마다 삼중대각계를 번갈아 푸는 방법)으로 풉니다.

function crank_nicolson_step(u, r)          // u: 내부 점 n개, 양끝 0
  for i = 1 .. n                             // 오른쪽: (I - r/2 T) u
    d[i] = u[i] + (r/2) * (u[i-1] - 2*u[i] + u[i+1])
  sub = -r/2;  diag = 1 + r;  sup = -r/2      // 왼쪽: (I + r/2 T)
  return thomas(sub, diag, sup, d)          // O(n): 12장

CN에도 약점이 있습니다. 가운데 0.4~0.6만 1이고 나머지가 0인 사각 펄스에서 $\Delta t = 0.05$($r = 20$)로 한 걸음 가면, $t = 0.05$의 참값 중앙 온도 0.244 대신 CN은 $-0.235$라는 음의 온도를 내고 이웃 점과 부호가 번갈아 바뀝니다. BTCS는 0.281로 조금 퍼졌지만 매끄럽습니다. 다음 절의 증폭 인자가 이유를 보여 줍니다.

폰 노이만 안정성 분석

ODEBook 16장은 FTCS의 증폭 인자 $G(\theta) = 1 - 4r\sin^2\frac\theta2$를 구해 $r \le \frac12$을 얻었습니다. 같은 계산을 모든 선형 상수 계수 차분법에 할 수 있습니다. 이것이 폰 노이만 안정성 분석(von Neumann stability analysis)입니다.

$\delta^2 e^{\mathrm ij\theta} = (e^{\mathrm i\theta} - 2 + e^{-\mathrm i\theta})e^{\mathrm ij\theta} = -4s^2e^{\mathrm ij\theta}$, $s = \sin\frac\theta2$를 쓰면 세 방법의 증폭 인자가 곧바로 나옵니다.

정리의 둘째 부분이 앞 절 CN의 지그재그를 설명합니다. $r = 20$이면 $G_{\text{CN}}(\pi) = -39/41 = -0.951$. 톱니 모드는 사라지지 않고 걸음마다 부호만 바꾸며 거의 그대로 남습니다. 참 해에서는 이 모드가 $e^{-20\pi^2}$배, 사실상 즉시 사라져야 합니다. 사각 펄스처럼 높은 모드를 많이 가진 자료에서 이것이 눈에 보입니다. BTCS처럼 $|G| \to 0$인 성질을 11장에서 L-안정이라 불렀습니다. 실무에서는 처음 몇 걸음을 BTCS로 걸어 높은 모드를 죽인 뒤 CN으로 바꾸는 요령(래너처 시작)을 씁니다.

폰 노이만 분석은 이류 방정식에도 그대로 쓰입니다(6절). 아래 시뮬레이터에서 두 방정식의 증폭 인자를 비교하세요.

랙스 동등성 정리: 일관성 + 안정성 = 수렴

수치 해가 참 해에 다가가는 것(수렴)은 우리가 원하는 성질이지만 직접 확인하기 어렵습니다. 참 해를 모르기 때문입니다. 반면 일관성은 테일러 전개로, 안정성은 폰 노이만 분석으로 확인할 수 있습니다. 다음 정리가 둘을 수렴에 잇습니다.

이 정리에서 실용적인 결론 두 가지가 나옵니다. 첫째, 안정한 방법의 수렴 차수는 절단 오차의 차수와 같습니다. 예제 16.3의 BTCS 1차, CN 2차가 그 예입니다. 그래서 새 코드를 검증할 때는 (1) 절단 오차의 차수를 손으로 구하고, (2) 폰 노이만으로 안정성을 확인한 뒤, (3) 격자 수렴 실험으로 그 차수가 실제로 나오는지 봅니다. 둘째, 일관적이지만 불안정한 방법은 어떤 격자에서도 쓸모가 없습니다. 격자를 줄이면 국소 오차는 줄지만 증폭 횟수가 늘어 결국 반올림 오차까지 키웁니다.

이류 방정식: 풍상 차분, 수치 확산, 랙스–벤드로프, CFL

$C = a\Delta t/\Delta x$를 쿠랑 수라 합니다. 시간 전진·공간 중앙 차분 $u_j^{k+1} = u_j^k - \frac C2(u_{j+1}^k - u_{j-1}^k)$의 증폭 인자는 $G = 1 - \mathrm iC\sin\theta$이고 $|G|^2 = 1 + C^2\sin^2\theta > 1$입니다. 어떤 $\Delta t$를 써도 불안정합니다. 열 방정식에는 확산이 높은 모드를 눌러 주지만, 이류에는 그런 힘이 없어 중앙 차분의 작은 증폭이 그대로 쌓입니다.

처방은 정보가 오는 쪽을 보는 것입니다. $a > 0$이면 해는 왼쪽에서 옵니다.

풍상 차분은 안정하지만 대가가 있습니다. 풍상 식에 참 해를 넣고 테일러 전개하면, 이 방법이 실제로 더 정확하게 푸는 방정식(수정 방정식(modified equation))은

랙스–벤드로프는 이 확산항을 지웠지만 다음 항인 3계 도함수 항(분산)이 남아, 파장마다 속도가 달라집니다. 짧은 파장이 뒤처지며 급한 경사 뒤에 물결 같은 진동이 생깁니다. 실제로 계산해 보면 다음과 같습니다.

이 딜레마는 우연이 아닙니다. 고두노프(Sergei Godunov)는 새 극값을 만들지 않는(단조) 선형 방법은 많아야 1차 정확이라는 것을 보였습니다. 현대의 고해상도 방법(TVD, 플럭스 제한자, WENO)은 매끄러운 곳에서는 2차 이상, 급한 경사 근처에서는 1차로 바뀌는 비선형 방법으로 이 벽을 넘습니다.

CFL 조건의 의미는 "한 걸음에 정보가 한 칸보다 멀리 가면 안 된다"입니다. 빠른 흐름($a$ 큼)이나 가는 격자($\Delta x$ 작음)는 시간 걸음을 직접 줄입니다. 3차원 CFD 코드가 비싼 이유의 상당 부분이 이것이고, 음적 방법이 여기서도 쓰입니다. 다만 이류에서는 큰 $C$의 음적 방법이 정확도를 크게 잃으므로, 열 방정식만큼 쉽게 걸음을 키우지는 못합니다.

파동 방정식: 립프로그와 수치 분산

ODEBook 16장은 파동 방정식 $u_{tt} = c^2u_{xx}$의 중앙 차분 $u_i^{k+1} = 2u_i^k - u_i^{k-1} + C^2\delta^2u_i^k$와 CFL 조건 $C = c\Delta t/\Delta x \le 1$을 보였습니다. 시간과 공간 모두 중앙 차분이라 앞뒤 두 시각을 "뛰어넘는" 꼴이어서 흔히 립프로그(leapfrog) 방법이라 부릅니다. 여기서는 그 다음 질문, 안정한 범위 안에서도 무엇이 틀리는가를 봅니다.

시작에는 두 시각이 필요합니다. $u^0 = f$와 $u_t(x, 0) = g$에서 테일러 전개로 $u_i^1 = u_i^0 + \Delta t\,g_i + \frac{C^2}{2}\delta^2u_i^0$을 만들면 2차 정확도가 유지됩니다. 모드 $u_j^k = e^{\mathrm i(\xi j\Delta x - \omega k\Delta t)}$를 넣으면 수치 분산 관계

가 나옵니다. 수치 파의 속도 $c_{\text{num}}/c = \frac{2}{C\theta}\arcsin\bigl(C\sin\frac\theta2\bigr)$, $\theta = \xi\Delta x$는 파장에 따라 달라집니다. 계산하면

짧은 파장일수록 늦고, $C$가 작을수록 더 늦습니다. 급한 펄스는 여러 파장의 합이므로, 시간이 지나면 짧은 성분이 뒤처져 펄스 뒤에 잔물결 꼬리가 생깁니다. 이것이 수치 분산(numerical dispersion)입니다. 흥미롭게도 1차원에서는 $C = 1$이 가장 정확하고(격자점에서 정확한 해), 안정 경계에서 가장 정확하다는 것은 다차원에서는 성립하지 않습니다. 지진파·음향 시뮬레이션은 "파장당 몇 칸"을 기준으로 격자를 고르고, 분산을 줄이려고 고차 스텐실을 씁니다.

유한 요소법 입문: 약형식, 모자 함수, 강성 행렬

$-u'' = f$, $u(0) = u(1) = 0$에 경계에서 0인 아무 시험 함수 $v$를 곱해 적분하고 부분 적분하면 경계항이 사라져

유한 요소법(finite element method, FEM)은 $u$와 $v$를 유한 차원 공간에서 고릅니다. 가장 단순한 공간은 마디 $x_i$에서 1, 다른 마디에서 0이고 그 사이는 직선인 모자 함수(hat function) $\varphi_i$들로 만든 조각별 선형 함수입니다. $u_h = \sum_j U_j\varphi_j$를 넣고 $v = \varphi_i$로 두면(갤러킨 방법)

2차원에서는 삼각형 요소 위의 조각별 선형 함수를 쓰고, 요소 강성 행렬은 3×3이 됩니다. 정사각 격자의 직각삼각형으로 라플라스 방정식을 풀면 5점 스텐실과 같은 행렬이 나옵니다. 유한 요소법은 차분법을 포함하면서 임의의 모양과 재료, 고차 요소, 오차 추정에 따른 적응 세분화로 확장됩니다. 구조 해석의 상용 코드, 전자기장 해석, 열 해석 대부분이 이 방법입니다.

유한 체적법 맛보기: 보존을 정확히

유체의 질량·운동량·에너지, 반도체의 전하는 보존됩니다. 칸의 양이 변하는 것은 칸의 경계를 드나드는 흐름(플럭스) 때문뿐입니다. 보존 법칙 $u_t + F(u)_x = 0$을 칸 $[x_{j-1/2}, x_{j+1/2}]$에서 적분하면 칸 평균 $\bar u_j$에 대해

이것이 유한 체적법(finite volume method)의 뼈대입니다. 무엇을 $F_{j+1/2}$로 쓰느냐가 방법을 정합니다. 이류 $F = au$에서 바람이 불어오는 쪽 칸의 값 $F_{j+1/2} = a\bar u_j$를 쓰면 정확히 앞 절의 풍상 차분이 되고, 예제 16.4의 표에서 모든 방법의 총량이 0.200000으로 보존된 것도 그 식들이 이런 플럭스 차이 꼴이기 때문입니다. 비선형 보존 법칙(충격파가 생기는 기체 역학)에서는 칸 경계마다 작은 리만 문제를 풀어 플럭스를 정하는 고두노프 계열 방법이 표준이고, 보존형으로 쓰지 않으면 충격파의 속도가 틀리게 나옵니다. 임의의 다각형 칸에서도 경계 면의 플럭스만 계산하면 되므로, 상용·공개 CFD 코드 대부분이 유한 체적법을 씁니다.

응용: 열 해석, 유체, 반도체 소자

• 전자 기기의 열 해석. 칩과 기판의 온도는 발열원이 있는 열 방정식 $\rho c\,u_t = \nabla\cdot(k\nabla u) + q$를 따릅니다. 재료마다 열전도도 $k$가 수백 배 다르고 층이 얇아 격자가 매우 불균일하므로 강성이 심합니다. 그래서 BTCS나 CN 같은 음적 시간 진행과, 매 걸음 큰 희소 선형계를 푸는 전처리 켤레 기울기법(13장)의 조합이 표준입니다. 정상 상태만 원하면 2절의 푸아송 문제 하나입니다.
• 전산 유체 역학. 나비에–스토크스 방정식은 이류(6절)·확산(3절)·압력의 푸아송 방정식(2절)을 모두 품고 있습니다. 많은 코드가 유한 체적법으로 보존을 지키고, 이류 항에는 풍상 계열의 고해상도 방법을, 시간 걸음은 CFL 수로 고르며, 매 걸음 압력의 푸아송 방정식을 다중 격자법으로 풉니다. 수치 확산이 실제 점성보다 커지지 않게 격자를 고르는 것이 결과를 믿기 위한 첫 점검입니다.
• 반도체 소자 시뮬레이션. 트랜지스터 안의 전위는 전하 밀도를 원천으로 하는 푸아송 방정식을, 전자·정공 농도는 표류–확산 연속 방정식을 따릅니다. 농도가 깊이에 따라 수십 자릿수 변하기 때문에 보통의 중앙 차분은 음의 농도를 만들기 쉽습니다. 1969년 샤르페터와 굼멜(Scharfetter, Gummel)이 제안한 이산화는 칸 사이의 전류를 지수 함수 꼴로 정확히 맞춰 이 문제를 피하며, 오늘날 소자 시뮬레이터의 기본 구성 요소로 알려져 있습니다. 결과 연립방정식은 비선형이라 4장의 뉴턴 방법으로, 그 안의 선형계는 희소 직접법이나 반복법으로 풉니다.

이 도구가 쓰이는 곳

맺으며: 오차를 예측하고, 통제하고, 비용과 바꾸기

이 책은 한 문장에서 출발했습니다. 컴퓨터는 유한한 수와 유한한 단계로 근사한다. 좋은 수치해석은 오차를 예측하고, 통제하고, 비용과 바꾸는 기술이다. 1장에서 그 "유한한 수"는 0.1 + 0.2가 0.3이 아니게 만드는 기계 엡실론 $2^{-52}$이었고, 이 장의 마지막에서 "유한한 단계"는 판을 자른 격자와 시간 걸음이었습니다. 무대가 숫자 하나에서 수백만 개의 미지수로 넓어졌을 뿐, 묻는 질문은 같았습니다.

우리는 그 질문을 몇 개의 공통 언어로 다루었습니다. 오차 차수 $O(h^p)$는 보간(5–6장), 미분과 적분(9–10장), ODE와 PDE(11장, 이 장)에서 로그–로그 그래프의 기울기로 되풀이해 나타났습니다. 수렴 차수는 이분법의 1차, 뉴턴의 2차, 할선법의 1.618(3–4장)에서 시작해 반복 해법과 최적화의 수렴률 $(\kappa - 1)/(\kappa + 1)$과 $\sqrt\kappa$ 가속(13장, 15장)으로 이어졌습니다. 조건수는 2장에서 문제의 민감도로 등장해 연립방정식의 오차 한계(12장), 최소제곱의 QR(7장), 경사 하강의 지그재그(15장), 격자를 줄일수록 나빠지는 푸아송 행렬(이 장)까지 따라왔습니다. 그리고 안정성은 2장의 점화식과 근의 공식에서 11장의 안정성 영역, 이 장의 증폭 인자와 랙스 동등성 정리까지, "작은 오차가 자라지 않는가"라는 한 질문의 여러 얼굴이었습니다.

여기서 길은 여러 갈래로 이어집니다. 다차원 유한 요소법과 적응 격자, 비선형 보존 법칙과 충격파, 확률 미분방정식과 몬테카를로, 큰 규모의 병렬 계산과 혼합 정밀도, 그리고 데이터로 학습한 모델과 물리 방정식을 결합하는 계산이 모두 이 책의 질문을 다시 묻습니다. 새로운 알고리즘을 만나면 먼저 물어보세요. 오차는 얼마나 빨리 줄어드는가, 작은 오차가 자라지는 않는가, 그 정확도에 얼마를 내는가. 세 질문에 숫자로 답할 수 있으면, 그 계산을 믿어도 됩니다.