얼마나 날았나: 수치 적분

입사 둘째 주가 끝나 갈 무렵, 민재 팀장이 비행 기록을 출력한 종이 두 장을 들고 들어옵니다. "어젯밤 4호기가 산마을 왕복을 했는데, 착륙할 때 배터리 표시가 25%였어요. 규정이 20% 이상 남기고 착륙하기니까 통과는 했죠. 그런데 정비 팀이 그 표시를 못 믿겠대요. 지난주에 계기가 한 번 틀렸거든요."

기록에는 이륙부터 착륙까지 960초 동안 15초마다 속도(m/s)와 전력(W)이 찍혀 있습니다. 서윤이 먼저 말합니다. "거리는 속도 곱하기 시간이잖아요. 속도가 17 m/s쯤이니까 곱하기 960초 하면…." 민재가 고개를 젓습니다. "중간에 1분 넘게 제자리에 떠 있었어요. 상자를 내려놓느라고요. 오는 길엔 맞바람이라 더 느렸고요."

오 박사가 기록을 받아 듭니다. "속도가 계속 바뀌면 곱셈 한 번으로는 안 돼. 그래프를 그리고 그 아래 넓이를 재야지. 넓이를 잘게 썰어서 더하는 게 수치 적분이야. 전력도 마찬가지고. 전력 그래프 아래 넓이가 쓴 에너지야."

이번 교시는 9교시(회귀)와 10교시(보간)에 이어 데이터를 다루는 세 번째 도구입니다. 표본 몇 개만 있는 곡선에서 넓이를 얼마나 정확하게 잴 수 있는지, 조각의 모양을 바꾸면 오차가 얼마나 빨리 줄어드는지 숫자로 확인합니다.

11-1. 거리는 그래프 아래 넓이 (8분)

속도가 일정하면 거리 = 속도 × 시간입니다. 속도-시간 그래프에서는 높이 × 폭, 즉 직사각형의 넓이입니다. 속도가 바뀌면 그래프가 구불구불해지지만 생각은 같습니다. 시간을 잘게 나누면 한 조각 안에서는 속도가 거의 일정하니, 조각마다 직사각형 넓이를 구해 더하면 됩니다. 조각을 한없이 잘게 할 때의 넓이를 정적분(definite integral) 이라고 하고 이렇게 적습니다.

$\int$은 합(Sum)의 S를 길게 늘인 기호이고, $dt$는 "아주 얇은 시간 조각"입니다. a초부터 b초까지 $v(t) \times dt$를 모두 더한다고 읽으면 됩니다. 전력 $P(t)$를 같은 방법으로 더하면 에너지가 됩니다.

조각 $n$개로 나누면 폭은 $h = (b - a)/n$입니다. 직사각형의 높이를 어디서 재느냐에 따라 이름이 붙습니다. 모두 리만 합(Riemann sum) 이라고 부릅니다.

정확한 답을 아는 함수로 먼저 연습해 봅시다. 이륙 직후 드론의 속도를 $v(t) = 20(1 - e^{-t/4})$ m/s로 모델링하면, 처음 16초 동안 날아간 거리는 미적분으로 정확히 구할 수 있습니다.

조각 4개($h = 4$초)로 세 가지 합을 구해 정확한 값과 견줍니다.

>>> f = lambda t: 20 * (1 - np.exp(-t / 4))
>>> a, b, n = 0, 16, 4
>>> h = (b - a) / n
>>> t = np.linspace(a, b, n + 1)
>>> left = h * f(t[:-1]).sum()
>>> right = h * f(t[1:]).sum()
>>> mid = h * f(t[:-1] + h / 2).sum()
>>> exact = 320 - 80 * (1 - np.exp(-4))
>>> np.round([left, right, mid, exact], 2)
array([195.76, 274.29, 244.64, 241.47])

속도가 계속 오르는 구간이라 왼쪽 직사각형은 넓이를 빼먹고(−45.7 m), 오른쪽은 더 칩니다(+32.8 m). 중점은 +3.2 m로 훨씬 낫습니다. 조각 수는 같은데 높이를 재는 곳만 바꿔도 오차가 열 배 넘게 차이 납니다. 조각의 모양을 잘 고르는 것이 이번 교시의 주제입니다.

11-2. 사다리꼴 공식 (10분)

비행 기록에는 15초마다의 값만 있어서 조각 한가운데의 값(중점)을 모릅니다. 가진 것은 조각 양 끝의 값입니다. 그렇다면 두 끝점을 직선으로 이어 사다리꼴을 만드는 것이 자연스럽습니다. 사다리꼴 하나의 넓이는 (윗변 + 아랫변) × 높이 ÷ 2, 그래프로는 이렇습니다.

조각 $n$개의 넓이를 모두 더하면 안쪽 점은 이웃한 사다리꼴 두 개에 한 번씩, 모두 두 번 쓰이고 양 끝점은 한 번씩만 쓰입니다. 그래서 복합 사다리꼴 공식(composite trapezoidal rule) 은 이렇게 정리됩니다.

식을 잘 보면 왼쪽 직사각형 합과 오른쪽 직사각형 합의 평균입니다. 하나는 모자라고 하나는 넘치던 두 합이 서로의 오차를 상당 부분 지웁니다.

손으로 한 번 해 봅시다. 4호기가 이륙한 첫 1분의 기록입니다.

같은 숫자로 왼쪽 직사각형은 $15 \times 19.4 = 291$ m, 오른쪽은 $15 \times 34.5 = 517.5$ m이고, 둘의 평균이 정확히 404.25 m입니다. 파이썬으로 옮기면 한 줄입니다. 배열 슬라이스 y[1:-1]이 양 끝을 뺀 안쪽 점들입니다.

>>> def trapezoid(y, h):
...     return h * (y[0] / 2 + y[1:-1].sum() + y[-1] / 2)
>>> v4 = np.array([0.0, 2.2, 6.4, 10.8, 15.1])
>>> round(trapezoid(v4, 15), 2)
np.float64(404.25)

이제 960초 기록 전체(점 65개, 조각 64개)입니다. 속도와 전력을 함께 적어 둡니다.

>>> v_log = np.array([0.0, 2.2, 6.4, 10.8, 15.1, 16.8, 17.2, 17.0, 17.3, 16.4, 17.5, 17.0, 17.2,
...     17.0, 16.9, 17.2, 17.3, 16.9, 16.9, 17.2, 16.7, 16.5, 17.1, 16.8, 13.7, 8.2, 2.5, 0.0, 0.0,
...     0.0, 0.0, 0.0, 1.2, 5.1, 9.3, 12.4, 14.2, 14.4, 13.9, 13.7, 14.1, 14.1, 14.4, 13.6, 13.8, 13.7,
...     13.4, 14.0, 14.2, 13.7, 14.5, 14.3, 14.2, 14.1, 14.3, 13.5, 14.2, 14.2, 13.4, 12.7, 9.0, 5.2,
...     1.3, 0.0, 0.0])
>>> P_log = np.array([1450, 1510, 1500, 1510, 1430, 1340, 1180, 1030, 1000, 1010, 1020, 940, 1000,
...     990, 980, 950, 1010, 950, 990, 1010, 940, 970, 950, 990, 1090, 1260, 1300, 1370, 1400, 1420,
...     1390, 1360, 1360, 1300, 1210, 1120, 1090, 1070, 1060, 1050, 1030, 1050, 1090, 1060, 1020,
...     1090, 1070, 1110, 1060, 1050, 1020, 1090, 1120, 1040, 1040, 1070, 1040, 1050, 1050, 1090,
...     1180, 1220, 1270, 1240, 0], dtype=float)
>>> t_log = np.arange(65) * 15.0
>>> len(v_log), t_log[-1]
(65, np.float64(960.0))
>>> round(trapezoid(v_log, 15), 1)
np.float64(11398.5)

4호기는 약 11.4 km를 날았습니다. 같은 일을 하는 NumPy 함수가 np.trapezoid입니다. 간격을 dx=로 주거나, 시각 배열을 두 번째 인자로 줍니다.

>>> round(np.trapezoid(v_log, dx=15), 1)
np.float64(11398.5)
>>> round(np.trapezoid(v_log, t_log), 1)
np.float64(11398.5)

옛 자료에서 np.trapz라는 이름을 자주 보게 됩니다. 같은 함수의 옛 이름인데, NumPy 2.0에서 np.trapezoid로 바뀌었고 지금 쓰는 NumPy에는 아예 없습니다. 예전 코드를 옮길 때는 이름만 바꾸면 됩니다.

11-3. 포물선으로 잇기: 심프슨 1/3 공식 (8분)

사다리꼴은 두 점을 직선으로 잇습니다. 곡선이 휘어 있으면 직선과 곡선 사이에 틈이 남습니다. 점을 세 개씩 묶으면 그 셋을 지나는 포물선(2차 곡선)이 딱 하나 있고, 포물선은 휘어진 모양을 따라갈 수 있습니다. 간격이 $h$로 같은 세 점 $f_0, f_1, f_2$를 지나는 포물선 아래 넓이는 놀랍게도 간단합니다.

가운데 점에 가중치 4, 양 끝에 1입니다. 조각 두 개씩을 한 묶음으로 써서 이어 붙이면 복합 심프슨 1/3 공식(Simpson's 1/3 rule) 이 됩니다.

홀수 번째 점은 4, 안쪽의 짝수 번째 점은 2(두 묶음의 경계라 두 번 쓰임), 양 끝은 1입니다. 두 조각씩 묶으므로 조각 수 $n$이 짝수여야 합니다. 이름의 1/3은 앞의 $h/3$에서 왔습니다.

이륙 첫 1분(조각 4개)으로 손 계산을 해 봅니다.

사다리꼴(404.25 m)보다 조금 작습니다. 가속이 점점 줄어드는(위로 볼록한) 구간이라 직선이 곡선보다 높게 지나간 부분을 포물선이 덜어 낸 것입니다. 코드에서는 y[1:-1:2]가 홀수 번째 점, y[2:-1:2]가 안쪽의 짝수 번째 점입니다.

>>> def simpson(y, h):
...     if (len(y) - 1) % 2:
...         raise ValueError("심프슨 1/3 공식은 조각 수가 짝수여야 합니다")
...     return h / 3 * (y[0] + 4 * y[1:-1:2].sum() + 2 * y[2:-1:2].sum() + y[-1])
>>> round(simpson(v4, 15), 2)
np.float64(399.5)
>>> round(simpson(f(t), h), 2)
np.float64(241.08)
>>> round(simpson(v_log, 15), 1)
np.float64(11386.0)

모델 함수에서는 조각 4개만으로 241.08 m, 정확한 값 241.47 m에 0.39 m까지 다가갔습니다. 같은 점 다섯 개로 사다리꼴은 235.03 m(6.4 m 모자람)였으니 오차가 16분의 1로 줄었습니다. 비행 기록 전체에서는 사다리꼴 11,398.5 m와 심프슨 11,386.0 m가 12.5 m(0.1%) 차이로 거의 같습니다.

아래 실습대에서 조각을 직접 채워 보세요. 회색 선이 참 곡선(또는 비행 기록), 남색 조각이 왼쪽 직사각형·사다리꼴, 청록 조각이 심프슨의 포물선입니다. n 슬라이더로 조각 수를 바꾸고, 'n 두 배씩 늘리기'를 누르면 n이 1 → 2 → 4 → … → 64로 두 배씩 늘며 조각이 촘촘해집니다. 아래 작은 그래프는 n에 따른 오차를 로그-로그 눈금으로 그린 것입니다. 남색 선이 사다리꼴, 청록 선이 심프슨, 주황 점이 지금 보고 있는 계산입니다.

11-4. 오차 차수: h²와 h⁴ (8분)

실습대의 작은 그래프에서 두 선은 모두 곧은 직선이고, 심프슨 선이 두 배쯤 가파릅니다. 로그-로그 그래프에서 직선이라는 것은 오차가 $h$의 거듭제곱에 비례한다는 뜻이고, 기울기가 그 지수입니다. 3교시의 테일러 급수로 유도한 결과만 적으면 이렇습니다.

$\xi$(크사이)는 구간 안 어딘가의 한 점입니다. 사다리꼴은 $O(h^2)$, 심프슨은 $O(h^4)$입니다. $f''$이 들어 있다는 것은 "곡선이 휜 정도만큼 틀린다"는 뜻입니다. 그래서 직선(1차식)은 사다리꼴이 정확히 맞히고, 심프슨은 3차식까지 정확히 맞힙니다.

숫자로 확인합니다. 정확한 값을 아는 가속 모델에서 조각 수를 두 배씩 늘립니다.

>>> ns = np.array([2, 4, 8, 16, 32, 64])
>>> hs = 16 / ns
>>> et = np.array([abs(trapezoid(f(np.linspace(0, 16, k + 1)), 16 / k) - exact) for k in ns])
>>> es = np.array([abs(simpson(f(np.linspace(0, 16, k + 1)), 16 / k) - exact) for k in ns])
>>> for k, a, b in zip(ns, et, es):
...     print(f"{k:3d}  {a:.3e}  {b:.3e}")
  2  2.458e+01  4.647e+00
  4  6.438e+00  3.893e-01
  8  1.629e+00  2.648e-02
 16  4.086e-01  1.692e-03
 32  1.022e-01  1.063e-04
 64  2.556e-02  6.654e-06
>>> np.round(et[:-1] / et[1:], 2)
array([3.82, 3.95, 3.99, 4.  , 4.  ])
>>> np.round(es[:-1] / es[1:], 1)
array([11.9, 14.7, 15.7, 15.9, 16. ])

조각을 두 배로 늘리면 사다리꼴 오차는 1/4, 심프슨 오차는 1/16에 다가갑니다. $2^2 = 4$, $2^4 = 16$입니다. 처음 몇 단계의 비율이 4와 16보다 작은 것은 $h$가 커서 아직 "$h$가 작을 때"의 규칙이 들어맞지 않기 때문입니다. 로그-로그 직선의 기울기는 9교시의 np.polyfit으로 잽니다. 심프슨은 $h$가 큰 앞의 두 점을 빼고 잽니다.

>>> round(np.polyfit(np.log(hs), np.log(et), 1)[0], 2)
np.float64(1.98)
>>> round(np.polyfit(np.log(hs[2:]), np.log(es[2:]), 1)[0], 2)
np.float64(3.99)

실무의 결론은 이렇습니다. 이 예에서 상대 오차를 백만분의 1 아래로 내리려면 사다리꼴은 조각이 660개쯤 필요하지만 심프슨은 28개면 됩니다. 점 하나에 드는 계산(함수 한 번)은 두 방법이 같습니다.

리처드슨 외삽(Richardson extrapolation). 사다리꼴 오차가 $h^2$에 비례한다는 사실을 거꾸로 이용할 수 있습니다. $T_n$의 오차를 $E$라 하면 $T_{2n}$의 오차는 약 $E/4$이므로, $\frac{4T_{2n} - T_n}{3}$에서는 $h^2$ 오차가 지워집니다. 계산해 보면 이것이 정확히 심프슨 공식입니다. 같은 일을 한 번 더 하면($h^4$ 오차 지우기) 더 높은 차수가 나오고, 이를 표로 반복하는 방법이 롬베르크 적분(Romberg integration) 입니다. 오차의 모양을 알면 오차를 지울 수 있다는 것, 수치해석에서 여러 번 다시 나오는 생각입니다.

>>> T4 = trapezoid(f(np.linspace(0, 16, 5)), 4)
>>> T8 = trapezoid(f(np.linspace(0, 16, 9)), 2)
>>> round((4 * T8 - T4) / 3, 6), round(simpson(f(np.linspace(0, 16, 9)), 2), 6)
(np.float64(241.438773), np.float64(241.438773))

11-5. 진짜 기록 다루기: 고르지 않은 간격, 누적 적분, 단위 (8분)

고르지 않은 간격. 정리된 15초 간격 기록 말고, 기체가 무선으로 보낸 원본 기록은 통신이 될 때마다 찍혀서 간격이 들쭉날쭉합니다. 사다리꼴은 조각마다 폭이 달라도 그대로 쓸 수 있습니다. 조각 $i$의 넓이가 $\frac{t_{i+1} - t_i}{2}(f_i + f_{i+1})$일 뿐입니다. np.trapezoid에 시각 배열을 x=로 주면 이 계산을 합니다.

>>> t_raw = np.array([0.0, 7.2, 15.0, 21.9, 30.0, 38.6, 45.0, 52.3, 60.0])
>>> v_raw = np.array([0.0, 0.8, 2.2, 4.1, 6.4, 9.1, 10.8, 13.2, 15.1])
>>> round(np.trapezoid(v_raw, x=t_raw), 2)
np.float64(405.72)
>>> round((np.diff(t_raw) * (v_raw[1:] + v_raw[:-1]) / 2).sum(), 2)
np.float64(405.72)
>>> round(np.trapezoid(v_raw), 2)
np.float64(54.15)

간격이 들쭉날쭉해도 15초 간격 기록으로 구한 404.25 m와 거의 같은 405.72 m가 나옵니다. 그런데 마지막 줄처럼 x=를 빼먹으면 NumPy는 간격을 1로 여기고 54.15를 돌려줍니다. 오류 없이 엉뚱한 값이 나오니 특히 조심해야 합니다. 심프슨 공식은 간격이 같다는 가정으로 만든 식이라 이런 기록에 그대로 쓸 수 없습니다. 원본 기록에는 사다리꼴, 또는 10교시처럼 먼저 고른 간격으로 보간한 뒤 심프슨을 씁니다.

누적 적분. "출발한 뒤 지금까지 몇 m 왔나"를 매 시각마다 알고 싶으면 조각 넓이를 차례로 누적하면 됩니다. np.cumsum은 앞에서부터의 누적 합입니다. 맨 앞에 0을 붙여 시각 배열과 길이를 맞춥니다.

>>> pieces = np.diff(t_log) * (v_log[1:] + v_log[:-1]) / 2
>>> dist = np.concatenate(([0.0], np.cumsum(pieces)))
>>> dist.shape
(65,)
>>> np.round(dist[[4, 27, 64]], 1)
array([  404.3,  5727. , 11398.5])

405초(27번째 기록)에 상자를 내려놓을 지점에 도착했고 그때까지 5,727 m, 돌아오는 길이 5,671.5 m로 왕복 거리가 거의 같습니다. 기록이 앞뒤가 맞는다는 좋은 신호입니다.

단위. 전력 W는 1초에 쓰는 에너지(J/s)입니다. W를 초 단위로 적분하면 W·s, 곧 J(줄)이 나옵니다. 배터리 용량은 Wh(와트시)로 적으니 1 Wh = 1 W × 3600 s = 3600 J로 나눠야 합니다.

>>> E_J = np.trapezoid(P_log, t_log)
>>> E_Wh = E_J / 3600
>>> round(E_Wh, 1)
np.float64(299.7)
>>> left = 400 - E_Wh
>>> round(left, 1), round(left / 400 * 100, 1)
(np.float64(100.3), np.float64(25.1))
>>> round((left - 80) / 1380 * 3600)
53

쓴 에너지는 약 299.7 Wh, 남은 것은 100.3 Wh(25.1%)입니다. 계기의 25%가 맞았습니다. 다만 규정선 20%(80 Wh)까지 여유는 20 Wh 남짓, 제자리 비행(약 1380 W)으로 치면 53초뿐입니다. 하강 지점에서 1분만 더 머물렀다면 규정을 어길 뻔했습니다.

실무에서는

SciPy의 scipy.integrate에 같은 일을 하는 함수가 모여 있습니다. 이 사이트의 서버에는 SciPy가 없어서 이름과 쓰임만 적습니다.

from scipy import integrate

# 함수를 알 때: 오차를 스스로 어림하며 필요한 곳만 촘촘하게 나눈다 (적응 구적법)
value, err_est = integrate.quad(lambda t: 20 * (1 - np.exp(-t / 4)), 0, 16)

# 표본만 있을 때: 심프슨 (간격이 고르지 않아도 x=로 처리, 조각 수가 홀수여도 동작)
energy_J = integrate.simpson(P_log, x=t_log)

# 누적 적분: 매 시각까지의 거리 (initial=0으로 길이를 맞춘다)
dist = integrate.cumulative_trapezoid(v_log, x=t_log, initial=0)

quad는 함수를 직접 부를 수 있을 때 가장 먼저 쓰는 도구입니다. 곡선이 급하게 변하는 곳만 골라 촘촘하게 나누고, 답과 함께 오차 어림값을 돌려줍니다. 기록 데이터처럼 표본만 있으면 simpson이나 trapezoid를 씁니다. 노이즈가 많은 기록이면 둘의 차이는 거의 없습니다.

책의 기호 읽는 법

자주 하는 실수 (3분)

점 개수와 조각 개수 헷갈리기. 기록 65개는 조각 64개입니다. 심프슨의 "짝수" 조건은 조각 수에 대한 것이라, 점 개수로는 홀수여야 합니다. 앞에서 만든 simpson이 len(y) - 1을 검사하는 이유입니다.

np.trapezoid(y)에 간격 빼먹기. 간격을 주지 않으면 1로 여깁니다. 15초 간격 기록이면 15배 작은 값이 오류 없이 나옵니다. 늘 np.trapezoid(y, t) 또는 dx=를 씁니다.

평균 × 전체 시간. v_log.mean() * 960은 양 끝점에 다른 점과 같은 무게를 줍니다. 끝점에 절반 무게를 주는 사다리꼴과 값이 다릅니다. 간격이 고르지 않으면 차이가 훨씬 커집니다.

>>> round(v_log.mean() * 960, 1), round(np.trapezoid(v_log, t_log), 1)
(np.float64(11223.1), np.float64(11398.5))

단위 빠뜨리기. W를 초로 적분한 값(J)을 Wh로 착각하면 1,078,875 Wh라는, 배터리 수천 개 분량의 답이 나옵니다.

연습 문제 (5분)

1. 전력 기록을 simpson으로 적분해 Wh로 바꾸세요. 사다리꼴 값(299.7 Wh)과의 차이가 400 Wh 배터리의 몇 %인지 구하고, 20 Wh 남짓한 여유와 견주어 보세요.
2. $\int_0^1 \frac{4}{1+x^2}dx = \pi$입니다. 조각 8개로 사다리꼴과 심프슨 값을 구하고 각각의 오차를 np.pi와 비교하세요.
3. 누적 적분으로 매 시각 배터리 잔량(%)을 구하고, 잔량이 처음으로 50% 아래로 내려간 기록 시각을 찾으세요.
4. 가속 모델에서 조각 8개로 중점 합과 사다리꼴 값의 오차를 구하세요. 두 오차의 부호와 크기에서 어떤 관계가 보이나요?

풀이

>>> E_simp = simpson(P_log, 15) / 3600
>>> round(E_simp, 2), round((E_simp - E_Wh) / 400 * 100, 2)
(np.float64(300.65), np.float64(0.24))
>>> g4 = lambda x: 4 / (1 + x**2)
>>> xs = np.linspace(0, 1, 9)
>>> print(f"{trapezoid(g4(xs), 1 / 8) - np.pi:.2e}  {simpson(g4(xs), 1 / 8) - np.pi:.2e}")
-2.60e-03  -1.51e-07
>>> E_cum = np.concatenate(([0.0], np.cumsum(np.diff(t_log) * (P_log[1:] + P_log[:-1]) / 2))) / 3600
>>> remain = (400 - E_cum) / 400 * 100
>>> k = np.argmax(remain < 50)
>>> k, t_log[k], round(remain[k], 1)
(np.int64(42), np.float64(630.0), np.float64(49.6))
>>> t8 = np.linspace(0, 16, 9)
>>> mid8 = 2 * f(t8[:-1] + 1).sum()
>>> print(f"{mid8 - exact:.4f}  {trapezoid(f(t8), 2) - exact:.4f}")
0.8121  -1.6294

1번은 약 0.97 Wh 차이로 배터리의 0.24%입니다. 적분 방법을 바꿔서 생기는 차이가 여유분(20 Wh)보다 스무 배 작으니, 이 판단에는 어느 방법을 써도 됩니다. 남은 불확실성은 공식이 아니라 전력 센서의 정확도에서 옵니다.

2번은 같은 점 아홉 개로 사다리꼴은 소수 둘째 자리, 심프슨은 여섯째 자리까지 맞습니다. 오차가 만 배 넘게 차이 납니다.

3번은 630초입니다. 돌아오는 길(맞바람)을 시작하고 얼마 지나지 않아 절반을 넘겼습니다. np.argmax(조건)은 조건이 처음 참이 되는 위치를 찾는 흔한 요령입니다(조건이 한 번도 참이 아니면 0을 돌려주니 확인이 필요합니다).

4번은 중점 오차가 사다리꼴 오차의 약 −1/2입니다. 위로 볼록한 곡선에서 사다리꼴은 모자라고 중점은 넘칩니다. 그래서 $\frac{2M + T}{3}$으로 섞으면 두 오차가 지워집니다. 이 값은 중점까지 포함한 점 17개(조각 16개)로 계산한 심프슨 값과 정확히 같습니다.

정비 팀에 계산을 넘기고 돌아오니, 다음 날 아침 일기 예보가 떠 있습니다. 산마을 쪽에 초속 5~8 m의 바람. 민재 팀장이 한숨을 쉽니다. "내일 아침에 큰 낙하산으로 비상 투하 연습이 잡혀 있는데, 바람이 불면 상자가 어디 떨어질지 모르겠어요." 마지막 교시에는 1교시의 낙하산 상자로 돌아가, 위치와 속도를 함께 푸는 미분방정식과 오일러보다 훨씬 정확한 RK4를 배웁니다.