처음 문제로 돌아가기: 미분방정식의 수치해

서윤이 입사한 지 꼭 2주가 되는 아침입니다. 관제실 창밖으로 바다 쪽 깃발이 팽팽하게 펄럭입니다. 민재 팀장이 들어오며 말합니다. "오늘 아침 투하 연습, 그대로 간대요. 4교시에 정한 큰 낙하산으로요. 그런데 바람이 초속 5 m예요. 상자가 옆으로 밀려갈 텐데, 지상 팀이 어디 서 있어야 해요?"

서윤이 웃으며 노트북을 엽니다. "첫날 했던 오일러 방법 그대로 쓰면 되죠. h를 아주 작게 잡고요." 오 박사가 고개를 끄덕입니다. "되긴 하지. 그런데 그 방법은 오차를 1/10로 줄이려면 계산을 열 배 해야 했어. 오늘은 같은 계산으로 훨씬 정확해지는 방법을 쓰자. 그리고 이번엔 속도 하나가 아니야. 가로 위치, 세로 위치, 가로 속도, 세로 속도를 함께 풀어야 해."

그때 정비 팀의 메시지가 옵니다. "오늘 낮 35도래요. 비행 중 배터리 온도가 얼마까지 오를지도 계산해 줄 수 있어요?"

이번 교시는 강의의 마지막입니다. 1교시의 질문 "공식이 없으면 어떻게 푸나"로 돌아가, 그동안 배운 도구(테일러 급수, 오차 차수, 근 찾기, 보간, 적분)를 모두 써서 미분방정식을 제대로 풉니다.

12-1. 다시 처음 문제로: 초깃값 문제와 오일러 (7분)

1교시의 식을 다시 봅시다. 속도가 변하는 빠르기(기울기)가 지금 속도에 따라 정해집니다.

이렇게 도함수가 들어 있는 식을 미분방정식(differential equation)이라고 하고, 변수가 시간 하나뿐이면 상미분방정식(ODE, ordinary differential equation) 이라고 합니다. 일반적인 모양은 이렇습니다.

$f(t, y)$는 "시각 $t$에 값이 $y$이면 기울기가 얼마인가"를 알려 주는 함수이고, $y(t_0) = y_0$은 출발점입니다. 출발점이 주어진 문제를 초깃값 문제(IVP, initial value problem) 라고 합니다. 수치해법은 모두 같은 일을 합니다. 지금 점 $(t_i, y_i)$에서 기울기를 재고, 한 걸음 $h$만큼 나아가 다음 점 $y_{i+1}$을 정합니다. 방법마다 다른 것은 기울기를 어떻게 재느냐 하나뿐입니다.

그래서 코드도 "한 걸음 함수"와 "반복하는 함수"로 나눠 두면 편합니다. solve는 걸음 함수를 인자로 받으므로, 뒤에서 방법을 바꿀 때 걸음 함수만 갈아 끼우면 됩니다. np.shape(y0)를 쓴 것은 $y$가 숫자일 때와 벡터일 때(12-4) 모두 쓰기 위해서입니다.

>>> g, m, cd = 9.81, 12, 0.48
>>> def f(t, v):
...     return g - cd / m * v**2
>>> def exact(t):
...     return np.sqrt(g * m / cd) * np.tanh(np.sqrt(g * cd / m) * t)
>>> def euler_step(f, t, y, h):
...     return y + h * f(t, y)
>>> def solve(step, f, y0, t_end, h):
...     t = np.linspace(0, t_end, round(t_end / h) + 1)
...     y = np.zeros((len(t),) + np.shape(y0))
...     y[0] = y0
...     for i in range(len(t) - 1):
...         y[i + 1] = step(f, t[i], y[i], h)
...     return t, y
>>> t, v = solve(euler_step, f, 0.0, 12, 1)
>>> round(np.abs(v - exact(t)).max(), 3)
np.float64(2.474)

1교시의 2.474가 그대로 나옵니다. 그 표에서 오차는 $h$에 비례했습니다($O(h)$). 그런데 한 걸음만 떼어 보면 이야기가 다릅니다. 정확한 값에서 출발해 딱 한 걸음 간 뒤의 오차를 국소 오차(local error) 라고 합니다.

>>> for h in [0.4, 0.2, 0.1]:
...     one = euler_step(f, 2.0, exact(2.0), h)
...     print(h, f"{abs(one - exact(2.0 + h)):.2e}")
0.4 2.07e-01
0.2 5.50e-02
0.1 1.41e-02

$h$를 반으로 하면 한 걸음의 오차는 약 1/4입니다. 3교시의 테일러 급수로 보면 당연합니다. 오일러는 $y(t+h) = y + hy' + \frac{h^2}{2}y'' + \cdots$에서 $h^2$ 항부터 버렸으니 한 걸음 오차가 $O(h^2)$입니다. 그런데 끝 시각까지 가려면 걸음을 $T/h$번 밟아야 하고, 한 걸음 오차가 쌓이면 $\frac{T}{h} \times h^2 = T h$, 곧 전역 오차(global error) 는 $O(h)$가 됩니다. 방법의 "차수"는 이 전역 오차의 차수로 말합니다.

12-2. 기울기를 두 번 재기: 호인과 중점 방법 (8분)

오일러는 걸음의 출발점 기울기를 걸음 내내 씁니다. 1교시에서 봤듯이 속도가 오를수록 기울기는 줄어드는데, 처음 기울기만 쓰니 늘 앞서 나갑니다. 11교시의 왼쪽 직사각형과 같은 실수입니다. 그렇다면 사다리꼴처럼 양 끝 기울기의 평균을 쓰면 됩니다. 문제는 끝점의 값을 아직 모른다는 것입니다. 그래서 두 단계로 나눕니다.

1. 예측(predictor): 오일러로 끝점을 어림한다. $\tilde{y} = y_i + h\,f(t_i, y_i)$
2. 수정(corrector): 출발점 기울기와 어림한 끝점의 기울기를 평균 내서 다시 간다.

이것이 호인 방법(Heun's method) 입니다. 기울기를 걸음 한가운데서 한 번 재는 방법도 있습니다. 반 걸음 가 보고 그곳의 기울기로 한 걸음 전체를 갑니다. 11교시의 중점 합과 같은 생각이라 중점 방법(midpoint method) 이라고 부릅니다.

$h = 1$로 첫 걸음을 손으로 가 봅니다. $c_d/m = 0.04$입니다.

정확한 값 8.7007에 오일러(9.81)는 1.11, 호인은 0.82, 중점은 0.15만큼 벗어났습니다. 기울기를 두 번 재는 대신 오차가 크게 줄었습니다.

>>> def heun_step(f, t, y, h):
...     k1 = f(t, y)
...     k2 = f(t + h, y + h * k1)
...     return y + h / 2 * (k1 + k2)
>>> def midpoint_step(f, t, y, h):
...     k1 = f(t, y)
...     k2 = f(t + h / 2, y + h / 2 * k1)
...     return y + h * k2
>>> round(heun_step(f, 0, 0.0, 1), 4), round(midpoint_step(f, 0, 0.0, 1), 4), round(exact(1), 4)
(7.8853, 8.8476, np.float64(8.7007))

두 방법 모두 테일러 급수의 $h^2$ 항까지 맞혀서 한 걸음 오차가 $O(h^3)$, 전역 오차가 $O(h^2)$입니다. 2차 방법입니다. 걸음을 반으로 줄이면 오차가 1/4이 됩니다.

12-3. 기울기를 네 번 재는 RK4 (10분)

기울기를 더 여러 번, 더 영리한 곳에서 재면 차수를 더 올릴 수 있습니다. 이런 방법을 통틀어 룽게-쿠타 방법(Runge–Kutta methods) 이라고 하고, 가장 널리 쓰이는 것이 기울기를 네 번 재는 고전적 4차 룽게-쿠타(RK4) 입니다.

가중치 1-2-2-1이 낯익지 않나요? 11교시 심프슨 공식의 1-4-1에서 가운데 4를 두 번의 중간 기울기($k_2$, $k_3$)로 나눠 가진 것입니다. 실제로 $f$가 $t$에만 기대는 경우 RK4는 심프슨 공식과 똑같아집니다. 그래서 오차도 심프슨처럼 $O(h^4)$입니다.

$h = 1$로 첫 걸음을 손으로 가 봅니다.

정확한 값 8.7007과 0.016 차이입니다. 오일러의 1.11보다 70배 넘게 정확합니다.

>>> def rk4_step(f, t, y, h):
...     k1 = f(t, y)
...     k2 = f(t + h / 2, y + h / 2 * k1)
...     k3 = f(t + h / 2, y + h / 2 * k2)
...     k4 = f(t + h, y + h * k3)
...     return y + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
>>> round(rk4_step(f, 0, 0.0, 1), 4)
8.685

이제 네 방법을 같은 문제(0~12초)에서 걸음을 반씩 줄이며 견줍니다. 오차 열은 12초 동안 가장 큰 오차이고, 마지막 열은 11교시처럼 로그-로그 기울기입니다.

>>> steps = {"euler": euler_step, "heun": heun_step, "mid": midpoint_step, "rk4": rk4_step}
>>> hs = np.array([1, 0.5, 0.25, 0.125])
>>> for name, step in steps.items():
...     errs = []
...     for h in hs:
...         t, v = solve(step, f, 0.0, 12, h)
...         errs.append(np.abs(v - exact(t)).max())
...     slope = np.polyfit(np.log(hs), np.log(errs), 1)[0]
...     print(f"{name:6s}", " ".join(f"{e:.2e}" for e in errs), f"{slope:.2f}")
euler  2.47e+00 9.82e-01 4.77e-01 2.32e-01 1.13
heun   1.47e+00 2.55e-01 5.44e-02 1.25e-02 2.29
mid    7.00e-01 1.26e-01 2.64e-02 6.07e-03 2.28
rk4    6.47e-02 3.01e-03 1.61e-04 9.30e-06 4.25

기울기가 1, 2, 4 언저리로 나옵니다(h = 1은 아직 걸음이 커서 기울기가 조금 큽니다). RK4는 걸음을 반으로 줄일 때마다 오차가 약 1/16, 걸음 8개만 더 쓰면 한 자리 이상 정확해집니다.

공정하게 견주려면 계산량도 봐야 합니다. 걸음 하나에 오일러는 $f$를 1번, 호인·중점은 2번, RK4는 4번 부릅니다. 12초를 $f$ 48번으로 푼다면 오일러는 $h = 0.25$(오차 0.477), 호인은 $h = 0.5$(0.255), RK4는 $h = 1$(0.0647)입니다. 같은 비용으로 RK4가 일곱 배 넘게 정확하고, 정확도를 더 요구할수록 차이는 기하급수로 벌어집니다. 실무의 기본 선택이 RK4 계열인 이유입니다.

12-4. 바람 부는 날: 연립 미분방정식 (11분)

먼저 낙하산부터 정합니다. 4교시에 이분법으로 "4초 뒤 속도가 10 m/s"가 되는 항력 계수를 구했습니다. 같은 계산을 짧게 다시 합니다.

>>> fc = lambda c: np.sqrt(g * m / c) * np.tanh(np.sqrt(g * c / m) * 4) - 10
>>> lo, hi = 0.5, 3.0
>>> for _ in range(40):
...     c = (lo + hi) / 2
...     if fc(lo) * fc(c) <= 0:
...         hi = c
...     else:
...         lo = c
>>> round(c, 4), round(np.sqrt(g * m / c), 2)
(1.1754, np.float64(10.01))

큰 낙하산의 $c_d \approx 1.1754$ kg/m, 종단 속도는 약 10.01 m/s입니다. 4초면 거의 종단 속도에 닿습니다.

이제 바람입니다. 상자는 가로($x$)와 세로($y$, 위가 +)로 움직이고, 공기 저항은 공기에 대한 상대 속도의 반대 방향으로 작용합니다. 바람이 $x$ 방향으로 $w$ m/s로 불면, 상자가 공기에 대해 움직이는 속도는 $\mathbf{u} = (v_x - w,\ v_y)$입니다. 저항의 크기는 1교시처럼 속력의 제곱 $c_d|\mathbf{u}|^2$이고 방향은 $-\mathbf{u}$이므로, 힘은 $-c_d|\mathbf{u}|\,\mathbf{u}$입니다. 위치의 변화율은 속도, 속도의 변화율은 가속도이니 식 네 개가 한 묶음이 됩니다.

이것이 연립 상미분방정식(system of ODEs) 입니다. 겁먹을 필요는 없습니다. 상태 네 개를 벡터 하나 $\mathbf{s} = [x, y, v_x, v_y]$로 묶으면 $\frac{d\mathbf{s}}{dt} = \mathbf{f}(t, \mathbf{s})$, 앞의 꼴과 똑같습니다. 오른쪽 식을 계산해 NumPy 배열로 돌려주는 함수만 만들면 됩니다. 드론은 바람을 거슬러 제자리에 떠 있다가 상자를 놓으므로, 처음 상태는 높이 60 m, 속도 0입니다.

>>> cd2 = 1.1754
>>> def wind_rhs(t, s, w):
...     x, y, vx, vy = s
...     ux, uy = vx - w, vy
...     k = cd2 / m * np.hypot(ux, uy)
...     return np.array([vx, vy, -k * ux, -g - k * uy])
>>> s0 = np.array([0.0, 60.0, 0.0, 0.0])
>>> wind_rhs(0, s0, 5.0)
array([ 0.     ,  0.     ,  2.44875, -9.81   ])

놓는 순간 상자는 가만히 있지만, 공기가 5 m/s로 스쳐 가니 바람 방향으로 2.45 m/s²의 가속도를 받습니다. np.hypot(a, b)는 $\sqrt{a^2 + b^2}$입니다.

놀라운 점은 12-3에서 만든 rk4_step을 한 글자도 고치지 않고 쓸 수 있다는 것입니다. y + h / 2 * k1 같은 식이 숫자에도, NumPy 배열에도 똑같이 동작하기 때문입니다. 남은 일은 땅에 닿는 순간을 찾는 것입니다. 걸음 단위로 가다 보면 $y$가 0을 건너뛰어 음수가 됩니다. 마지막 두 걸음 사이를 직선으로 잇고 $y = 0$인 곳을 찾습니다. 10교시의 구간별 선형 보간, np.interp 그대로입니다(np.interp는 두 번째 인자가 커지는 순서여야 해서 음수인 새 값을 앞에 둡니다).

>>> def drop(step, h, w, y0=60.0):
...     f2 = lambda t, s: wind_rhs(t, s, w)
...     t, s = 0.0, np.array([0.0, y0, 0.0, 0.0])
...     while s[1] > 0:
...         prev, s = s, step(f2, t, s, h)
...         t += h
...     t_land = np.interp(0.0, [s[1], prev[1]], [t, t - h])
...     x_land = np.interp(0.0, [s[1], prev[1]], [s[0], prev[0]])
...     return round(float(t_land), 3), round(float(x_land), 3)
>>> drop(rk4_step, 0.5, 5.0)
(6.752, 26.789)
>>> for h in [1.0, 0.5, 0.25, 0.1]:
...     print(h, drop(euler_step, h, 5.0), drop(heun_step, h, 5.0), drop(rk4_step, h, 5.0))
1.0 (7.027, 27.565) (6.895, 27.459) (6.762, 26.844)
0.5 (6.856, 27.022) (6.76, 26.893) (6.752, 26.789)
0.25 (6.799, 26.874) (6.752, 26.809) (6.752, 26.787)
0.1 (6.769, 26.815) (6.752, 26.79) (6.752, 26.787)

상자는 약 6.75초 뒤, 드론 바로 아래에서 바람 방향으로 26.79 m 떨어진 곳에 닿습니다. RK4는 $h = 0.5$(함수 계산 약 56번)로 이미 2 mm 안에 들었고, 오일러는 $h = 0.1$(약 68번)로도 3 cm 가까이 벗어납니다. 참값을 모를 때는 1교시에서 말한 대로 걸음을 반으로 줄여 두 번 풀고 차이를 봅니다.

>>> drop(rk4_step, 0.1, 5.0), drop(rk4_step, 0.05, 5.0)
((6.752, 26.787), (6.752, 26.787))
>>> [drop(rk4_step, 0.1, w)[1] for w in [0, 2, 4, 6, 8]]
[0.0, 10.347, 21.176, 32.519, 44.3]

두 번 풀어 소수 셋째 자리까지 같으니 믿어도 됩니다. 바람이 세질수록 착지점은 바람 세기에 거의 비례해 멀어집니다. 큰 낙하산 덕분에 상자가 금세 바람을 따라잡아 거의 바람과 같은 빠르기로 흘러가기 때문입니다. 민재 팀장에게 보낼 답은 "바람 5 m/s면 바람 방향으로 27 m 앞, 8 m/s면 44 m 앞"입니다.

아래 실습대에는 두 가지 실습이 있습니다. '속도 비교 (1교시 식)'에서는 회색 곡선이 정확한 해, 남색 점과 선이 고른 방법의 수치해입니다. 방법(오일러·호인·RK4)과 걸음 크기 h를 바꾸면 오른쪽 표에 같은 h에서 세 방법의 함수 계산 횟수와 최대 오차가 나란히 나옵니다. '바람 속 낙하'로 바꾸면 옆에서 본 낙하 경로가 나옵니다. 회색 점선이 아주 작은 걸음(RK4, h = 0.01)으로 구한 기준 경로, 남색이 고른 방법의 경로, 주황 점이 착지점입니다. 바람 슬라이더로 0~8 m/s를 고르고 '상자 떨어뜨리기'를 누르면 상자가 계산된 경로를 따라 떨어집니다.

12-5. 안정성과 걸음 크기 (6분)

1교시에서 오일러를 $h = 2$로 돌리면 종단 속도를 넘어 출렁이다 튀어 나갔습니다. 왜 그럴까요? 종단 속도 $v_T$ 근처에서 차이 $e = v - v_T$만 보면, 식은 거의 이렇게 됩니다(기울기를 $v_T$에서 직선으로 어림, 3교시의 테일러 1차 근사).

차이가 지수적으로 줄어드는 얌전한 식입니다. 그런데 오일러로 한 걸음 가면 $e_{i+1} = e_i - h\lambda e_i = (1 - h\lambda)\,e_i$입니다. 걸음마다 차이에 $(1 - h\lambda)$가 곱해지니, 이 수의 절댓값이 1보다 크면 차이가 줄기는커녕 커집니다. 음수면 부호가 번갈아 바뀌며 출렁입니다.

이것이 오일러의 안정 조건(stability condition) 입니다.

>>> lam = 2 * cd / m * np.sqrt(g * m / cd)
>>> round(lam, 3), round(2 / lam, 3)
(np.float64(1.253), np.float64(1.596))
>>> for h in [0.5, 1.0, 1.5, 2.0]:
...     vs = [0.0]
...     for i in range(40):
...         vs.append(euler_step(f, 0, vs[-1], h))
...     print(h, round(1 - lam * h, 2), np.round(vs[-3:], 2))
0.5 0.37 [15.66 15.66 15.66]
1.0 -0.25 [15.66 15.66 15.66]
1.5 -0.88 [15.67 15.65 15.67]
2.0 -1.51 [ 6.44 22.74  0.99]

1교시 상자는 $h < 1.6$초여야 안정합니다. $h = 1.5$는 출렁이지만($1 - h\lambda < 0$) 점점 가라앉고, $h = 2$는 출렁임이 커지며 엉뚱한 값을 오갑니다. 큰 낙하산은 $\lambda = 2g/v_T \approx 1.96$이라 오일러가 $h < 1.02$초여야 합니다. 실습대에서 오일러 $h = 1$의 세로 속도가 종단 속도 위아래로 출렁인 이유입니다. $1 - h\lambda \approx -0.96$이라 출렁임이 걸음마다 거의 줄지 않습니다. 저항이 클수록(값이 빨리 자리 잡는 문제일수록) $\lambda$가 커지고 허용되는 걸음은 작아집니다.

불안정이 늘 요란하게 드러나지는 않습니다. 실습대에서 호인 방법을 $h = 2$로 돌리면 속도가 3.07 m/s에서 딱 멈춥니다. 출렁이지도 튀지도 않으니 그럴듯해 보이지만, 종단 속도 15.66 m/s와는 전혀 다른 가짜 평형입니다. 걸음이 너무 크면 수치해법이 원래 식에 없는 답을 만들어 내기도 합니다. 이때도 해법은 같습니다. $h$를 반으로 줄여 다시 풀면 답이 크게 바뀌니 바로 들킵니다.

RK4도 한계가 있지만 넉넉합니다($h\lambda < 2.79$ 정도). 그래도 $\lambda$가 아주 큰 문제, 예를 들어 몇 밀리초 만에 자리 잡는 전기 회로와 몇 시간에 걸쳐 변하는 온도가 한 식에 섞인 문제는 정확도가 아니라 안정성 때문에 걸음을 아주 작게 잡아야 합니다. 이런 문제를 경직된(stiff) 문제라고 하고, 다음 값을 방정식으로 풀어 구하는 음함수 방법(implicit method) 을 씁니다.

걸음 크기는 어떻게 정할까요? 실무 프로그램은 걸음마다 오차를 어림해서 걸음을 스스로 조절합니다. 가장 쉬운 방법은 이 교시에서 한 그대로입니다. 한 걸음을 $h$로 한 번, $h/2$로 두 번 가 보고 두 답의 차이를 오차로 봅니다. 차이가 허용치보다 크면 걸음을 줄이고, 아주 작으면 다음 걸음을 키웁니다. 곡선이 급하게 변하는 곳(상자를 놓은 직후)에서는 촘촘하게, 잔잔한 곳(종단 속도)에서는 성큼성큼 갑니다. 이것이 적응 걸음(adaptive step size) 이고, 4차와 5차 공식을 함께 계산해 차이를 보는 RK45가 대표입니다.

실무에서는

SciPy의 scipy.integrate.solve_ivp가 이 교시의 모든 일을 합니다. 이 사이트의 서버에는 SciPy가 없어서 이름과 쓰임만 적습니다.

from scipy.integrate import solve_ivp

def ground(t, s, w):
    return s[1]              # 높이 y. 이 값이 0이 되는 순간이 '사건'
ground.terminal = True       # 사건이 일어나면 계산을 멈춘다
ground.direction = -1        # 위에서 아래로 지날 때만

sol = solve_ivp(wind_rhs, (0, 30), [0.0, 60.0, 0.0, 0.0], args=(5.0,),
                method="RK45", events=ground, rtol=1e-8, atol=1e-8)
t_land = sol.t_events[0][0]      # 착지 시각
x_land = sol.y_events[0][0][0]   # 착지 x

method="RK45"가 기본값이고 적응 걸음을 씁니다. 걸음 크기 대신 허용 오차 rtol(상대), atol(절대)를 줍니다. events는 "y가 0이 되는 순간"처럼 사건을 정확히 찾아 주므로, 12-4에서 직접 한 선형 보간보다 정확합니다. 경직된 문제에는 method="Radau"나 "BDF"(음함수 방법)를 씁니다. 매끈한 곡선이 필요하면 dense_output=True로 걸음 사이 값도 얻습니다.

책의 기호 읽는 법

자주 하는 실수 (3분)

벡터 상태를 리스트로 돌려주기. rhs가 파이썬 리스트를 돌려주면 y + h * k1이 엉뚱하게 됩니다. 리스트끼리 더하면 이어 붙이고, 리스트에 실수를 곱하면 오류가 납니다. 늘 np.array(...)로 돌려줍니다.

>>> [1.0, 2.0] + [3.0, 4.0]
[1.0, 2.0, 3.0, 4.0]
>>> np.array([1.0, 2.0]) + np.array([3.0, 4.0])
array([4., 6.])

기울기를 재는 도중에 상태 바꾸기. s[2] += h * k1[2]처럼 배열 칸을 그 자리에서 고치면, 뒤에서 k2를 계산할 때 이미 바뀐 값을 쓰게 됩니다. 한 걸음은 늘 새 배열(y + h * k)을 만들어 돌려줍니다.

땅을 지나친 마지막 점을 그대로 쓰기. 루프가 끝났을 때의 s는 이미 땅속입니다. 그 x를 착지점으로 쓰면 오차가 $h$에 비례합니다(오일러 수준). RK4로 공들여 풀고 마지막에 1차 오차를 붙이는 셈이니, 마지막 두 점을 보간하거나 events를 씁니다.

"그래프가 괜찮아 보이니 됐다". 걸음을 반으로 줄여 다시 풀고 답이 몇 자리까지 같은지 확인하기 전에는 답의 정확도를 말할 수 없습니다.

연습 문제 (5분)

1. 배터리 온도는 바깥 공기로 식고 모터 전류로 데워집니다: $\frac{dT}{dt} = -k(T - T_{\text{air}}) + q$, $k = 0.1$/분, $T_{\text{air}} = 35$°C, $q = 1.5$°C/분, 처음 $T(0) = 35$°C. rk4_step과 solve로 h = 5분, 0~30분을 풀고 정확한 해 $T(t) = 50 - 15e^{-kt}$와 비교하세요. 배터리 상한 55°C를 넘나요?
2. 같은 식을 오일러로 h = 25분, 0~100분 풀어 보세요. 무슨 일이 일어나고, 안정 조건 $h < 2/k$로 어떻게 설명되나요?
3. 바람 8 m/s일 때 착지점을 RK4 h = 0.1과 h = 0.05로 구하고, 두 값이 몇 자리까지 같은지 확인하세요.

풀이

>>> k, T_air, q = 0.1, 35.0, 1.5
>>> dT = lambda t, T: -k * (T - T_air) + q
>>> t, T = solve(rk4_step, dT, 35.0, 30, 5)
>>> np.round(T, 2)
array([35.  , 40.9 , 44.48, 46.65, 47.97, 48.77, 49.25])
>>> print(f"{np.abs(T - (50 - 15 * np.exp(-k * t))).max():.1e}")
4.4e-03
>>> t, T = solve(euler_step, dT, 35.0, 100, 25)
>>> np.round(T, 1)
array([ 35. ,  72.5,  16.2, 100.6, -25.9])
>>> drop(rk4_step, 0.1, 8.0), drop(rk4_step, 0.05, 8.0)
((6.793, 44.3), (6.793, 44.3))

1번은 30분 뒤 49.25°C, 오래 날아도 $T_{\text{air}} + q/k = 50$°C에 다가갈 뿐이라 55°C를 넘지 않습니다. 걸음 5분짜리 RK4가 정확한 해와 0.005°C 안에서 맞습니다.

2번은 $1 - hk = 1 - 2.5 = -1.5$라서 50°C와의 차이가 걸음마다 −1.5배가 됩니다. 72.5°C, 16.2°C, 100.6°C로 출렁이며 커지다 영하로 떨어집니다. 실제 배터리가 이렇게 될 리는 없고 계산이 불안정한 것입니다. $h < 2/k = 20$분이어야 하고, 출렁임 없이 가라앉게 하려면 $h < 1/k = 10$분이어야 합니다.

3번은 두 걸음 크기에서 착지 시각 6.793초, 착지 x 44.3 m가 소수 셋째 자리까지 같습니다. 바람 8 m/s면 지상 팀은 바람 방향으로 44 m쯤 떨어진 곳에서 기다리면 됩니다.

강의를 마치며

12교시 동안 하늘상자 관제실의 문제를 따라 수치해석 입문서의 큰 갈래를 모두 지났습니다. 13교시에서 컴퓨터 속 숫자와 근사의 대가(반올림 오차, 절단 오차)를 재는 법을 익혔고, 46교시에서 방정식의 근과 최적점을 반복으로 찾았습니다. 78교시에서는 연립방정식을 소거와 반복으로 풀며 조건수를 만났고, 911교시에서는 데이터에 선을 맞추고, 빈 곳을 채우고, 넓이를 쟀습니다. 그리고 오늘, 1교시에 오일러 한 걸음으로 시작한 낙하산 상자를 RK4와 벡터 상태로 다시 풀었습니다. 방법은 모두 달랐지만 할 일은 늘 세 가지였습니다. 식을 세우고, 덧셈과 곱셈의 반복으로 근사해를 구하고, 그 답이 얼마나 믿을 만한지 숫자로 말하기. 다음에 처음 보는 문제를 만나도 이 세 단계와 "h를 바꿔 다시 풀어 보기" 하나만 들고 가면 됩니다.

그날 오전, 상자는 바람 방향으로 27 m쯤 떨어진 풀밭에 내려앉았습니다. 미리 그 자리에 서 있던 지상 팀이 손을 흔들었습니다.

2주 뒤, 관제실에 새 신입이 들어왔습니다. 첫날부터 민재 팀장이 태블릿을 흔들며 뛰어 들어옵니다. "충전소 대기 줄이 이상해요!" 신입이 공식부터 찾는 사이, 서윤은 노트북을 열며 말합니다. "공식이 없으면 한 걸음씩 가면 돼요. 그 전에, 오차부터 정하죠." 창가 자리의 오 박사가 커피잔 너머로 슬쩍 웃습니다.