테마 선택

📚 지식 노트

AI·하드웨어·주변기기·게임·컴퓨터 시스템·이론 컴퓨터 과학·수학·개발과 운영·보안·응용 분야에 대해 알아 두면 좋은 내용을 정리한 노트다. 위에서 분야를 고르고 아래에서 노트를 고른다.

계산 과학

풀 수 없는 방정식을 컴퓨터로 근사해 푸는 일 — 날씨 예보, 다리의 응력, 별의 궤도 — 의 바탕에는 유한한 비트로 수를 적는 방식과, 근사가 얼마나 틀리고 언제 폭발하는지에 대한 이해가 있다. 부동소수점에서 시작해 근 찾기·선형 시스템·미분 방정식·몬테카를로·열 방정식을 작은 예로 직접 계산하며 따라간다.

목차 (13)
  1. 컴퓨터의 수: IEEE 754 부동소수점
  2. 오차: 반올림, 절단, 그리고 조건수
  3. 방정식의 근: 이분법
  4. 선형 시스템: 직접법과 반복법
  5. 미분 방정식 시뮬레이션: 오일러와 RK4
  6. 안정성과 시간 간격: 강성 문제
  7. 몬테카를로: 무작위로 적분하기
  8. 편미분 방정식: 열 방정식과 유한 차분
  9. 유한 요소법
  10. N체 시뮬레이션과 날씨·기후 모델
  11. 검증과 확인(V&V)
  12. 도구: NumPy, SciPy, Julia, Fortran
  13. 핵심 자료

컴퓨터의 수: IEEE 754 부동소수점

컴퓨터는 실수를 유한한 비트로 적는다. 거의 모든 CPU·GPU·언어가 따르는 규칙이 IEEE 754다. 1985년에 처음 나왔고 2008년, 2019년에 개정됐다. 과학 계산의 기본인 64비트 binary64(C의 double, Python의 float, JavaScript의 수)는 부호 1비트, 지수 11비트, 가수 52비트로 ±1.xxx… × 2지수를 적는다. 맨 앞의 1은 늘 1이니 적지 않아서 실제 정밀도는 53비트다.

문제는 10진수로 깔끔한 수가 2진수로는 끝나지 않는다는 것이다. 1/3이 10진수로 0.333…인 것처럼, 1/10은 2진수로 0.000110011…이 끝없이 되풀이된다. 그래서 0.1은 53비트에서 반올림된 가장 가까운 수로 저장된다. Python 문서가 보여 주듯 그 값은 정확히 3602879701896397 / 255이고, 0.1보다 아주 조금 크다. 0.1과 0.2가 모두 조금씩 크게 저장되니 더한 결과도 0.3에 가장 가까운 double이 아니라 그 다음 수가 된다. 0.1 + 0.2 == 0.3이 거짓인 이유다. 결과는 언어와 관계없이 같다. IEEE 754가 연산 결과를 "참값을 기본 규칙(가장 가까운 수, 같으면 짝수 쪽)으로 반올림한 것"으로 정해 두었기 때문이다.

1과 그다음으로 큰 double 사이의 간격을 기계 엡실론 ε라 한다. 1 + ε/2가 1과 같아질 때까지 ε를 반으로 줄여 보면 ε = 2−52 ≈ 2.22e-16이 나온다. 반올림 한 번의 상대 오차는 ε/2를 넘지 않는다. 그래서 double로는 대략 유효 숫자 15~16자리를 믿을 수 있다. 이웃한 double 사이의 간격은 수의 크기에 비례해 커진다. 253(약 9×1015)을 넘으면 간격이 2가 되어 홀수 정수를 정확히 적지 못한다.

0.1의 binary64 비트 구성과 0.1, 0.2, 0.1 더하기 0.2, 0.3이 실제로 저장한 값지수 11비트가수 52비트(앞의 1은 저장하지 않음)↓ 부호 1비트부호 0 · 지수 01111111011 = 1019 − 1023 = -41001 1001 1001 1001 1001 1001 1001 1001 1001 1001 1001 1001 10101.1001 1001 …(이진 순환 소수)을 52비트에서 반올림 → 끝이 …1010으로 올라감실제로 저장된 값0.10.10000000000000000555111512312578270211815834045410156250.20.2000000000000000111022302462515654042363166809082031250.1 + 0.20.30000000000000004440892098500626161694526672363281250.30.299999999999999988897769753748434595763683319091796875기계 엡실론 ε = 2⁻⁵² ≈ 2.22 × 10⁻¹⁶: 1과 그다음 double 사이의 간격
0.1의 binary64 표현과 실제로 저장된 값. 0.1 + 0.2의 결과(색칠)와 0.3을 적은 값(빨강)이 마지막 자리에서 다르다. 정확한 10진수 값은 비트에서 직접 풀어 계산했다.

그래서 실수 계산의 결과는 ==로 비교하지 않는다. 두 값의 차이가 허용 오차 안에 드는지를 본다. 허용 오차는 값의 크기에 비례하는 상대 오차와 0 근처를 위한 절대 오차를 함께 쓴다. 값이 0에 가까우면 상대 오차가 아무 뜻이 없어지기 때문이다. Python의 math.isclose는 기본 상대 오차가 10−9, 절대 오차가 0이라, 0과 비교할 때는 절대 오차를 직접 줘야 한다. 정밀도가 낮은 형식(32비트 float, 학습용 BF16)은 ε가 훨씬 커서 같은 계산도 더 일찍 틀어진다. 형식별 비트 배치는 AI 인프라 › 혼합 정밀도에 있다.

오차: 반올림, 절단, 그리고 조건수

계산 과학의 오차는 크게 두 갈래다. 반올림 오차는 앞 절처럼 유한한 비트 때문에 생긴다. 절단 오차는 무한한 과정을 유한하게 끊어서 생긴다. 도함수를 (f(x+h) − f(x))/h로, 적분을 직사각형 몇 개의 합으로, 무한급수를 몇 항으로 바꾸는 순간 생기는 오차다. h를 줄이면 절단 오차는 줄지만 거의 같은 두 수를 빼는 일이 늘어 반올림 오차가 커진다. 그래서 가장 좋은 h는 0이 아니라 둘이 맞서는 중간 어딘가에 있다.

가장 무서운 반올림 오차는 거의 같은 두 수의 뺄셈에서 나온다(자릿수 손실, catastrophic cancellation). 앞자리는 같아서 지워지고, 남는 것은 원래 반올림 오차가 들어 있던 뒷자리뿐이다. (1 − cos x)/x²는 x가 0으로 가면 1/2에 다가가지만, x = 10−8에서 그대로 계산하면 cos x가 1로 반올림돼 결과가 0이 된다. 같은 식을 2 sin²(x/2)/x²로 바꾸면 뺄셈이 사라져 결과가 0.5 그대로다. 이차 방정식 근의 공식에서 b와 √(b² − 4ac)가 비슷할 때 생기는 문제도 같은 원리이고, Goldberg의 글이 식을 바꿔 피하는 방법을 정리했다.

알고리즘과 별개로 문제 자체가 예민할 수도 있다. 입력이 조금 바뀔 때 답이 얼마나 크게 바뀌는지를 재는 값이 조건수다. 거의 평행한 두 직선의 교점을 구하는 연립방정식은 오른쪽 값을 0.0001만 바꿔도 답이 (1, 1)에서 (0, 2)로 뛴다. 이 행렬의 조건수는 약 40,004이다. 대략 "조건수가 10k이면 유효 숫자 k자리를 잃는다"고 본다. 조건수가 큰 문제는 어떤 좋은 알고리즘으로도 정확히 풀 수 없다. 할 수 있는 것은 문제를 다시 세우는 것이다.

1 빼기 cos x의 자릿수 손실 표와 조건수가 큰 연립방정식의 민감도자릿수 손실: 참값은 x → 0에서 0.5x(1 − cos x)/x²2sin²(x/2)/x²10−20.49999583330.499995833310−40.49999999700.499999999610−60.50004445030.500000000010−800.5000000000cos x가 1에 너무 가까워 1 − cos x에서앞자리가 모두 지워진다. 식을 바꾸면뺄셈 없이 같은 값을 구한다조건수: 거의 평행한 두 직선x + y = 2x + 1.0001y = 2.0001답 (1, 1)오른쪽 값을 0.0001만 바꿔x + 1.0001y = 2.0002답 (0, 2)조건수 κ(A) ≈ 40,004입력의 상대 오차가 최대 이만큼 커진다
왼쪽: 같은 값을 두 식으로 계산했다. x가 작아질수록 그대로 쓴 식(빨강)은 자릿수를 잃다가 0이 된다. 오른쪽: 조건수가 큰 2×2 연립방정식. 수치는 모두 double로 직접 계산했다.

덧셈 순서도 결과를 바꾼다. 부동소수점 덧셈은 결합 법칙이 성립하지 않아 (0.1 + 0.2) + 0.3은 0.6000000000000001, 0.1 + (0.2 + 0.3)은 0.6이다. 큰 합계에 작은 값을 하나씩 더하면 작은 값의 아랫자리가 매번 잘려 나가고, 0.1을 열 번 더하면 1이 아니라 0.9999999999999999가 된다. Goldberg의 글은 작은 수부터 더하는 정렬보다 훨씬 싼 방법으로 카한 합산(Kahan summation)을 소개한다. 매 덧셈에서 잘려 나간 부분을 따로 모아 두었다가 다음 덧셈에 돌려 넣는다. NumPy의 sum은 많은 경우 부분 쌍별 합산으로 오차를 줄인다. 같은 이유로 병렬 계산은 스레드 수나 GPU 종류가 바뀌어 더하는 순서가 달라지면 마지막 자리가 달라질 수 있다. 결과를 비트 단위로 재현해야 한다면 합산 순서를 고정한다.

방정식의 근: 이분법

f(x) = 0을 푸는 가장 튼튼한 방법은 이분법이다. f가 연속이고 f(a)와 f(b)의 부호가 다르면 그 사이에 근이 있다(중간값 정리). 가운데 m을 찍어 f(m)의 부호를 보고, 부호가 바뀌는 반쪽만 남긴다. 걸음마다 구간이 정확히 반으로 주니, 폭 1에서 시작해 10−6 아래로 줄이려면 20걸음이면 된다. 미분이 필요 없고 절대 발산하지 않는다는 것이 장점이다.

x² − 2의 근을 이분법으로 좁혀 가는 단계 그림11.52f(x) = x² − 21단계구간 [1, 2]가운데 m = 1.5f(m) = 0.25000→ 왼쪽 반을 남긴다폭 12단계구간 [1, 1.5]가운데 m = 1.25f(m) = -0.43750→ 오른쪽 반을 남긴다폭 0.53단계구간 [1.25, 1.5]가운데 m = 1.375f(m) = -0.10938→ 오른쪽 반을 남긴다폭 0.254단계구간 [1.375, 1.5]가운데 m = 1.4375f(m) = 0.06641→ 왼쪽 반을 남긴다폭 0.1255단계구간 [1.375, 1.4375]가운데 m = 1.40625f(m) = -0.02246→ 오른쪽 반을 남긴다폭 0.06256단계구간 [1.40625, 1.4375]가운데 m = 1.421875f(m) = 0.02173→ 왼쪽 반을 남긴다폭 0.031257단계구간 [1.40625, 1.421875]가운데 m = 1.414063f(m) = -0.00043→ 오른쪽 반을 남긴다폭 0.0156258단계구간 [1.414063, 1.421875]가운데 m = 1.417969f(m) = 0.01064→ 왼쪽 반을 남긴다폭 0.0078138단계 뒤[1.414063, 1.417969]폭 1/256√2 ≈ 1.414214
1 / 9

1단계: [1, 2]의 가운데 1.5에서 f = 0.2500. 부호가 +이니 근은 왼쪽 반에 있다.

이분법. f(a)와 f(b)의 부호가 다르면 그 사이에 근이 있다(중간값 정리). 가운데를 찍어 부호가 바뀌는 반쪽만 남긴다. 한 걸음마다 오차가 정확히 반으로 준다.

대신 느리다. 한 걸음에 이진 자리 하나씩만 얻는다(1차 수렴). 접선을 따라 근으로 뛰는 뉴턴 방법은 근 가까이에서 맞는 자릿수가 걸음마다 두 배로 늘어(2차 수렴) 훨씬 빠르지만, 도함수가 필요하고 출발점이 나쁘면 엉뚱한 곳으로 튄다. 뉴턴 방법의 원리와 예는 미적분 노트의 미분의 쓰임 절에 있다. 실무 라이브러리(SciPy의 brentq 등)는 이분법으로 근을 늘 감싸 두면서 가능할 때 빠른 보간 걸음을 섞는 방식을 쓴다.

선형 시스템: 직접법과 반복법

시뮬레이션의 대부분은 결국 Ax = b를 푸는 일로 바뀐다. 미분 방정식을 격자 위에서 이산화하면 격자점 수만큼의 미지수를 가진 연립방정식이 나온다. 크기가 작거나 중간이면 직접법을 쓴다. 가우스 소거로 A를 삼각 행렬 둘의 곱(LU 분해)으로 바꾸고 앞뒤로 대입한다. 계산량은 n³에 비례한다. 가우스 소거의 원리는 선형대수 노트의 연립방정식 절에 있다. 수치 계산에서는 0에 가까운 수로 나누지 않도록 행을 바꾸는 피벗팅이 꼭 필요하다.

격자가 3차원으로 커지면 미지수가 수백만 개를 넘고 n³은 감당할 수 없다. 다행히 이런 행렬은 거의 모든 칸이 0인 희소 행렬이다(각 격자점은 이웃 몇 개와만 엮인다). 이때는 답을 추측하고 고쳐 나가는 반복법이 낫다. 가장 단순한 야코비 반복은 i번째 식을 xᵢ에 대해 풀고, 오른쪽에는 직전 반복의 값을 넣는다. 대각 원소가 같은 행의 나머지보다 크면(대각 우세) 수렴이 보장된다. 그림의 3×3 예에서 10번 반복하면 잔차가 1.8e-4까지 준다. 실제로는 방금 고친 값을 바로 쓰는 가우스–자이델, 대칭 양의 정부호 행렬에 강한 켤레 기울기법(CG), 격자를 굵게·가늘게 오가는 멀티그리드를 쓴다.

야코비 반복의 반복 값 표와 잔차가 반복마다 줄어드는 로그 그래프4x₁ − x₂ = 2−x₁ + 4x₂ − x₃ = 4−x₂ + 4x₃ = 10반복x₁x₂x₃00.00000.00000.000010.50001.00002.500020.75001.75002.750030.93751.87502.937550.99221.98442.9922100.99991.99992.9999답 (1, 2, 3). 각 줄은 바로 위 줄의 값만 써서 계산한다10110-110-310-5051015반복 횟수잔차 ‖b − Ax‖ (로그)
야코비 반복. (0, 0, 0)에서 출발해 반복할 때마다 답 (1, 2, 3)에 다가가고, 잔차 ‖b − Ax‖가 로그 눈금에서 거의 직선으로(일정한 비율로) 준다.

미분 방정식 시뮬레이션: 오일러와 RK4

운동, 화학 반응, 전염병, 회로의 시간 변화는 y′ = f(t, y) 꼴의 미분 방정식으로 적힌다. 대부분은 식으로 풀리지 않으므로 시간을 h씩 잘라 한 걸음씩 나아간다. 가장 단순한 오일러 방법 yn+1 = yn + h·f(tn, yn)의 원리와 그림은 미적분 노트의 미분 방정식 맛보기 절에 있다. 여기서는 그 오차가 얼마나 빨리 주는지를 비교한다.

고전 4차 룽게–쿠타(RK4)는 한 걸음 안에서 기울기를 네 번(처음, 가운데 두 번, 끝) 재서 가중 평균으로 간다. y′ = y, y(0) = 1을 t = 1까지 풀 때 h = 0.1이면 오일러의 오차는 0.125, RK4는 2.1e-6다. h를 반으로 줄이면 오일러 오차는 약 1/2(1차), RK4는 약 1/16(4차)로 준다. RK4가 한 걸음에 네 배 비싸도, 같은 정확도를 훨씬 큰 걸음으로 얻으니 이득이다. 실무 풀이기(SciPy의 solve_ivp 등)는 차수가 다른 두 방법의 차이로 오차를 추정해, 해가 완만한 곳에서는 걸음을 늘리고 급한 곳에서는 줄이는 적응형 걸음을 쓴다.

걸음 크기에 따른 오일러 방법과 RK4의 오차를 로그 눈금으로 비교한 그래프와 표10010-310-610-90.20.10.050.025걸음 크기 h (오른쪽일수록 작음)t = 1의 오차오일러(1차)RK4(4차)h오일러 오차RK4 오차0.22.3e-13.1e-50.11.2e-12.1e-60.056.5e-21.4e-70.0253.3e-28.7e-9h를 반으로 → 오일러 오차 약 1/2,RK4 오차 약 1/16. RK4는 한 걸음에기울기를 네 번 계산해 비싸지만같은 정확도를 훨씬 큰 h로 얻는다
걸음 크기 h에 따른 t = 1의 오차(로그–로그). 기울기가 곧 차수다. 오일러(빨강)는 1, RK4(색칠)는 4. 수치는 두 방법을 직접 돌려 계산했다.

오래 돌리는 시뮬레이션에서는 정확도 차수보다 무엇을 보존하느냐가 중요할 때가 있다. 행성 궤도를 RK4로 오래 돌리면 에너지가 조금씩 새어 궤도가 서서히 줄거나 커진다. 위치와 속도를 엇갈려 갱신하는 심플렉틱 방법(리프프로그·벨레)은 차수가 낮아도 에너지를 오래 지켜, 분자 동역학과 천체 계산의 표준이다.

안정성과 시간 간격: 강성 문제

정확도와 다른 문제가 안정성이다. y′ = −λy(λ > 0)는 참값 e−λt가 빠르게 0으로 줄어드는 단순한 감쇠다. 오일러를 쓰면 매 걸음 y에 (1 − λh)가 곱해진다. 이 배율의 절댓값이 1을 넘으면, 즉 h > 2/λ이면 해가 부호를 바꾸며 기하급수로 커진다. λ = 50이면 h가 0.04보다 작아야 하고, h = 0.05이면 배율이 -1.5가 되어 걸음마다 부호를 바꾸며 1.5배씩 커진다. 참값은 0인데 계산은 폭발한다.

이런 제약이 정말 문제가 되는 것은 강성(stiff) 문제다. 아주 빠르게 사라지는 성분(빠른 화학 반응, 단단한 용수철)과 천천히 변하는 성분이 섞여 있으면, 우리가 보고 싶은 것은 느린 변화인데 빠른 성분 때문에 양함수 방법은 아주 작은 h를 강요받는다. 해법은 암시적 방법이다. 후방 오일러 yn+1 = yn + h·f(tn+1, yn+1)은 미래 값을 식에 넣어 매 걸음 방정식을 풀어야 해서 비싸지만, 같은 문제에서 배율이 1/(1 + λh)로 h가 아무리 커도 1보다 작다. 회로 시뮬레이터나 화학 반응 계산은 이런 강성용 풀이기를 쓴다.

y' = −50y를 큰 걸음의 오일러로 풀면 폭발하고 작은 걸음이나 암시적 방법은 안정한 그래프-4-2240.10.20.30.40.5시간 t↑↓ 화면 밖으로 폭발y′ = −50y, y(0) = 1참값 e^(−50t)는 금방 0이 된다오일러 h = 0.05배율 1 − λh = -1.50오일러 h = 0.01배율 0.50암시적 오일러 h = 0.05배율 1/(1 + λh) = 0.29|배율| > 1이면 걸음마다 오차가 커진다.h = 0.05로 10걸음 가면y = 57.67
y′ = −50y. 굵은 회색이 참값. 오일러를 h = 0.05로 쓰면(빨강) 배율 −1.5로 부호를 바꾸며 폭발하고, h = 0.01(색칠)이나 같은 h의 암시적 오일러(점선)는 0으로 줄어든다.

몬테카를로: 무작위로 적분하기

몬테카를로 방법은 답을 무작위 표본의 평균으로 추정한다. 한 변이 1인 정사각형에 점을 고르게 던지면 사분원 안에 들어갈 확률이 π/4이니, 안에 든 비율에 4를 곱하면 π가 나온다. 이 노트의 고정된 난수열로 점 100개를 던지면 3.2, 100,000개면 3.13732가 나왔다.

핵심 성질은 오차가 1/√N에 비례한다는 것이다. 표본 평균의 표준 편차가 그렇게 줄기 때문이다. 오차를 1/10로 줄이려면 표본이 100배 필요하니 느리다. 대신 이 속도가 차원과 무관하다. 격자로 적분하면 차원마다 칸 수가 곱해져 10차원만 돼도 감당할 수 없지만, 몬테카를로는 몇 차원이든 같은 비율로 준다. 그래서 입자 물리의 검출기 시뮬레이션, 금융의 옵션 가격, 수많은 빛의 경로를 평균하는 렌더링(컴퓨터 그래픽스 노트의 레이 트레이싱 절)에 쓰인다. 표본을 고르게 퍼뜨리는 준난수, 중요한 곳에 표본을 몰아 주는 중요도 표집이 같은 N으로 오차를 줄이는 기법이다.

정사각형에 던진 무작위 점으로 π를 추정하는 그림과 점 수에 따른 오차 그래프점 400개 중 사분원 안 318개π ≈ 4 × 318/400 = 3.1810010-110-210-310²10³10⁴10⁵던진 점 수 N|추정 − π|표준 오차 ∝ 1/√NN을 100배 → 오차는 약 1/10
몬테카를로 π 추정. 왼쪽은 처음 400개의 점. 오른쪽은 점 수에 따른 실제 오차(색칠)와 표준 오차 4√(p(1−p)/N)(점선). 한 번의 실행이라 들쭉날쭉하지만 1/√N 선을 따라간다.

편미분 방정식: 열 방정식과 유한 차분

공간과 시간에 함께 걸친 현상(열의 확산, 파동, 유체)은 편미분 방정식(PDE)으로 적힌다. 가장 단순한 것이 1차원 열 방정식 ut = α uxx다. "온도가 이웃보다 높으면 내려가고 낮으면 올라간다"는 뜻이다. 유한 차분법은 막대를 Δx 간격의 점으로, 시간을 Δt 걸음으로 자르고 미분을 이웃 값의 차로 바꾼다. 시간은 앞으로 한 걸음, 공간은 가운데 차분을 쓰는 FTCS 방식이면 매 걸음 ui ← ui + r(ui−1 − 2ui + ui+1), r = αΔt/Δx²다.

가운데가 뜨거운 막대의 온도가 유한 차분으로 퍼지는 단계 그림과 r이 클 때의 불안정01왼쪽 끝 0오른쪽 끝 00걸음r = 0.4최고 1.00최저 0.00총열량 3.001걸음r = 0.4최고 1.00최저 0.00총열량 3.002걸음r = 0.4최고 0.68최저 0.00총열량 3.004걸음r = 0.4최고 0.56최저 0.00총열량 2.958걸음r = 0.4최고 0.44최저 0.00총열량 2.6316걸음r = 0.4최고 0.31최저 0.00총열량 1.9332걸음r = 0.4최고 0.16최저 0.00총열량 1.026걸음r = 0.6 (불안정)최고 1.40최저 -0.66총열량 2.38아래 띠: 막대 색이 진할수록 뜨겁다
1 / 8

막대를 칸 10개(격자점 11개)로 나누고 가운데 세 점만 온도 1로 시작한다. 양 끝은 늘 0이다.

1차원 열 방정식 u_t = αu_xx의 유한 차분. 매 걸음 uᵢ ← uᵢ + r(uᵢ₋₁ − 2uᵢ + uᵢ₊₁), r = αΔt/Δx². 이웃보다 높은 점은 내려가고 낮은 점은 올라간다.

앞 절의 안정성 문제가 여기서도 나타난다. 이 방식은 r ≤ 1/2일 때만 안정하다. 그림의 r = 0.4는 매끄럽게 퍼지지만, r = 0.6는 이웃끼리 톱니처럼 엇갈리며 커진다. 공간을 두 배로 잘게 나누면 Δt는 네 배로 작아져야 하니, 세밀한 격자에서는 걸음 수가 빠르게 늘어난다. 그래서 크랭크–니콜슨 같은 암시적 방식으로 매 걸음 (희소한) 연립방정식을 풀어 큰 Δt를 쓰거나, 문제에 맞춰 둘을 고른다. 파동 방정식에서는 같은 역할을 CFL 조건(한 걸음에 정보가 한 칸보다 멀리 가면 안 됨)이 한다.

유한 요소법

유한 차분은 반듯한 격자에서 편하지만, 자동차 차체나 뼈처럼 모양이 복잡하면 격자를 맞추기 어렵다. 유한 요소법(FEM)은 영역을 삼각형·사면체 같은 작은 요소로 나누고(메시), 해를 요소마다 단순한 함수(대개 1차식)를 이어 붙인 것으로 근사한다. 1차원으로 보면 마디마다 그 마디에서 1, 이웃 마디에서 0인 "모자 함수"를 두고, 해를 이 함수들의 가중합으로 적는다. 가중치를 정하는 방법은 방정식을 그대로 맞추는 대신 "적분했을 때 평균적으로 맞게"(약형식) 요구하는 것이다. 결과는 큰 희소 연립방정식이 되어 앞 절의 반복법으로 푼다.

장점은 모양에 맞춰 요소를 자유롭게 배치하고, 응력이 몰리는 모서리·구멍 둘레만 잘게 나눌 수 있다는 것이다. 구조 해석(다리·건물·비행기 날개), 열 전달, 전자기장 해석의 표준 도구다. 메시의 질이 결과를 좌우해서, 지나치게 찌그러진 요소는 정확도를 크게 떨어뜨린다. 메시를 점점 줄여 가며 답이 더 변하지 않는지 확인하는 수렴 시험이 기본 절차다.

삼각형 요소로 나눈 메시와 모자 함수의 합으로 곡선을 근사하는 유한 요소법 개념 그림메시: 영역을 작은 요소로 나눔모서리처럼 값이 급히 변하는 곳은 잘게(색칠)1차원: 모자 함수의 합아래 ∧: 마디마다 하나씩, 이웃 마디에서 0굵은 회색: 참 해 · 꺾은선: 근사
유한 요소법의 두 그림. 왼쪽은 영역을 삼각형 요소로 나눈 메시로, 값이 급히 변하는 모서리 근처(색칠)만 더 잘게 나눴다. 오른쪽은 마디 값을 꺾은선으로 이은 근사와, 그것을 이루는 모자 함수들이다.

N체 시뮬레이션과 날씨·기후 모델

N체 문제는 별·은하·분자처럼 서로 끌어당기는 입자 N개의 움직임을 계산한다. 모든 쌍의 힘을 더하면 N²에 비례하는 일이 필요해, 입자가 수백만 개면 감당할 수 없다. 반스–헛(Barnes–Hut) 같은 트리 방법은 공간을 나무처럼 나누고 멀리 있는 무리를 질량 중심 하나로 묶어 대략 N log N으로 줄인다. 시간 적분에는 앞에서 본 심플렉틱 방법을 써서 에너지를 오래 지킨다.

날씨·기후 모델은 대기(와 바다)를 가로·세로 격자로 나누고, 칸마다 온도·기압·바람·습도를 두어 유체 역학과 열역학 방정식을 시간에 따라 푼다. 구름이나 난류처럼 칸보다 작은 현상은 직접 풀 수 없어 근사식(매개변수화)으로 넣는데, 이것이 모델 간 차이의 큰 원인이다. 대기는 초기값의 작은 차이가 커지는 혼돈계라, 날씨 예보는 초기값을 조금씩 바꾼 여러 번의 실행(앙상블)으로 불확실성까지 함께 내놓는다. 기후 모델은 개별 날씨가 아니라 수십 년 평균의 통계를 묻는다는 점이 다르다. 두 경우 모두 격자를 여러 프로세서에 나눠 맡기고 경계 값을 주고받는 병렬 계산이 필수이며, 그 방법은 병렬 프로그래밍 노트에 있다.

N체 시뮬레이션의 쌍별 힘과 트리 근사, 대기를 격자로 나눈 날씨와 기후 모델 개념 그림N체: 모든 쌍의 힘을 더한다먼 무리는 질량 중심 하나로그대로 하면 N²번, 트리로 묶으면 약 N log N번.시간 적분은 에너지를 지키는 방법을 쓴다날씨·기후: 대기를 격자로칸위땅칸마다 온도·바람·습도를 두고 이웃과주고받는다. 구름처럼 칸보다 작은 일은근사식(매개변수화)으로 넣는다
왼쪽: N체 문제에서 가까운 입자는 쌍마다 힘을 계산하고, 먼 무리는 질량 중심 하나(점선)로 묶는다. 오른쪽: 대기를 칸으로 나눈 격자. 칸 하나가 이웃과 열·운동량·수증기를 주고받는다.

검증과 확인(V&V)

시뮬레이션이 그럴듯한 그림을 내놓는다고 맞는 것은 아니다. 계산 과학은 두 질문을 구별한다. 검증(verification)은 "코드가 방정식을 맞게 푸는가"다. 정확한 해가 알려진 문제로 시험하고, 격자나 걸음을 줄일 때 오차가 이론 차수대로 줄어드는지 본다(앞의 오일러 1차, RK4 4차가 바로 그런 시험이다). 정확한 해가 없으면 답을 먼저 정하고 그에 맞는 원천 항을 거꾸로 만들어 넣는 "만든 해" 방법을 쓴다. 확인(validation)은 "맞는 방정식을 푸는가"다. 결과를 실험·관측과 비교한다. 코드가 완벽해도 모델의 가정(마찰 없음, 2차원 근사)이 틀리면 답은 틀린다.

여기에 불확실성 정량화(UQ)가 더해진다. 입력 매개변수(재료 강도, 초기 조건)에 오차 범위가 있으면 답도 하나의 수가 아니라 범위여야 한다. 입력을 표본으로 바꿔 여러 번 돌리는 몬테카를로가 가장 단순한 방법이다. 재현성도 같은 맥락이다. 난수 씨앗, 라이브러리 버전, 컴파일러 최적화(부동소수점 덧셈 순서가 바뀌면 마지막 자리가 달라진다)를 기록해 두어야 같은 결과를 다시 낼 수 있다.

현실, 수학 모델, 코드 사이의 검증과 확인 관계현실실험·관측 자료수학 모델방정식·가정·매개변수코드이산화한 프로그램검증(verification): 방정식을 맞게 푸는가?맞는 방정식을 푸는가?확인(validation): 결과를 실험과 비교검증 도구: 정확한 해가 있는 문제, 격자를 줄일 때 오차가 이론 차수대로 주는지 확인
검증과 확인. 검증(위)은 모델 → 코드 사이의 간극, 확인(아래)은 코드의 결과 → 현실 사이의 간극을 본다. 둘 다 통과해야 결과를 믿을 수 있다.

도구: NumPy, SciPy, Julia, Fortran

  • NumPy: Python 과학 계산의 바닥이다. N차원 배열과 벡터화(반복문 대신 배열 전체 연산)를 제공하고, 핵심은 최적화된 C 코드다. 파이썬 반복문으로 원소를 하나씩 다루면 느리고, 배열 연산으로 쓰면 빠르다.
  • SciPy: NumPy 위에서 적분, 최적화, 보간, 고윳값, 미분 방정식, 통계, 희소 행렬 같은 표준 알고리즘을 제공한다. 이 노트의 거의 모든 절에 해당하는 함수가 있다.
  • Julia: 과학 계산을 위해 만든 언어로, LLVM으로 기계어로 컴파일돼 반복문을 그대로 써도 빠르다. 다중 디스패치로 같은 알고리즘을 여러 수 타입에 재사용한다.
  • Fortran: 가장 오래된 과학 계산 언어이자 여전히 현역이다. 배열 연산과 병렬(코어레이)을 언어 차원에서 지원하고, 최신 개정판은 Fortran 2023이다. 날씨·기후 모델과 오래된 대형 시뮬레이션 코드의 상당수가 Fortran이다.

층을 내려가 보면 모두 같은 바닥을 공유한다. 행렬 곱과 분해는 BLAS와 LAPACK이 하고, LAPACK은 Fortran 90으로 짜여 있으며 계산을 최대한 BLAS에 맡긴다. 그래서 CPU 제조사가 BLAS를 최적화하면 NumPy·SciPy·Julia·MATLAB이 함께 빨라진다. 그 아래에서는 SIMD·멀티코어·GPU·여러 노드(MPI)가 일을 나눈다. 직접 수치 알고리즘을 짜기 전에 검증된 라이브러리 함수가 있는지 찾는 것이 첫 번째 원칙이다.

사용자 코드, 배열 라이브러리, BLAS와 LAPACK, 병렬 하드웨어로 쌓인 계산 과학 도구의 층내 코드Python 스크립트 · Jupyter 노트북 · Julia배열·알고리즘 라이브러리NumPy(배열) · SciPy(적분·최적화·희소 행렬) · Julia 패키지선형대수 핵심BLAS · LAPACK — Fortran·C로 짜고 CPU마다 최적화병렬 하드웨어SIMD · 멀티코어(OpenMP) · GPU · 클러스터(MPI)
계산 과학 도구의 층. 위의 편한 언어가 아래의 최적화된 선형대수 라이브러리와 병렬 하드웨어를 부른다.

핵심 자료

여러 코어·GPU·노드로 계산을 나누는 방법은 병렬 프로그래밍 노트, 행렬과 연립방정식의 원리는 선형대수 노트에 있다.

자료종류요지링크
What Every Computer Scientist Should Know About Floating-Point Arithmetic논문Goldberg(1991). 반올림 오차, ulp와 기계 엡실론, 자릿수 손실과 그것을 피하는 식 변형, IEEE 754의 설계 이유를 다룬 고전.Oracle 문서
Floating-Point Arithmetic: Issues and Limitations튜토리얼Python 공식 문서. 0.1이 실제로는 3602879701896397 / 2⁵⁵로 저장된다는 것을 직접 확인하는 짧은 글.python.org
Fundamentals of Numerical Computation교재(무료)Driscoll·Braun. 근 찾기·선형 시스템·미분 방정식·편미분 방정식을 Julia·MATLAB·Python 판으로 나란히 보여 준다.fncbook.com
LAPACK라이브러리연립방정식·최소제곱·고윳값·특잇값을 푸는 Fortran 90 라이브러리. 계산은 최대한 BLAS에 맡긴다. NumPy·SciPy·Julia 밑에서 돈다.netlib
NumPy · SciPy라이브러리NumPy는 N차원 배열과 벡터화, SciPy는 그 위의 적분·최적화·보간·미분 방정식·희소 행렬 알고리즘.scipy.org
Julia언어LLVM으로 기계어로 컴파일되는 과학 계산용 언어. 다중 디스패치로 일반적인 코드를 빠르게 쓴다.julialang.org
Fortran언어과학·공학 계산을 위해 처음부터 설계된 언어. 배열 연산과 코어레이 병렬을 언어 차원에서 지원한다.fortran-lang.org