계산 과학
풀 수 없는 방정식을 컴퓨터로 근사해 푸는 일 — 날씨 예보, 다리의 응력, 별의 궤도 — 의 바탕에는 유한한 비트로 수를 적는 방식과, 근사가 얼마나 틀리고 언제 폭발하는지에 대한 이해가 있다. 부동소수점에서 시작해 근 찾기·선형 시스템·미분 방정식·몬테카를로·열 방정식을 작은 예로 직접 계산하며 따라간다.
목차 (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 근처를 위한 절대 오차를 함께 쓴다. 값이 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자리를 잃는다"고 본다. 조건수가 큰 문제는 어떤 좋은 알고리즘으로도 정확히 풀 수 없다. 할 수 있는 것은 문제를 다시 세우는 것이다.
덧셈 순서도 결과를 바꾼다. 부동소수점 덧셈은 결합 법칙이 성립하지 않아 (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걸음이면 된다. 미분이 필요 없고 절대 발산하지 않는다는 것이 장점이다.
1단계: [1, 2]의 가운데 1.5에서 f = 0.2500. 부호가 +이니 근은 왼쪽 반에 있다.
대신 느리다. 한 걸음에 이진 자리 하나씩만 얻는다(1차 수렴). 접선을 따라 근으로 뛰는 뉴턴 방법은 근 가까이에서 맞는 자릿수가 걸음마다 두 배로 늘어(2차 수렴) 훨씬 빠르지만, 도함수가 필요하고 출발점이 나쁘면 엉뚱한 곳으로 튄다. 뉴턴 방법의 원리와 예는 미적분 노트의 미분의 쓰임 절에 있다. 실무 라이브러리(SciPy의 brentq 등)는 이분법으로 근을 늘 감싸 두면서 가능할 때 빠른 보간 걸음을 섞는 방식을 쓴다.
선형 시스템: 직접법과 반복법
시뮬레이션의 대부분은 결국 Ax = b를 푸는 일로 바뀐다. 미분 방정식을 격자 위에서 이산화하면 격자점 수만큼의 미지수를 가진 연립방정식이 나온다. 크기가 작거나 중간이면 직접법을 쓴다. 가우스 소거로 A를 삼각 행렬 둘의 곱(LU 분해)으로 바꾸고 앞뒤로 대입한다. 계산량은 n³에 비례한다. 가우스 소거의 원리는 선형대수 노트의 연립방정식 절에 있다. 수치 계산에서는 0에 가까운 수로 나누지 않도록 행을 바꾸는 피벗팅이 꼭 필요하다.
격자가 3차원으로 커지면 미지수가 수백만 개를 넘고 n³은 감당할 수 없다. 다행히 이런 행렬은 거의 모든 칸이 0인 희소 행렬이다(각 격자점은 이웃 몇 개와만 엮인다). 이때는 답을 추측하고 고쳐 나가는 반복법이 낫다. 가장 단순한 야코비 반복은 i번째 식을 xᵢ에 대해 풀고, 오른쪽에는 직전 반복의 값을 넣는다. 대각 원소가 같은 행의 나머지보다 크면(대각 우세) 수렴이 보장된다. 그림의 3×3 예에서 10번 반복하면 잔차가 1.8e-4까지 준다. 실제로는 방금 고친 값을 바로 쓰는 가우스–자이델, 대칭 양의 정부호 행렬에 강한 켤레 기울기법(CG), 격자를 굵게·가늘게 오가는 멀티그리드를 쓴다.
미분 방정식 시뮬레이션: 오일러와 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로 오래 돌리면 에너지가 조금씩 새어 궤도가 서서히 줄거나 커진다. 위치와 속도를 엇갈려 갱신하는 심플렉틱 방법(리프프로그·벨레)은 차수가 낮아도 에너지를 오래 지켜, 분자 동역학과 천체 계산의 표준이다.
안정성과 시간 간격: 강성 문제
정확도와 다른 문제가 안정성이다. 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보다 작다. 회로 시뮬레이터나 화학 반응 계산은 이런 강성용 풀이기를 쓴다.
몬테카를로: 무작위로 적분하기
몬테카를로 방법은 답을 무작위 표본의 평균으로 추정한다. 한 변이 1인 정사각형에 점을 고르게 던지면 사분원 안에 들어갈 확률이 π/4이니, 안에 든 비율에 4를 곱하면 π가 나온다. 이 노트의 고정된 난수열로 점 100개를 던지면 3.2, 100,000개면 3.13732가 나왔다.
핵심 성질은 오차가 1/√N에 비례한다는 것이다. 표본 평균의 표준 편차가 그렇게 줄기 때문이다. 오차를 1/10로 줄이려면 표본이 100배 필요하니 느리다. 대신 이 속도가 차원과 무관하다. 격자로 적분하면 차원마다 칸 수가 곱해져 10차원만 돼도 감당할 수 없지만, 몬테카를로는 몇 차원이든 같은 비율로 준다. 그래서 입자 물리의 검출기 시뮬레이션, 금융의 옵션 가격, 수많은 빛의 경로를 평균하는 렌더링(컴퓨터 그래픽스 노트의 레이 트레이싱 절)에 쓰인다. 표본을 고르게 퍼뜨리는 준난수, 중요한 곳에 표본을 몰아 주는 중요도 표집이 같은 N으로 오차를 줄이는 기법이다.
편미분 방정식: 열 방정식과 유한 차분
공간과 시간에 함께 걸친 현상(열의 확산, 파동, 유체)은 편미분 방정식(PDE)으로 적힌다. 가장 단순한 것이 1차원 열 방정식 ut = α uxx다. "온도가 이웃보다 높으면 내려가고 낮으면 올라간다"는 뜻이다. 유한 차분법은 막대를 Δx 간격의 점으로, 시간을 Δt 걸음으로 자르고 미분을 이웃 값의 차로 바꾼다. 시간은 앞으로 한 걸음, 공간은 가운데 차분을 쓰는 FTCS 방식이면 매 걸음 ui ← ui + r(ui−1 − 2ui + ui+1), r = αΔt/Δx²다.
막대를 칸 10개(격자점 11개)로 나누고 가운데 세 점만 온도 1로 시작한다. 양 끝은 늘 0이다.
앞 절의 안정성 문제가 여기서도 나타난다. 이 방식은 r ≤ 1/2일 때만 안정하다. 그림의 r = 0.4는 매끄럽게 퍼지지만, r = 0.6는 이웃끼리 톱니처럼 엇갈리며 커진다. 공간을 두 배로 잘게 나누면 Δt는 네 배로 작아져야 하니, 세밀한 격자에서는 걸음 수가 빠르게 늘어난다. 그래서 크랭크–니콜슨 같은 암시적 방식으로 매 걸음 (희소한) 연립방정식을 풀어 큰 Δt를 쓰거나, 문제에 맞춰 둘을 고른다. 파동 방정식에서는 같은 역할을 CFL 조건(한 걸음에 정보가 한 칸보다 멀리 가면 안 됨)이 한다.
유한 요소법
유한 차분은 반듯한 격자에서 편하지만, 자동차 차체나 뼈처럼 모양이 복잡하면 격자를 맞추기 어렵다. 유한 요소법(FEM)은 영역을 삼각형·사면체 같은 작은 요소로 나누고(메시), 해를 요소마다 단순한 함수(대개 1차식)를 이어 붙인 것으로 근사한다. 1차원으로 보면 마디마다 그 마디에서 1, 이웃 마디에서 0인 "모자 함수"를 두고, 해를 이 함수들의 가중합으로 적는다. 가중치를 정하는 방법은 방정식을 그대로 맞추는 대신 "적분했을 때 평균적으로 맞게"(약형식) 요구하는 것이다. 결과는 큰 희소 연립방정식이 되어 앞 절의 반복법으로 푼다.
장점은 모양에 맞춰 요소를 자유롭게 배치하고, 응력이 몰리는 모서리·구멍 둘레만 잘게 나눌 수 있다는 것이다. 구조 해석(다리·건물·비행기 날개), 열 전달, 전자기장 해석의 표준 도구다. 메시의 질이 결과를 좌우해서, 지나치게 찌그러진 요소는 정확도를 크게 떨어뜨린다. 메시를 점점 줄여 가며 답이 더 변하지 않는지 확인하는 수렴 시험이 기본 절차다.
N체 시뮬레이션과 날씨·기후 모델
N체 문제는 별·은하·분자처럼 서로 끌어당기는 입자 N개의 움직임을 계산한다. 모든 쌍의 힘을 더하면 N²에 비례하는 일이 필요해, 입자가 수백만 개면 감당할 수 없다. 반스–헛(Barnes–Hut) 같은 트리 방법은 공간을 나무처럼 나누고 멀리 있는 무리를 질량 중심 하나로 묶어 대략 N log N으로 줄인다. 시간 적분에는 앞에서 본 심플렉틱 방법을 써서 에너지를 오래 지킨다.
날씨·기후 모델은 대기(와 바다)를 가로·세로 격자로 나누고, 칸마다 온도·기압·바람·습도를 두어 유체 역학과 열역학 방정식을 시간에 따라 푼다. 구름이나 난류처럼 칸보다 작은 현상은 직접 풀 수 없어 근사식(매개변수화)으로 넣는데, 이것이 모델 간 차이의 큰 원인이다. 대기는 초기값의 작은 차이가 커지는 혼돈계라, 날씨 예보는 초기값을 조금씩 바꾼 여러 번의 실행(앙상블)으로 불확실성까지 함께 내놓는다. 기후 모델은 개별 날씨가 아니라 수십 년 평균의 통계를 묻는다는 점이 다르다. 두 경우 모두 격자를 여러 프로세서에 나눠 맡기고 경계 값을 주고받는 병렬 계산이 필수이며, 그 방법은 병렬 프로그래밍 노트에 있다.
검증과 확인(V&V)
시뮬레이션이 그럴듯한 그림을 내놓는다고 맞는 것은 아니다. 계산 과학은 두 질문을 구별한다. 검증(verification)은 "코드가 방정식을 맞게 푸는가"다. 정확한 해가 알려진 문제로 시험하고, 격자나 걸음을 줄일 때 오차가 이론 차수대로 줄어드는지 본다(앞의 오일러 1차, RK4 4차가 바로 그런 시험이다). 정확한 해가 없으면 답을 먼저 정하고 그에 맞는 원천 항을 거꾸로 만들어 넣는 "만든 해" 방법을 쓴다. 확인(validation)은 "맞는 방정식을 푸는가"다. 결과를 실험·관측과 비교한다. 코드가 완벽해도 모델의 가정(마찰 없음, 2차원 근사)이 틀리면 답은 틀린다.
여기에 불확실성 정량화(UQ)가 더해진다. 입력 매개변수(재료 강도, 초기 조건)에 오차 범위가 있으면 답도 하나의 수가 아니라 범위여야 한다. 입력을 표본으로 바꿔 여러 번 돌리는 몬테카를로가 가장 단순한 방법이다. 재현성도 같은 맥락이다. 난수 씨앗, 라이브러리 버전, 컴파일러 최적화(부동소수점 덧셈 순서가 바뀌면 마지막 자리가 달라진다)를 기록해 두어야 같은 결과를 다시 낼 수 있다.
도구: 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)가 일을 나눈다. 직접 수치 알고리즘을 짜기 전에 검증된 라이브러리 함수가 있는지 찾는 것이 첫 번째 원칙이다.
핵심 자료
여러 코어·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 |
- What Every Computer Scientist Should Know About Floating-Point Arithmetic — Oracle 문서
- Floating-Point Arithmetic: Issues and Limitations — python.org
- Fundamentals of Numerical Computation — fncbook.com
- LAPACK — netlib
- NumPy · SciPy — scipy.org
- Julia — julialang.org
- Fortran — fortran-lang.org
- NumPy
- math.isclose — Python 문서
- numpy.sum (부분 쌍별 합산) — NumPy 문서
- IEEE 754 — Wikipedia