공식이 없으면 한 걸음씩: 수치해석이란

서윤의 첫 출근 날, 관제실 문이 벌컥 열렸습니다. 운항 팀장 민재가 태블릿을 흔들며 들어옵니다. "어젯밤 3호기가 비상 투하를 했는데, 상자가 땅에 닿을 때 얼마나 빨랐는지 아무도 몰라요. 지상 팀이 받을 수 있는 속도인지 알아야 하거든요."

서윤은 자신 있게 노트북을 엽니다. "물리 공식 있잖아요. 속도는 중력 가속도 곱하기 시간, v = gt." 계산해 보니 12초 뒤 속도가 117.72 m/s, 시속 420 km입니다. 민재가 고개를 젓습니다. "낙하산을 달았는데 그렇게 빠를 리가 없어요."

창가 자리의 오 박사가 커피잔을 내려놓습니다. "공기 저항을 뺐군. 저항을 넣으면 식이 갑자기 어려워지지. 공식으로 안 풀리면 한 걸음씩 가면 돼. 그게 수치해석이야."

이번 교시는 강의 전체의 예고편입니다. 하나의 문제를 공식으로 한 번, 한 걸음씩 한 번 풀어 보고, 두 답이 얼마나 다른지 재 봅니다. 앞으로 11교시 동안 배울 내용이 모두 이 세 단계(식 세우기, 근사해 풀기, 오차 재기)의 변주입니다.

1-1. 세상을 식으로 적기: 수학 모델 (8분)

상자에는 힘이 두 개 걸립니다. 아래로 당기는 중력 $mg$와, 위로 밀어 올리는 공기 저항입니다. 공기 저항은 빨리 떨어질수록 커지는데, 낙하산처럼 큰 물체는 속도의 제곱에 비례한다고 보면 잘 맞습니다. 비례 상수를 항력 계수(drag coefficient) $c_d$라고 부르겠습니다.

뉴턴의 운동 법칙(질량 × 가속도 = 힘의 합)을 적으면 이렇습니다.

$\dfrac{dv}{dt}$는 "속도가 1초에 얼마나 변하는가", 즉 속도 그래프의 기울기입니다. 식을 말로 읽으면 이렇습니다. 속도가 느릴 때는 거의 중력 가속도 $g$만큼 빨라지고, 빨라질수록 저항이 커져 덜 빨라진다.

이렇게 현실을 식으로 옮긴 것을 수학 모델(mathematical model) 이라고 합니다. 모델에는 언제나 가정이 숨어 있습니다.

$v$가 커져서 $g = \frac{c_d}{m}v^2$이 되면 기울기가 0, 즉 더는 빨라지지 않습니다. 이 속도를 종단 속도(terminal velocity) 라고 합니다.

시속 약 56 km입니다. 서윤의 시속 420 km보다는 훨씬 현실적이지만, 상자를 받기에는 여전히 빠릅니다. 이 문제는 4교시에 다시 풉니다.

1-2. 공식으로 풀기: 해석해 (8분)

이 식은 운 좋게도 미적분으로 풀립니다. 풀이 과정은 건너뛰고 결과만 보면 이렇습니다.

$\tanh$(하이퍼볼릭 탄젠트)는 0에서 출발해 1에 점점 다가가는 S자 함수입니다. 그래서 속도도 0에서 출발해 종단 속도에 점점 다가갑니다. 이렇게 식을 정리해 정확한 답을 공식으로 얻은 것을 해석해(analytical solution) 라고 합니다.

NumPy로 0초부터 12초까지 2초 간격으로 계산해 봅시다. np.linspace(0, 12, 7)은 0과 12 사이를 고르게 나눈 숫자 7개(양 끝 포함)를 만듭니다. 배열에 np.sqrt, np.tanh를 쓰면 원소마다 한꺼번에 계산됩니다.

>>> import numpy as np
>>> g, m, cd = 9.81, 12, 0.48
>>> t = np.linspace(0, 12, 7)
>>> t
array([ 0.,  2.,  4.,  6.,  8., 10., 12.])
>>> v_exact = np.sqrt(g * m / cd) * np.tanh(np.sqrt(g * cd / m) * t)
>>> np.round(v_exact, 2)
array([ 0.  , 13.3 , 15.45, 15.64, 15.66, 15.66, 15.66])

2초 만에 13.3 m/s, 4초면 거의 종단 속도입니다. np.round(배열, 2)는 보기 좋게 소수 둘째 자리까지 반올림한 것이고, 계산 자체는 더 많은 자릿수로 이루어집니다(2교시).

그런데 해석해는 아주 운 좋은 경우에만 있습니다. 저항이 속도의 1.5제곱에 비례한다면, 바람이 시간에 따라 바뀐다면, 낙하산이 2초 뒤에야 펴진다면 이런 깔끔한 공식은 없습니다. 현실의 문제는 대부분 그쪽입니다.

1-3. 한 걸음씩 풀기: 오일러 방법 (12분)

공식이 없을 때도 쓸 수 있는 방법을 만들어 봅시다. 출발점은 "기울기는 1초에 변하는 양"이라는 뜻 그대로입니다. 지금 속도 $v_i$를 알면, 지금 기울기 $g - \frac{c_d}{m}v_i^2$도 계산할 수 있습니다. 그 기울기가 $h$초 동안 그대로 유지된다고 치고, 한 걸음 나아갑니다.

새 속도 = 지금 속도 + 지금 기울기 × 걸음 크기. 이것이 오일러 방법(Euler's method) 입니다. 18세기 수학자 오일러의 이름이 붙었지만, 하는 일은 "지금 빠르기로 잠깐 가 보고, 도착한 곳에서 빠르기를 다시 재기"의 반복일 뿐입니다.

$h = 1$초로 손으로 두 걸음만 가 봅시다. $c_d/m = 0.48/12 = 0.04$입니다.

공식으로 구한 값은 $v(1) = 8.70$, $v(2) = 13.30$입니다. 오일러 방법이 조금 앞서 나갑니다. 걸음을 내딛는 동안 실제로는 기울기가 점점 줄어드는데, 처음 기울기를 그대로 썼기 때문입니다.

같은 계산을 반복문으로 적으면 이렇습니다. np.zeros_like(t)는 t와 같은 모양의 0 배열을 만들어 결과를 담을 자리를 마련합니다.

>>> h = 1
>>> t = np.linspace(0, 12, 13)
>>> v = np.zeros_like(t)
>>> for i in range(len(t) - 1):
...     v[i + 1] = v[i] + (g - cd / m * v[i]**2) * h
>>> np.round(v[:5], 2)
array([ 0.  ,  9.81, 15.77, 15.63, 15.67])
>>> v_exact = np.sqrt(g * m / cd) * np.tanh(np.sqrt(g * cd / m) * t)
>>> np.round(v_exact[:5], 2)
array([ 0.  ,  8.7 , 13.3 , 14.95, 15.45])

반복문 안에서 v[i]로 지금 값을 읽고 v[i + 1]에 다음 값을 씁니다. 배열 전체를 한 번에 계산하는 NumPy다운 방식이 아니라 한 칸씩 채우는 이유는, 다음 값이 이전 값에 기대기 때문입니다. 이런 계산을 점화식(recurrence) 이라고 하고, 수치해석에는 아주 많이 나옵니다.

아래 실습대에서 걸음 크기 h를 바꿔 보세요. 회색 곡선이 공식(해석해), 남색 점과 선이 오일러 방법입니다. 질량과 항력 계수 슬라이더를 움직이면 종단 속도(점선)가 바뀌고, '상자 떨어뜨리기'를 누르면 12초 동안의 낙하가 4초로 재생되며 두 계산이 나란히 그려집니다. 왼쪽 기둥은 그동안 떨어진 거리입니다.

1-4. 얼마나 틀렸나: 오차 재기 (8분)

근사해를 얻었으면 반드시 물어야 할 질문이 있습니다. 얼마나 틀렸나? 수치해석에서 답만 내고 오차를 말하지 않는 것은 측정값에 단위를 빼먹는 것과 같습니다.

$t = 2$초에서 오일러 방법($h = 1$)은 15.77, 참값은 13.30이니 $E_t = -2.47$ m/s, $\varepsilon_t \approx -18.6\%$입니다. 꽤 큽니다.

걸음을 줄이면 어떻게 될까요? 걸음 크기를 반씩 줄이며 12초 동안의 가장 큰 오차를 재 봅니다. 오일러 방법을 함수로 묶어 두면 여러 번 부르기 편합니다.

>>> def euler_fall(h, t_end=12):
...     t = np.linspace(0, t_end, round(t_end / h) + 1)
...     v = np.zeros_like(t)
...     for i in range(len(t) - 1):
...         v[i + 1] = v[i] + (g - cd / m * v[i]**2) * h
...     return t, v
>>> for h in [1, 0.5, 0.25, 0.125]:
...     t, v = euler_fall(h)
...     exact = np.sqrt(g * m / cd) * np.tanh(np.sqrt(g * cd / m) * t)
...     print(h, round(np.abs(v - exact).max(), 3))
1 2.474
0.5 0.982
0.25 0.477
0.125 0.232

0.5부터는 걸음을 반으로 줄일 때마다 오차도 거의 정확히 반이 됩니다. 이런 방법을 1차 정확도(first-order accurate) 라고 하고 오차가 $O(h)$라고 적습니다. "오차가 $h$에 비례한다"는 뜻입니다. 반대로 걸음 수는 두 배가 되니, 정확도 한 자리(오차 1/10)를 더 얻으려면 계산을 10배 더 해야 합니다. 12교시에서 걸음 수를 조금만 늘려도 오차가 확 줄어드는 방법(RK4, 오차 $O(h^4)$)을 배웁니다.

참값을 모르는 진짜 문제에서는 어떻게 할까요? 걸음 크기 $h$와 $h/2$로 두 번 계산해서 둘의 차이를 오차의 어림값으로 씁니다. 두 답이 거의 같아지면 충분히 작은 걸음이라고 판단합니다. 4교시의 정지 조건이 바로 이 생각입니다.

1-5. 수치해석이 하는 일 (6분)

오 박사가 화이트보드에 세 줄을 적습니다.

1. 식 세우기: 현실을 수학 모델로 옮긴다. 가정을 적어 둔다.
2. 근사해 풀기: 공식이 없으면 덧셈·곱셈의 반복으로 답에 다가간다.
3. 오차 재기: 그 답이 얼마나 믿을 만한지 숫자로 말한다.

컴퓨터는 덧셈, 뺄셈, 곱셈, 나눗셈과 비교만 빠르게 할 줄 압니다. 미분도, 적분도, 방정식 풀기도 결국 이 네 가지 연산의 반복으로 바꿔야 컴퓨터가 풀 수 있습니다. 그렇게 바꾸는 방법과, 바꾸면서 생기는 오차를 다루는 학문이 수치해석(numerical analysis) 입니다.

오차는 크게 네 군데서 생깁니다.

이 강의의 지도입니다. 하늘상자 관제실에서 하루하루 생기는 문제를 따라가다 보면 수치해석 입문서의 큰 갈래를 모두 지나게 됩니다.

책의 기호 읽는 법

자주 하는 실수 (3분)

참값이 없는데 오차를 말하기. 실제 문제에는 공식이 없으니 참 오차도 계산할 수 없습니다. 이때는 걸음 크기를 바꿔 두 번 계산하고 차이를 봅니다. "한 번 돌려서 나온 숫자"를 그대로 믿지 마세요.

배열 대신 숫자에 덮어쓰기. 반복문에서 v = v + (g - cd / m * v**2) * h처럼 쓰면 마지막 값 하나만 남습니다. 시간마다 값을 남기려면 v[i + 1] = ...처럼 배열의 칸을 채우세요.

np.linspace의 개수 착각. 0초부터 12초까지 1초 간격이면 칸은 12개지만 점은 13개입니다. np.linspace(0, 12, 12)라고 쓰면 간격이 약 1.09초가 되어 버립니다.

>>> np.linspace(0, 12, 12)[:3]
array([0.        , 1.09090909, 2.18181818])
>>> np.linspace(0, 12, 13)[:3]
array([0., 1., 2.])

연습 문제 (5분)

1. 질량이 두 배(24 kg)가 되면 종단 속도는 몇 배가 될까요? 먼저 식으로 예상한 뒤 NumPy로 확인해 보세요.
2. euler_fall로 h = 0.25일 때 t = 4초의 속도를 구하고, 공식 값과의 상대 참 오차(%)를 계산해 보세요.
3. 실습대에서 h를 2로 두면 출렁이는데, 항력 계수를 0.1로 줄이면 출렁임이 사라집니다. 왜 그럴지 종단 속도와 기울기의 크기로 생각해 보세요.

풀이

>>> np.sqrt(g * 24 / cd) / np.sqrt(g * 12 / cd)
np.float64(1.4142135623730951)
>>> t, v = euler_fall(0.25)
>>> t[16], round(v[16], 4)
(np.float64(4.0), np.float64(15.5577))
>>> exact4 = np.sqrt(g * m / cd) * np.tanh(np.sqrt(g * cd / m) * 4)
>>> round(exact4, 4)
np.float64(15.4532)
>>> round((exact4 - v[16]) / exact4 * 100, 2)
np.float64(-0.68)

1번은 종단 속도가 $\sqrt{m}$에 비례하므로 $\sqrt{2} \approx 1.41$배입니다. 무거운 상자일수록 큰 낙하산이 필요합니다.

2번은 h = 0.25초로 16걸음 가면 t = 4초입니다(t[16]). 오일러 방법이 약 0.68% 빠르게 나옵니다. 같은 시점에서 h = 1이면 15.67 m/s로 약 1.39% 차이였으니, 걸음을 1/4로 줄여 오차가 절반쯤으로 줄었습니다. 1-4의 표처럼 딱 1/4이 되지 않는 것은 한 시점의 오차가 그 앞 걸음들에서 쌓인 오차와 서로 상쇄되기도 하기 때문입니다. 방법끼리 견줄 때는 한 시점보다 전체 구간의 가장 큰 오차를 보는 편이 안전합니다.

3번은 항력 계수가 작으면 종단 속도 근처에서 기울기가 바뀌는 빠르기도 느려집니다. 오일러 방법은 "기울기가 h초 동안 그대로"라고 치는데, 기울기가 천천히 바뀌면 이 가정이 크게 틀리지 않습니다. 반대로 기울기가 빠르게 바뀌는 문제(항력이 큰 문제)일수록 작은 걸음이 필요합니다. 12교시에서 이 조건을 식으로 정리합니다.

퇴근길, 서윤은 정산 팀에서 온 메일을 봅니다. "어제 배송 요금 합계를 두 프로그램으로 냈더니 끝자리가 서로 다릅니다." 다음 교시에는 컴퓨터가 0.1이라는 숫자조차 정확히 저장하지 못한다는 사실에서 출발합니다.