Chapter 10

컴퓨터는 어떻게 틀리는가

컴퓨터는 계산을 틀리지 않는다고들 합니다. 실제로는 거의 모든 계산에서 조금씩 틀리고, 그 '조금'이 어떤 조건에서 커지는지 아는 사람이 시뮬레이터를 믿을 만하게 만듭니다. 이 장은 숫자가 비트로 저장되는 방식에서 출발해 뺄셈 한 번이 유효숫자를 날리는 과정, 미분을 차분으로 바꿀 때 생기는 두 종류의 오차, 시간 간격을 조금만 크게 잡아도 시뮬레이션이 폭발하는 이유를 차례로 봅니다. 끝에서는 같은 원리가 이미지 센서 픽셀을 설계하는 전자기 시뮬레이터 FDTD와 RCWA, 그리고 연립방정식의 조건수로 이어지는 것을 확인합니다.

IEEE 754 배정밀도가 수를 부호·지수·가수로 저장하는 방식을 알고, 0.1 + 0.2 ≠ 0.3의 이유를 비트 수준에서 설명한다.
부동소수점 간격이 크기에 비례한다는 것(상대 정밀도 일정)과 기계 엡실론의 의미, AI에서 쓰는 bf16·fp8의 절충을 안다.
상쇄 오차를 알아보고 1 − cos x, 근의 공식을 안정적인 형태로 바꿔 쓴다.
수치 미분에서 절단 오차와 반올림 오차의 줄다리기(U자 곡선)와 최적 간격 √ε, ∛ε을 설명한다.
유한 차분으로 파동 방정식을 풀고 CFL 조건을 직접 깨 본다. FDTD의 Yee 격자, RCWA의 푸리에 전개 개념을 안다.
조건수가 입력 오차를 얼마나 키우는지 기하적으로 이해한다.

이 수식이 없었다면1991년 다란의 패트리어트 미사일

1991년 2월 25일 걸프전, 사우디아라비아 다란에 배치된 패트리어트 포대가 날아오는 스커드 미사일을 요격하지 못했고, 미사일은 미군 막사에 떨어져 28명이 숨졌습니다. 미국 회계감사원(GAO)의 1992년 보고서가 짚은 원인은 계산 오차였습니다. 시스템 시계는 시간을 0.1초 단위로 세는데, 0.1은 2진수로 끝나지 않는 소수입니다. 24비트 레지스터에 잘려 들어간 0.1은 진짜보다 약 0.000000095초 작았고, 포대가 약 100시간 연속 가동되는 동안 이 오차가 360만 번 쌓여 약 0.34초가 되었습니다. 초속 1.6 km가 넘는 스커드에게 0.34초는 수백 m입니다. 레이더가 다음에 표적을 찾아야 할 '추적 창'이 엉뚱한 곳으로 옮겨 갔고, 시스템은 표적이 없다고 판단했습니다.

수치 오차는 교과서 연습 문제가 아닙니다. 이 장의 도구는 '얼마나 틀릴 수 있는가'를 계산하는 수학입니다.

0.1초 시계→2진 무한소수→잘림 오차 누적→오차 해석·부동소수점 모델→믿을 수 있는 시뮬레이션

0.1 + 0.2 ≠ 0.3: 수는 비트로 저장된다

계산기가 틀린 답을 낸다

브라우저 콘솔(F12)에 0.1 + 0.2 === 0.3을 입력하면 false가 나옵니다. 0.1 + 0.2의 값은 0.30000000000000004입니다. 파이썬, C, 엑셀도 같습니다. 버그가 아니라 설계입니다. 왜 이렇게 설계했을까요?

10진수에서 1/3 = 0.333…이 끝나지 않듯, 2진수에서는 분모에 2 말고 다른 소인수(예: 5)가 있는 분수가 끝나지 않습니다. 0.1 = 1/10 = 1/(2·5)이므로 2진수로 쓰면 무한 반복 소수입니다.

$$0.1_{10} = 0.0\,0011\,0011\,0011\,0011\ldots_2 = \sum_{k=1}^{\infty}\left(2^{-4k}+2^{-4k-1}\right)$$

'0011'이 영원히 반복된다. 유한한 비트에 담으려면 어디선가 잘라야(또는 반올림해야) 하고, 그 순간 저장값은 0.1이 아니다.

컴퓨터의 실수형은 대부분 IEEE 754 부동소수점(Floating-point)입니다. 과학적 표기법 $6.02\times10^{23}$의 2진 버전으로, 64비트 배정밀도(double, binary64)는 부호 1비트, 지수 11비트, 가수 52비트로 나뉩니다.

그림 10-1. IEEE 754 배정밀도(binary64)의 비트 배치. 지수는 1023을 더한 '편향(bias)' 형태로 저장한다. 0.1은 $1.6\times2^{-4}$이므로 e = 1019, 가수에는 0.6을 2진으로 52자리까지 쓴 값이 반올림되어 들어간다.

아래 시뮬레이터에 아무 수나 넣어 보세요. 입력한 10진수가 아니라 실제로 메모리에 들어간 값을 정확한 10진수로 풀어 보여 줍니다. 2진 유한소수는 언제나 10진 유한소수로 정확히 쓸 수 있으므로(2⁻ᵏ = 5ᵏ/10ᵏ) 이 값은 근사가 아니라 완전한 전개입니다.

SIMULATOR

배정밀도 해부기: 저장된 값은 정확히 얼마인가

수 또는 식 (숫자, + − * / ( ) 와 e 표기)

예시

부호지수 11비트가수 52비트

입력과의 오차—

이웃 간격 (ULP)—

float32로 저장하면—

해볼 것: ① 0.1, 0.2, 0.3을 차례로 넣고 각 저장값이 참값보다 큰지 작은지 본다. 0.1과 0.2는 모두 위로 반올림되어 합이 0.3의 저장값(아래로 반올림)과 한 칸 어긋난다. ② 2⁵³+1을 넣으면 입력과 다른 수(2⁵³)가 저장된다. 정수도 2⁵³을 넘으면 끊긴다. ③ 1e−310은 지수 칸이 모두 0인 비정규수다. 모델: 입력 식은 브라우저의 JavaScript(배정밀도)로 계산한다. '입력과의 오차'는 입력 문자열을 BigInt 유리수로 다시 읽어 저장값과 비교한 값이며, 식(연산이 있는 입력)일 때는 각 단계 반올림이 누적된 결과다.

0.1과 0.2는 각각 저장될 때 참값보다 살짝 크게 반올림되고, 둘을 더한 뒤 결과를 다시 53비트로 반올림하면 0.3의 저장값보다 정확히 한 칸 위에 떨어집니다. 그래서 같지 않습니다. 이 '한 칸'의 크기를 단위로 쓰는데, 이것이 ULP(Unit in the Last Place)입니다.

실무 규칙 부동소수점끼리 ==로 비교하지 않는다. |a − b| ≤ tol·max(|a|,|b|)처럼 상대 허용오차로 비교한다. 돈처럼 10진 단위가 정확해야 하는 값은 정수(원, 센트 단위)나 10진 소수형(Decimal)으로 다룬다.

부동소수점 수직선: 간격은 크기에 비례한다

1.0과 10억에 같은 정밀도가 필요할까

센서 신호는 nV에서 V까지, 거리는 nm에서 km까지 쓰입니다. 고정소수점(소수점 위치 고정)이라면 작은 수는 유효숫자가 거의 없고 큰 수는 넘칩니다. 필요한 것은 '절대 오차'가 아니라 '상대 오차'가 일정한 수 체계입니다.

부동소수점은 지수가 같은 구간 $[2^k, 2^{k+1})$(한 '옥타브', 이진 구간(binade)) 안에 $2^{52}$개의 수를 등간격으로 깝니다. 다음 구간으로 가면 간격이 두 배가 됩니다. 그래서 수직선 위의 점들은 0 근처에 빽빽하고 멀어질수록 듬성듬성하지만, 간격 ÷ 크기는 거의 일정합니다. 그 상대 간격의 상한이 기계 엡실론(Machine Epsilon)입니다.

$$\mathrm{fl}(x) = x(1+\delta),\quad |\delta| \le u = \tfrac12\varepsilon,\qquad \varepsilon = 2^{-p+1}$$

$p$: 가수 정밀도(숨은 비트 포함). 배정밀도 $p=53$이면 $\varepsilon = 2^{-52}\approx2.22\times10^{-16}$, 단정밀도 $p=24$이면 $\varepsilon=2^{-23}\approx1.19\times10^{-7}$. $\varepsilon$은 '1과 1 다음 수의 간격'이고, 가장 가까운 수로 반올림하면 상대 오차는 그 절반 $u$ 이하다. 이 한 줄이 모든 반올림 오차 해석의 출발점이다.

실제 64비트 수직선은 너무 촘촘해 그릴 수 없으니, 지수와 가수 비트를 몇 개만 쓰는 장난감 형식(미니플로트(minifloat))으로 모든 수를 찍어 봅시다.

SIMULATOR

미니플로트의 모든 수를 수직선에

지수 비트 E 가수 비트 M 찍어 볼 값 x 비정규수(subnormal) 포함

양수 개수—

기계 엡실론 ε—

최댓값—

x → fl(x), 상대 오차—

위 줄은 선형 눈금, 아래 줄은 로그 눈금(각 이진 구간이 같은 폭)이다. 해볼 것: ① 위 줄에서 오른쪽으로 갈수록 점 간격이 두 배씩 벌어지는 것, 아래 줄에서는 구간마다 점 수가 $2^M$개로 같은 것을 확인한다. ② x를 옮기며 상대 오차가 항상 ε/2 이하인지 본다(비정규수 영역 제외). ③ 비정규수를 끄면 0 근처에 구멍(언더플로 틈)이 생긴다. 모델: IEEE 754식 형식(편향 $2^{E-1}-1$, 최상위 지수는 ∞·NaN 예약, 가장 가까운 짝수로 반올림). 실제 fp8 E4M3(OCP 규격)는 최상위 지수도 일부 수에 써서 최댓값이 448로 이 모델(240)보다 크다.

가수 비트 M개 형식에서 정수는 $2^{M+1}$까지만 빠짐없이 표현됩니다. 그 위에서는 간격이 2 이상이라 홀수가 사라집니다. 단정밀도(가수 23비트)라면 $2^{24}=16\,777\,216$이고, 그래서 float32에서는 16777216 + 1 == 16777216입니다. float32 누적 카운터나 누적 합이 어느 순간 더 이상 늘지 않는 버그가 이것입니다. 배정밀도에서는 $2^{53}\approx9.007\times10^{15}$이 그 경계입니다.

형식	지수/가수 비트	ε (1 다음 수 간격)	대략 최댓값	주 용도
binary64 (double)	11 / 52	2.2×10⁻¹⁶	1.8×10³⁰⁸	과학 계산, 시뮬레이터, JavaScript의 모든 수
binary32 (float)	8 / 23	1.2×10⁻⁷	3.4×10³⁸	그래픽스, 신호 처리, 신경망 학습의 기준
binary16 (fp16)	5 / 10	9.8×10⁻⁴	65 504	GPU 추론, 이미지 HDR 저장
bfloat16 (bf16)	8 / 7	7.8×10⁻³	3.4×10³⁸	신경망 학습: float32와 같은 범위, 정밀도는 포기
fp8 E4M3 / E5M2	4 / 3, 5 / 2	0.125 / 0.25	448 / 57 344	대형 모델 학습·추론의 행렬 곱

표의 아래쪽으로 갈수록 비트를 아낍니다. 신경망 가중치 하나가 1% 틀려도 수백만 개가 평균되면 결과는 크게 변하지 않지만, 값이 범위를 넘으면(오버플로) 학습이 망가집니다. 그래서 bf16은 가수를 깎고 지수는 float32만큼 남긴 형식입니다. fp8은 더 극단적이어서 텐서마다 배율(scale factor)을 따로 곱해 값을 표현 가능한 범위 한가운데로 옮겨 놓고 씁니다(AIBook 참고). 미니플로트 시뮬레이터에서 E=4, M=3을 놓으면 fp8의 점들이 얼마나 성긴지 보입니다.

패트리어트 계산 다시 하기: 작은 오차가 쌓일 때

이제 첫머리의 사고를 숫자로 재현합니다. GAO 보고서에 따르면 시스템은 0.1초마다 1씩 늘어나는 정수 카운터를 두고, 이를 24비트 고정소수점 레지스터에 담긴 0.1과 곱해 초 단위 시간을 만들었습니다. 저장된 0.1은 2진 소수점 아래 23자리에서 잘린 값이었습니다.

$$0.1 \;\to\; 0.000\,1100\,1100\,1100\,1100\,1100_2 = \frac{\lfloor 0.1\cdot2^{23}\rfloor}{2^{23}} = \frac{838\,860}{8\,388\,608}$$

한 번 셀 때마다 $0.1 - 838\,860/2^{23} \approx 9.54\times10^{-8}$초가 모자란다. 100시간이면 $100\times3600\times10 = 3.6\times10^{6}$번 세므로 누적 오차는 약 0.3433초.

그림 10-2. 레이더는 탐지한 표적의 다음 위치를 예측해 그 주변의 '레인지 게이트'만 살핀다. 시간 계산이 0.34초 어긋나면 예측 위치가 표적의 이동 거리만큼 밀리고, 창 안에 아무것도 없으니 시스템은 추적을 포기한다.

SIMULATOR

0.1초 시계의 누적 오차

연속 가동 시간 0.1을 담는 소수 비트 수 표적 속도 GAO 보고서 표 값 표시

한 번 셀 때 오차—

센 횟수—

누적 시간 오차—

속도 × 오차—

해볼 것: ① 23비트, 100시간에서 누적 오차가 0.3433초인지 확인한다(GAO 표의 점과 겹친다). ② 비트 수를 1씩 늘려 보면 오차가 대략 절반씩 준다. 단 0.1의 2진 전개 '0011' 반복 때문에 잘리는 위치에 따라 들쭉날쭉하다. ③ 8시간이면 오차가 0.03초 남짓이다. 짧게 가동하고 재부팅하면 피할 수 있는 버그였다. 모델: 오차 = (0.1 − 잘린 0.1) × 센 횟수, 거리 = 표적 속도 × 시간 오차의 단순 추정. 1676 m/s는 GAO가 인용한 스커드 속도다. GAO 표의 '레인지 게이트 이동'(100시간에 687 m)은 시스템 내부 계산 방식을 반영한 값이라 이 단순 곱(약 575 m)과 다르다.

교훈은 두 가지입니다. 첫째, 한 번의 오차가 $10^{-7}$이라도 같은 방향으로 백만 번 쌓이면 $10^{-1}$이 됩니다. 무작위 반올림 오차는 서로 상쇄되어 $\sqrt{N}$으로 자라지만, 잘림처럼 한쪽으로 치우친 오차는 $N$으로 자랍니다. 둘째, 같은 시간 값을 일부 루틴은 개선된 변환으로, 일부는 옛 변환으로 계산하면서 두 값의 차이를 쓰는 바람에 오차가 상쇄되지 않았다는 분석도 있습니다(R. Skeel, 1992). 시간처럼 오래 누적되는 양은 정수 틱으로 세고 마지막에 한 번만 변환하는 것이 정석입니다.

상쇄 오차: 비슷한 두 수의 뺄셈

정확한 공식이 틀린 답을 낸다

렌즈 처짐(sag), 빛의 경로차, 작은 각도의 위상차를 계산하다 보면 $1-\cos x$ 같은 식이 자주 나옵니다. $x = 10^{-8}$이면 참값은 $5\times10^{-17}$인데, 배정밀도로 1 - Math.cos(1e-8)를 계산하면 정확히 0이 나옵니다. 수식은 맞는데 답은 100% 틀렸습니다.

원인은 뺄셈 자체가 아니라 그 전 단계입니다. $\cos(10^{-8}) = 1 - 5\times10^{-17}$인데, 1 근처의 간격은 $\varepsilon\approx2.2\times10^{-16}$이라 이 값은 저장되는 순간 1로 반올림됩니다. 뺄셈은 그 반올림 오차를 드러낼 뿐입니다. 일반적으로 $a\approx b$일 때 $a-b$의 상대 오차는

$$\frac{|\Delta(a-b)|}{|a-b|} \lesssim u\,\frac{|a|+|b|}{|a-b|}$$

$a$와 $b$가 각각 상대 오차 $u$를 품고 있으면, 결과의 상대 오차는 $\tfrac{|a|+|b|}{|a-b|}$배로 커진다. 두 수가 앞 $k$자리까지 같으면 유효숫자 약 $k$자리를 잃는다. 이를 파국적 상쇄(Catastrophic Cancellation)라 한다.

치료법은 대수적으로 같은 식으로 바꿔 뺄셈을 없애는 것입니다.

$$1-\cos x = 2\sin^2\frac{x}{2},\qquad x_{1,2}=\frac{-b\mp\sqrt{b^2-4ac}}{2a}\;\Rightarrow\; x_1 = \frac{-b-\operatorname{sgn}(b)\sqrt{b^2-4ac}}{2a},\;\; x_2 = \frac{c}{a\,x_1}$$

근의 공식: $b^2\gg4ac$이면 $\sqrt{b^2-4ac}\approx|b|$라 $-b+\sqrt{\cdot}$ 쪽 근에서 상쇄가 일어난다. 부호가 같은 두 수를 더하는 쪽 근을 먼저 구하고, 다른 근은 근과 계수의 관계 $x_1x_2=c/a$로 얻는다. 이 밖에 $\sqrt{x+1}-\sqrt{x} = 1/(\sqrt{x+1}+\sqrt{x})$, 표준 라이브러리의 expm1, log1p, hypot도 같은 목적의 도구다.

SIMULATOR

순진한 식 vs 안정적인 식

문제

정밀도

관찰점

참값—

순진한 식—

안정적인 식—

그래프는 상대 오차(로그 눈금)다. 10⁰ = 100% 틀림. 해볼 것: ① 1 − cos x에서 x를 줄이면 순진한 식의 오차가 x⁻²에 비례해 커지다가(기울기 −2) 결국 100%(결과 0)에 붙는다. ② float32로 바꾸면 같은 일이 x ≈ 10⁻³에서 이미 일어난다. ③ 근의 공식(x² + bx + 1 = 0의 작은 근)에서 b가 커질수록 순진한 식이 무너지고, 안정적인 식은 끝까지 ε 수준을 유지한다. 모델: 참값은 배정밀도의 안정적인 식(float32 모드에서는 이것이 충분히 정확한 기준)으로 계산. float32는 각 연산 결과를 Math.fround로 반올림해 흉내 낸다. 오차가 0(완전 일치)이면 그래프 맨 아래에 붙여 그린다.

상쇄는 어디에나 숨어 있다 분산을 $\overline{x^2}-\bar{x}^2$로 계산하면(큰 평균 위의 작은 산포) 음수 분산이 나오기도 한다. 웰퍼드(Welford) 알고리즘처럼 평균을 먼저 빼는 방식으로 계산해야 한다. 이미지 센서의 다크 프레임 차감, 큰 오프셋 위의 작은 신호 측정도 같은 구조다.

수치 미분의 U자 곡선: h를 줄이면 더 정확해질까

시뮬레이터에 기울기가 필요하다

최적화(9장)는 목적함수의 기울기가 필요하지만, 상용 시뮬레이터의 출력은 공식이 아니라 숫자입니다. 가장 쉬운 방법은 $f'(x)\approx\frac{f(x+h)-f(x)}{h}$입니다. 극한의 정의대로라면 h를 작게 할수록 정확해져야 합니다. 정말 그럴까요?

테일러 전개(1장)로 오차를 쓰면 두 항이 싸웁니다.

$$\left|\frac{\mathrm{fl}(f(x+h))-\mathrm{fl}(f(x))}{h}-f'(x)\right| \;\lesssim\; \underbrace{\frac{|f''|}{2}\,h}_{\text{절단 오차}} \;+\; \underbrace{\frac{2u|f|}{h}}_{\text{반올림 오차}}$$

절단 오차(truncation error)는 테일러 급수를 1차에서 자른 대가로 h에 비례해 줄고, 반올림 오차(round-off error)는 거의 같은 두 함수값의 상쇄를 h로 나누므로 1/h로 커진다. 합을 최소화하면 $h^*\approx2\sqrt{u|f|/|f''|}\sim\sqrt{\varepsilon}\approx10^{-8}$, 그때 오차 $\sim\sqrt{\varepsilon}$: 배정밀도의 16자리 중 8자리만 건진다.

중심 차분 $\frac{f(x+h)-f(x-h)}{2h}$은 짝수 차수 항이 상쇄되어 절단 오차가 $\frac{|f'''|}{6}h^2$입니다. 반올림 항 $u|f|/h$과 균형을 맞추면 $h^*\sim\sqrt[3]{\varepsilon}\approx6\times10^{-6}$, 오차 $\sim\varepsilon^{2/3}\approx10^{-11}$입니다.

SIMULATOR

미분 오차 vs 간격 h (로그–로그)

함수 (x = 1에서)

정밀도

h = 10^ 이론 오차 상한(점선)

전진 차분 오차—

중심 차분 오차—

최적 h (전진 / 중심)—

해볼 것: ① 오른쪽(큰 h)에서 전진 차분은 기울기 1, 중심 차분은 기울기 2로 내려간다(절단 오차 ∝ h, h²). ② 왼쪽에서는 둘 다 기울기 −1로 다시 올라가며 들쭉날쭉하다. 반올림 오차는 무작위처럼 보인다. ③ x³ + 1000은 함수값이 커서(|f| 큼) 반올림 항이 커지고 최적 h가 오른쪽으로 옮겨 간다. float32에서는 최적 h가 10⁻³ 근처, 최선의 오차도 훨씬 크다. 모델: 오차 = |차분 근사 − 해석적 도함수|, 0이면 10⁻¹⁹에 표시. 점선은 위 공식의 상한 $\tfrac{|f''|}{2}h + \tfrac{2u|f|}{h}$, $\tfrac{|f'''|}{6}h^2 + \tfrac{u|f|}{h}$.

오차 없는 미분 반올림 문제를 피하는 길도 있다. 복소 스텝 미분 $f'(x)\approx\mathrm{Im}\,f(x+ih)/h$(2장의 복소수)는 뺄셈이 없어 $h=10^{-200}$도 쓸 수 있다. 자동 미분(Automatic Differentiation)은 프로그램의 각 연산에 연쇄 법칙을 적용해 기계 정밀도로 정확한 기울기를 준다. 딥러닝 프레임워크의 역전파가 바로 역방향 자동 미분이다.

유한 차분: 미분방정식을 격자 위의 덧셈으로

해석해가 없는 방정식

칩의 열 분포, 막대의 진동, 공정 중 도펀트 확산은 모두 편미분방정식입니다. 경계 모양이 조금만 복잡해도 손으로 푼 해는 없습니다. 컴퓨터는 덧셈과 곱셈밖에 못 하는데, 미분을 어떻게 시킬까요?

답은 앞 절의 차분을 공간과 시간 양쪽에 쓰는 것입니다. 연속 함수 $u(x,t)$를 간격 $\Delta x,\Delta t$의 격자점 값 $u_i^n = u(i\Delta x, n\Delta t)$으로 바꾸고(6장의 표본화), 도함수를 이웃 값의 가중합으로 근사합니다. 어떤 이웃을 어떤 가중치로 쓰는지 그린 그림을 스텐실(stencil)이라 합니다.

그림 10-3. 유한 차분 스텐실. 가중치 (1, −2, 1)은 $u(x\pm\Delta x)$를 테일러 전개해 더하면 1계 항이 상쇄되고 $\Delta x^2 u''$가 남는다는 데서 나온다. 오른쪽은 파동 방정식에서 시간 n−1, n의 값 4개로 n+1의 값 1개를 만드는 스텐실이다.

파동 방정식 $u_{tt} = c^2u_{xx}$의 양변에 2계 중심 차분을 쓰고 미래 값을 풀면 갱신식이 나옵니다.

$$u_i^{n+1} = 2u_i^n - u_i^{n-1} + C^2\left(u_{i+1}^n - 2u_i^n + u_{i-1}^n\right),\qquad C = \frac{c\,\Delta t}{\Delta x}$$

$C$: 쿠랑 수(Courant number). 한 시간 스텝 동안 파동이 몇 칸 가는지를 뜻한다. 프로그램은 이 한 줄을 모든 i에 대해 반복할 뿐이다.

여기에 함정이 있습니다. 스텐실은 한 스텝에 이웃 한 칸까지만 봅니다. 그런데 실제 파동은 $\Delta t$ 동안 $c\Delta t$만큼 갑니다. $c\Delta t > \Delta x$이면 실제로는 영향을 받아야 할 정보가 격자 알고리즘의 '시야' 밖에 있게 됩니다. 이것이 1928년 쿠랑·프리드리히스·레비가 밝힌 CFL 조건(Courant–Friedrichs–Lewy condition)입니다.

그림 10-4. 한 점(주황)의 값이 과거의 어느 범위에 의존하는가. 실제 해는 기울기 $\pm c$인 원뿔 안의 초기값에, 격자 해는 기울기 $\pm\Delta x/\Delta t$ 안의 값에 의존한다. 격자의 영역이 실제 영역을 덮어야 수렴할 수 있다(필요조건).

SIMULATOR

1D 파동 방정식과 CFL 조건

쿠랑 수 C = cΔt/Δx

초기 모양

초기값에 10⁻¹² 잡음 섞기

시간 스텝—

max |u|—

최고 주파수 모드 증폭률—

해볼 것: ① C = 0.9에서 펄스가 둘로 갈라져 양 끝에서 뒤집혀 반사되는 것을 본다. ② C를 1.02로 올리고 실행하면 처음에는 멀쩡하다가 수십~백여 스텝 뒤 톱니 모양 진동이 폭발한다. 반올림 수준(10⁻¹⁶)의 고주파 성분이 매 스텝 증폭된 것이다. ③ C = 1.00은 1차원 파동에서 차분 해가 정확해지는 '마법의 시간 간격'이다. 0.5로 낮추면 안정하지만 뾰족한 삼각 모양 뒤에 잔물결(수치 분산)이 생긴다. 모델: 격자 200칸, 양 끝 고정(u = 0), 정지 상태에서 출발. 증폭률은 폰 노이만 분석으로 구한 파수 π/Δx 모드의 |λ|.

시뮬레이터의 폭발은 폰 노이만 안정성 분석(von Neumann stability analysis)으로 정확히 예측됩니다. 오차를 푸리에 모드 $e^{ik i\Delta x}$(5장)로 분해해 한 스텝에 몇 배가 되는지(증폭률 $\lambda$) 계산하면, 위 스킴은 $\lambda^2 - 2(1-2C^2\sin^2\tfrac{k\Delta x}{2})\lambda + 1 = 0$을 만족합니다. $C\le1$이면 모든 k에서 $|\lambda|=1$이고, $C>1$이면 가장 짧은 파장($k\Delta x=\pi$)에서 $|\lambda|>1$인 근이 생깁니다. 반올림 오차는 모든 주파수 성분을 조금씩 갖고 있으므로, 증폭되는 모드가 하나라도 있으면 반드시 폭발합니다.

열 방정식은 더 엄격하다 확산 방정식 $u_t = D u_{xx}$을 전진 오일러(명시적 방법)로 풀면 안정 조건이 $D\Delta t/\Delta x^2\le\tfrac12$이다. 격자를 2배 촘촘히 하면 시간 간격은 4배 줄여야 한다. 그래서 확산·열 해석(TCAD의 도펀트 확산, 4장의 열 방정식)은 조건 없이 안정한 암시적 방법(크랭크–니콜슨 등)으로 매 스텝 연립방정식을 푸는 경우가 많다.

FDTD: 맥스웰 방정식을 격자에서 직접 풀기

1 µm 픽셀 속의 빛

스마트폰 이미지 센서의 픽셀은 한 변이 1 µm 안팎으로, 가시광 파장과 비슷합니다. 이 크기에서는 광선 추적이 맞지 않습니다. 마이크로렌즈, 컬러 필터, 금속 격벽, 실리콘 속 깊은 트렌치 격리(DTI)를 지나는 빛의 회절과 간섭을 계산하려면 맥스웰 방정식을 그대로 풀어야 합니다.

1966년 케인 이(Kane Yee)는 유한 차분을 맥스웰 방정식에 쓰는 우아한 방법을 냈습니다. 1차원($x$로 진행하는 $E_z, H_y$)으로 쓰면 두 방정식이 서로를 갱신합니다.

$$\varepsilon\frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x},\qquad \mu\frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x}$$

(부호는 장의 방향 규약에 따라 다르게 쓰기도 한다.) 핵심은 E의 시간 변화가 H의 공간 변화로, H의 시간 변화가 E의 공간 변화로 정해진다는 것이다. Yee는 E와 H를 공간에서 반 칸, 시간에서 반 스텝 엇갈려 놓았다. 그러면 모든 미분이 자동으로 정확도 2차의 중심 차분이 된다.

그림 10-5. 1차원 Yee 격자. 원(E)은 정수 격자점과 정수 시간에, 사각형(H)은 반 칸 위치와 반 스텝 시간에 산다. 3차원에서는 정육면체 셀의 모서리에 E의 세 성분, 면 중심에 H의 세 성분을 놓아 회전(curl) 연산이 저절로 맞물린다.

FDTD(Finite-Difference Time-Domain)도 앞 절의 파동 방정식과 같은 CFL 조건을 따릅니다. 3차원 균일 격자라면 $c\Delta t\le\Delta x/\sqrt3$입니다. 또 계산 영역은 유한하므로 바깥으로 나가는 파동이 경계에서 되돌아오지 않게 해야 합니다. 아래 시뮬레이터는 가장 단순한 흡수 경계(1차 무어(Mur) 조건의 특수한 경우)를 쓰고, 실무에서는 경계에 인공 손실층을 두는 PML(Perfectly Matched Layer)을 씁니다.

SIMULATOR

1D FDTD: 유전체 판에 부딪히는 펄스

판의 굴절률 n 판 두께 흡수 경계(끄면 완전 도체 벽) H 함께 표시

시간 스텝—

프레넬 반사 계수 r—

FDTD로 잰 r (첫 면)—

해볼 것: ① n = 2에서 첫 면 반사 펄스가 뒤집혀(음수) 돌아오고, 측정한 r이 이론값 (1−n)/(1+n) = −1/3에 맞는지 본다. ② 판 안에서 펄스가 n배 느려지고 짧아지는 것, 뒷면에서 반사된 두 번째 펄스(부호 +)를 찾는다. ③ 흡수 경계를 끄면 양 끝에서 펄스가 뒤집혀 영원히 돌아다닌다. 모델: 정규화 단위의 1D FDTD(Sullivan 형식), 쿠랑 수 0.5, 무손실 비자성 유전체(ε = n²), 소스 위치에 가우시안 펄스를 더하는 소프트 소스. r은 소스와 판 사이 탐침에서 입사 펄스 최댓값 대비 첫 반사 펄스의 극값.

FDTD의 장점은 한 번의 시간 영역 계산에 넓은 대역의 펄스를 넣고 푸리에 변환하면 모든 파장의 응답이 한꺼번에 나온다는 것, 구조가 아무리 복잡해도 격자에 재료값만 칠하면 된다는 것입니다. 대가는 계산량입니다. 파장당 10~20칸 이상이 필요하고 3차원이면 격자 수가 세제곱으로 늘어납니다. 그래서 다음 절의 방법이 함께 쓰입니다.

RCWA: 주기 구조는 푸리에 고조파로 푼다

같은 픽셀이 수천만 개 반복된다

이미지 센서의 픽셀 배열, 반도체 노광 마스크의 라인·스페이스 패턴, 회절 격자, 메타표면은 같은 단위 셀이 주기적으로 반복됩니다. 전체를 FDTD 격자로 덮을 필요가 있을까요? 한 주기만 알면 나머지는 같지 않을까요?

주기 $\Lambda$인 구조에 평면파가 들어오면 산란된 장도 같은 주기성(블로흐–플로케 조건)을 가집니다. 그래서 장을 이산적인 방향들의 합으로 쓸 수 있고(5장의 푸리에 급수), 각 방향이 회절 차수(diffraction order)입니다. 방향은 다음 격자 방정식(grating equation)으로 정해집니다.

$$n_{\text{out}}\sin\theta_m = n_{\text{in}}\sin\theta_i + m\frac{\lambda}{\Lambda},\qquad m = 0,\pm1,\pm2,\ldots$$

$\lambda$: 진공 파장, $\Lambda$: 주기. 격자는 입사파의 가로 파수에 $2\pi m/\Lambda$만 더할 수 있다. 우변의 절댓값이 $n_{\text{out}}$보다 크면 그 차수는 퍼져 나가지 못하고 표면에 붙어 지수적으로 사라지는 소멸파(evanescent wave)가 된다.

SIMULATOR

회절 차수의 방향

파장 λ 주기 Λ 입사각 θᵢ 아래 매질 굴절률

반사 전파 차수—

투과 전파 차수—

λ/Λ—

해볼 것: ① Λ를 줄여 λ/Λ > 2가 되면(수직 입사 기준 Λ < λ/2) 반사와 투과 모두 0차만 남는다. 이것이 '서브파장 구조'로, 반사 방지 나방눈 구조와 메타렌즈의 조건이다. ② 아래 매질을 실리콘(n ≈ 4)으로 바꾸면 투과 쪽에는 공기 쪽보다 많은 차수가 산다. 1 µm 픽셀 배열이 실리콘 안으로 빛을 여러 방향으로 흩뜨리는 이유다. ③ 파장을 바꾸면 ±1차의 각도가 바뀐다. 분광기의 원리다. 모델: 위 매질은 공기(n = 1), 각 차수의 방향만 계산하며 세기(회절 효율)는 계산하지 않는다. 세기는 아래에 설명할 RCWA 같은 엄밀 해석이 필요하다.

격자 방정식은 방향만 알려 줍니다. 각 차수에 에너지가 얼마씩 가는지는 격자 안의 장을 풀어야 합니다. RCWA(Rigorous Coupled-Wave Analysis, 엄밀 결합파 해석)는 이렇게 합니다.

그림 10-6. RCWA의 흐름. z 방향으로 굴절률이 일정한 층들로 근사하고, 각 층에서 x 방향은 푸리에 급수로 전개한다. 층 안의 해는 행렬 고유값 문제(3장)로 정확히 풀리고, 층과 층은 경계 조건으로 이어 붙인다. 결과로 각 회절 차수의 복소 진폭이 나온다.

근사는 단 하나, 푸리에 고조파를 유한 개 $2N+1$에서 자른다는 것입니다. 고조파를 늘리면 결과가 수렴하고, 실무에서는 N을 키워 가며 회절 효율이 더 변하지 않는 지점을 찾습니다(수렴 시험). 날카로운 금속 모서리처럼 유전율 대비가 큰 구조는 깁스 현상(5장) 때문에 수렴이 느려서, 리(L. Li, 1996)가 정리한 '푸리에 인수분해 규칙'으로 행렬을 만드는 방식이 표준이 되었습니다. 2차원 주기(픽셀 배열)이면 고조파가 $(2N_x+1)(2N_y+1)$개라 행렬 크기가 빠르게 커지고, 고유값 계산 비용은 그 세제곱에 비례합니다.

	FDTD	RCWA
영역	시간 영역, 공간 격자	주파수 영역, 푸리에 공간(층별)
잘 맞는 구조	임의 모양, 비주기, 비선형	주기 구조, 층으로 나눌 수 있는 단면
한 번의 계산으로	넓은 대역(펄스 + 푸리에 변환)	한 파장·한 입사각
주된 수치 오차	격자 분산, CFL, 경계 반사	고조파 절단, 계단 근사
쓰이는 곳	픽셀 광학, 안테나, 나노포토닉스	노광 마스크·계측(OCD), 격자, 픽셀 배열

조건수: 문제 자체가 오차를 키울 때

측정값을 0.1% 바꿨더니 답이 50% 바뀌었다

공정 계측에서 두 측정값으로 막 두께 두 개를 역산하거나, 색 보정에서 RGB 응답으로 스펙트럼을 추정할 때, 우리는 $A\mathbf{x}=\mathbf{b}$를 풉니다. 계산은 완벽하게 했는데 측정 잡음이 조금만 바뀌어도 답이 크게 흔들린다면, 잘못은 알고리즘이 아니라 문제에 있습니다.

2×2 연립방정식은 두 직선의 교점입니다. 측정값 $\mathbf{b}$의 오차는 직선을 평행 이동시킵니다. 두 직선이 거의 평행하면 조금만 옮겨도 교점은 멀리 미끄러집니다. 이 증폭률의 최악값이 조건수(Condition Number)입니다.

$$\frac{\|\Delta\mathbf{x}\|}{\|\mathbf{x}\|} \;\le\; \kappa(A)\,\frac{\|\Delta\mathbf{b}\|}{\|\mathbf{b}\|},\qquad \kappa(A) = \|A\|\,\|A^{-1}\| = \frac{\sigma_{\max}}{\sigma_{\min}}$$

$\sigma$: 특이값(3장의 SVD). $A$가 공간을 어떤 방향으로는 $\sigma_{\max}$배, 어떤 방향으로는 $\sigma_{\min}$배 늘린다면, 역산은 그 반대로 $1/\sigma_{\min}$배 키운다. 경험 법칙: 조건수가 $10^k$이면 답에서 유효숫자 약 $k$자리를 잃는다. 배정밀도로 $\kappa\approx10^{16}$이면 답의 자릿수가 하나도 믿을 수 없다.

SIMULATOR

거의 평행한 두 직선의 교점

측정 오차 δ (직선 이동 폭) 무작위 오차로 푼 교점 200개

두 직선 사이 각—

조건수 κ—

교점 최대 이동 / δ—

점 네 개를 끌어 두 직선의 방향과 위치를 바꾼다. 해볼 것: ① 두 직선을 직교로 두면 교점이 흔들리는 영역(평행사변형)이 작은 마름모, κ = 1. ② 거의 평행하게 하면 영역이 직선 방향으로 길게 늘어나고 κ가 수십~수백으로 커진다. ③ 각을 절반으로 줄이면 교점 이동이 약 두 배가 되는 것을 확인한다(작은 각 θ에서 κ ≈ 2/θ). 모델: 각 직선을 $\mathbf{n}_i\cdot\mathbf{x} = c_i$(단위 법선)로 쓰고 $c_i$에 ±δ 균일 오차. κ는 법선을 행으로 하는 행렬의 2-노름 조건수 $\sqrt{(1+|\cos\theta|)/(1-|\cos\theta|)}$.

조건수는 알고리즘이 아니라 문제의 성질입니다. 어떤 알고리즘도 $\kappa\cdot u$보다 정확한 답을 보장할 수 없습니다. 좋은 알고리즘(후방 안정(backward stable), 예: 부분 피벗팅 가우스 소거, QR 분해)은 '조금 다른 문제의 정확한 답'을 주며, 그 이상은 문제를 바꿔야 얻습니다. 유명한 예는 힐베르트 행렬 $H_{ij}=1/(i+j-1)$로, 5×5에서 $\kappa\approx5\times10^{5}$, 10×10에서 $\kappa\approx1.6\times10^{13}$입니다. 다항식 피팅에서 $1, x, x^2,\ldots$를 기저로 쓰면 이런 행렬이 나오기 때문에, 고차 피팅에는 직교 다항식이나 QR 기반 최소제곱을 씁니다.

조건수를 낮추는 세 가지 방법 ① 측정을 다시 설계한다(두 직선이 덜 평행하도록, 예: 계측 파장·각도를 다양하게). ② 변수의 단위·스케일을 맞춘다(nm와 m를 한 벡터에 섞지 않는다). ③ 정칙화(티호노프, 9장): 작은 특이값 방향의 해를 억제해 잡음 증폭을 막고, 대신 약간의 편향을 받아들인다.

이 도구가 쓰이는 곳

시리즈의 다른 책에서 이 장의 도구가 등장하는 곳입니다. 시뮬레이터의 결과를 볼 때 "격자는 충분히 촘촘한가, 시간 간격은 안정 조건 안인가, 고조파는 수렴했는가, 문제의 조건수는 괜찮은가"를 먼저 묻는 습관이 이 장의 목표입니다.

TCADBOOK소자·공정 시뮬레이션
확산·포아송 방정식의 유한 차분/유한 요소, 암시적 시간 적분, 뉴턴 반복과 조건수 OPTICSBOOK파동 광학 시뮬레이션
FDTD, RCWA, 회절 차수와 격자 방정식, 박막 전달 행렬 SENSORBOOK픽셀 광학과 신호 처리
마이크로렌즈·DTI 구조의 FDTD 해석, 다크 프레임 차감의 상쇄, ADC 양자화 LITHOBOOK마스크와 계측
마스크 3D 효과의 RCWA, OCD 계측의 역문제와 조건수 COMPUTERBOOK컴퓨터 산술
IEEE 754 하드웨어, FPU와 융합 곱셈–덧셈(FMA), 정수 오버플로 AIBOOK저정밀도 학습·추론
bf16·fp8, 손실 스케일링, 혼합 정밀도, 자동 미분과 역전파 SOCBOOK고정소수점 DSP
비트 폭 결정, 잘림과 반올림, 누적기 오버플로 MONEYBOOK돈의 계산
금액은 왜 부동소수점이 아닌 정수·10진 소수로 다루는가, 복리의 누적 반올림

핵심 정리

부동소수점은 $(-1)^s\times1.f\times2^{e-\text{bias}}$. 0.1처럼 분모에 2 이외의 소인수가 있는 수는 2진수로 끝나지 않아 저장 순간 반올림된다. 그래서 0.1 + 0.2 ≠ 0.3.
표현 가능한 수의 간격은 크기에 비례하고 상대 오차는 $u=\varepsilon/2$ 이하다(배정밀도 $\varepsilon=2^{-52}\approx2.2\times10^{-16}$). 가수 비트 수가 정밀도를, 지수 비트 수가 범위를 정한다. bf16은 범위를, fp16은 정밀도를 택했다.
한쪽으로 치우친 작은 오차는 횟수에 비례해 쌓인다. 패트리어트의 0.1초 시계는 틱당 $9.5\times10^{-8}$초 오차가 100시간 동안 0.34초가 되었다.
비슷한 두 수의 뺄셈은 앞자리가 같은 만큼 유효숫자를 잃는다(상쇄). 식을 대수적으로 바꿔 뺄셈을 없앤다: $2\sin^2(x/2)$, 근과 계수의 관계, expm1·log1p.
수치 미분은 절단 오차(∝ h 또는 h²)와 반올림 오차(∝ 1/h)의 합이라 U자다. 최적 간격은 전진 차분 $\sim\sqrt\varepsilon$, 중심 차분 $\sim\sqrt[3]\varepsilon$.
명시적 유한 차분은 CFL 조건 $c\Delta t/\Delta x\le1$(3차원 FDTD는 $1/\sqrt3$)을 어기면 반올림 수준의 고주파가 매 스텝 증폭되어 폭발한다.
FDTD는 E와 H를 반 칸·반 스텝 엇갈린 Yee 격자에서 번갈아 갱신하고, RCWA는 주기 구조를 층별 푸리에 고조파의 고유값 문제로 푼다. 회절 차수 방향은 $\sin\theta_m=\sin\theta_i+m\lambda/\Lambda$.
조건수 $\kappa=\sigma_{\max}/\sigma_{\min}$는 문제 자체의 오차 증폭률이다. $\kappa\approx10^k$이면 약 k자리를 잃는다.

확인 퀴즈

1. 다음 중 IEEE 754 배정밀도로 정확히 저장되는 수는?

0.375 = 3/8 = 0.011₂로 분모가 2의 거듭제곱이라 유한 2진 소수다. 0.1 = 1/10과 0.2 = 1/5는 분모에 5가, 1/3은 3이 있어 2진수로 끝나지 않는다.

2. float32로 1부터 1씩 계속 더하는 카운터가 있다. 결국 값이 더 이상 늘지 않는 지점은?

float32는 가수 23비트 + 숨은 비트 1 = 24비트 정밀도다. 2²⁴ 이상에서는 이웃 수 간격이 2라서 2²⁴ + 1은 짝수 쪽(2²⁴)으로 반올림되고 카운터가 멈춘다. 3.4×10³⁸은 최댓값(오버플로 경계)이다.

3. 작은 x에서 $\sqrt{1+x}-1$을 정확하게 계산하려면?

분자를 유리화하면 $\sqrt{1+x}-1 = x/(\sqrt{1+x}+1)$이 되어 비슷한 두 수의 뺄셈이 사라진다. 정밀도를 높이면 상쇄가 일어나는 x의 범위가 작은 쪽으로 밀릴 뿐, 충분히 작은 x에서는 같은 문제가 생긴다.

4. 배정밀도에서 중심 차분 $\frac{f(x+h)-f(x-h)}{2h}$로 미분할 때 h를 $10^{-12}$로 잡았다. 결과는?

중심 차분의 최적 h는 $\sqrt[3]{\varepsilon}\approx6\times10^{-6}$ 근처다. h = 10⁻¹²에서는 절단 오차(∝ h² ≈ 10⁻²⁴)는 무시할 만하지만 반올림 오차 ε/h ≈ 10⁻⁴가 지배한다. U자 곡선의 왼쪽 벽이다.

5. 1D 파동 시뮬레이션이 Δx = 1 mm, c = 340 m/s(공기 중 음속)에서 안정했다. 격자를 Δx = 0.5 mm로 촘촘히 하면 명시적 리프프로그 스킴의 Δt는?

CFL 조건 cΔt/Δx ≤ 1에서 Δt의 상한은 Δx에 비례한다. Δx를 절반으로 하면 Δt도 절반으로 줄여야 한다(Δx = 0.5 mm이면 Δt ≤ 약 1.47 µs). 1/4로 줄여야 하는 것은 확산 방정식(Δt ∝ Δx²)의 명시적 방법이다.

6. 주기 Λ = 400 nm인 격자에 λ = 550 nm 빛이 공기 중에서 수직으로 입사한다. 공기 쪽으로 반사되는 전파 차수는?

sin θ₁ = λ/Λ = 1.375 > 1이므로 ±1차는 소멸파가 된다. 공기 쪽으로는 0차(거울 반사)만 퍼져 나간다. 다만 아래 매질이 n = 1.5 유리라면 투과 쪽에서는 1.375 < 1.5라서 ±1차가 산다(약 66°).