근사의 대가: 테일러 급수와 절단 오차

셋째 날 아침, 민재 팀장이 손바닥만 한 기판을 들고 왔습니다. "새로 들인 소형 드론의 비행 제어 칩이에요. 배터리 남은 에너지를 1초에 50번씩 계산해야 하는데, 제조사 말로는 이 칩에서 exp() 한 번이 곱셈 수백 번만큼 오래 걸린대요. 그리고 관제 화면에 '지금 1분에 몇 Wh씩 닳는지'도 띄워 달래요."

배터리 제조사가 준 모델은 드론 한 대의 400 Wh 배터리가 비행하는 동안 이렇게 줄어든다는 경험식입니다. $t$는 이륙 뒤 지난 분입니다.

$$E(t) = 400, e^{-0.012, t}\ \text{Wh}$$

서윤이 말합니다. "기울기는 미분하면 되잖아요. $E'(t) = -4.8, e^{-0.012 t}$. 근데 exp를 못 쓰면 이것도 못 쓰네요." 오 박사가 덧붙입니다. "그리고 실제 비행에서는 모델이 아니라 1분마다 찍힌 숫자가 들어와. 숫자에는 미분 공식이 없지. 두 문제 모두 같은 도구로 풀려. 함수를 그 점 근처에서 다항식으로 바꿔 보는 거야. 대신 공짜는 아니야. 근사에는 늘 대가가 있어."

1교시 오차 표에서 절단 오차(truncation error) 를 3교시에 다룬다고 했습니다. 무한한 계산을 유한한 곳에서 잘라 생기는 오차입니다. 이번 교시에는 테일러 다항식으로 그 오차의 크기를 재고, 같은 생각으로 수치미분 공식을 만들고, 마지막으로 2교시의 반올림 오차와 절단 오차가 줄다리기하는 모습을 봅니다.

3-1. 함수를 다항식으로: 테일러 다항식 (11분)

다항식은 덧셈과 곱셈만으로 계산됩니다. 그래서 어떤 함수든 다항식으로 바꿀 수 있다면 작은 칩에서도 계산할 수 있습니다. 문제는 어떤 다항식이냐입니다.

한 점 $a$ 근처에서 함수를 흉내 내는 다항식을 단계별로 만들어 봅시다.

이렇게 $a$에서 도함수 $n$개까지 함수와 똑같이 맞춘 다항식이 테일러 다항식(Taylor polynomial) 입니다.

$f^{(k)}$는 $k$번 미분한 함수, $k! = 1 \times 2 \times \cdots \times k$는 계승(factorial) 입니다($0! = 1$). $k!$로 나누는 이유는, $(x-a)^k$를 $k$번 미분하면 $k!$이 튀어나오기 때문입니다. 나눠 두어야 $P_n$의 $k$번째 도함수가 정확히 $f^{(k)}(a)$가 됩니다. 항의 개수를 끝없이 늘린 것을 테일러 급수(Taylor series) 라고 합니다.

$e^x$는 몇 번 미분해도 $e^x$이고 $e^0 = 1$입니다. 그래서 $a = 0$에서 모든 도함수가 1이고, 계수는 $1/k!$뿐입니다.

손으로 한 번 계산해 봅시다. 30분 비행 뒤 $x = -0.012 \times 30 = -0.36$입니다. 3차까지 더하면 $1 - 0.36 + 0.0648 - 0.007776 = 0.697024$입니다. 참값 $e^{-0.36} = 0.697676$과 0.00065 차이이고, 400 Wh를 곱하면 0.26 Wh입니다. 항을 늘리면 어떻게 될까요?

>>> from math import factorial
>>> def taylor_exp(x, n):
...     return sum(x**k / factorial(k) for k in range(n + 1))
>>> for n in [1, 2, 3, 4, 6, 8]:
...     approx = taylor_exp(-0.36, n)
...     print(n, round(approx, 10), f"{abs(approx - np.exp(-0.36)):.1e}")
1 0.64 5.8e-02
2 0.7048 7.1e-03
3 0.697024 6.5e-04
4 0.69772384 4.8e-05
6 0.6976764748 1.5e-07
8 0.6976763263 2.7e-10

항을 하나 더할 때마다 오차가 8배, 11배, 13배 … 로 줄고, 차수가 높을수록 더 크게 줄어듭니다. 계수의 분모 k!이 빠르게 커지기 때문입니다. 칩이 실제로 쓰는 범위는 이륙부터 60분까지, $x \in [-0.72, 0]$입니다. 그 범위 전체에서 가장 큰 오차를 Wh로 재 봅시다.

>>> t = np.linspace(0, 60, 61)
>>> x = -0.012 * t
>>> for n in [2, 4, 6]:
...     approx = sum(x**k / factorial(k) for k in range(n + 1))
...     print(n, f"{400 * np.abs(approx - np.exp(x)).max():.1e}")
2 2.1e+01
4 5.7e-01
6 7.3e-03

6차면 60분 내내 오차가 0.01 Wh 밑입니다. 곱셈 열 번 남짓이면 되니 exp() 한 번보다 훨씬 쌉니다. 칩에서는 $1 + x(1 + \frac{x}{2}(1 + \frac{x}{3}(\cdots)))$처럼 안쪽부터 묶어 계산하는 호너 방법(Horner's method) 을 쓰면 곱셈이 더 줄어듭니다.

아래 실습대에서 회색 곡선이 함수, 남색 곡선이 테일러 다항식입니다. 검은 손잡이(전개점 a)를 함수 위에서 끌면 다항식이 그 점을 기준으로 다시 만들어집니다. 호박색 손잡이(x)를 끌면 진홍 막대가 그 점에서의 오차를 보여 줍니다. '항 하나씩 더하기'의 재생을 누르면 차수가 하나씩 올라가며 다항식이 함수에 붙어 가고, 오른쪽 표에는 계수 cₖ가 k!로 나뉘며 빠르게 작아지는 모습이 보입니다.

3-2. 잘라 낸 값: 절단 오차와 O(hⁿ) (8분)

테일러 급수를 $n$차에서 자르면, 버린 항들의 합이 곧 오차입니다. 이것이 절단 오차입니다. 수학자 라그랑주는 버린 항 전체를 한 항으로 쓸 수 있다는 것을 보였습니다. 나머지 항(remainder) 이라고 부릅니다.

$\xi$(크시)가 정확히 어디인지는 모릅니다. 하지만 그 사이에서 $|f^{(n+1)}|$이 가장 클 때의 값으로 오차의 상한을 잴 수 있습니다. 말로 하면 이렇습니다. 오차는 버린 첫 항과 비슷한 크기다.

$e^x$에서 $x = -0.36$, $n = 3$이면 $\xi \le 0$이라 $e^{\xi} \le 1$이고, 상한은 $0.36^4 / 4! = 0.000700$입니다. 실제 오차 0.000652가 그 안에 들어옵니다.

나머지 항에서 더 중요한 것은 $(x-a)^{n+1}$입니다. 전개점에서의 거리를 $h = x - a$라고 두면, $n$차 다항식의 오차는 $h^{n+1}$에 비례합니다. 이것을 $O(h^{n+1})$ 이라고 적고 "오 에이치 엔 플러스 일"이라고 읽습니다. 뜻은 간단합니다. $h$를 반으로 줄이면 오차가 $2^{n+1}$분의 1로 준다. 2차 다항식으로 확인해 봅시다.

>>> for h in [0.4, 0.2, 0.1, 0.05]:
...     err = abs(taylor_exp(h, 2) - np.exp(h))
...     print(h, f"{err:.3e}")
0.4 1.182e-02
0.2 1.403e-03
0.1 1.709e-04
0.05 2.110e-05

$h$를 반으로 줄일 때마다 오차가 약 8분의 1($2^3$)이 됩니다. 2차 다항식의 오차가 $O(h^3)$이라는 뜻입니다. 1교시의 오일러 방법이 $O(h)$였던 것도 같은 언어입니다. 지수가 클수록 $h$를 조금만 줄여도 오차가 확 줄어듭니다.

테일러 다항식이 늘 좋아지기만 하지는 않습니다. $\ln x$를 $a = 1$에서 전개하면 $x$가 0과 2 사이에 있을 때만 항을 더할수록 참값에 다가갑니다. 그 밖에서는 항을 더할수록 오히려 멀어집니다. 급수가 함수에 다가가는 범위를 수렴 반지름(radius of convergence) 이라고 합니다. $e^x$, $\sin x$, $\cos x$는 수렴 반지름이 무한대라 어디서나 쓸 수 있지만, 실습대의 $\ln x$처럼 범위가 정해진 함수도 많습니다. 근사식을 쓸 때는 어디까지 써도 되는지를 함께 적어 두세요.

3-3. 기울기를 숫자로: 전진·후진·중앙 차분 (9분)

이제 기울기입니다. 미분의 정의는 "아주 작은 $h$만큼 갔을 때의 변화율"입니다. 컴퓨터는 극한을 계산할 수 없으니, 작은 $h$를 하나 정해 그대로 계산합니다. 어떤 공식이 얼마나 정확한지는 테일러 급수가 알려 줍니다. $x$에서 $h$만큼 떨어진 값을 전개하면 이렇습니다.

$f(x)$를 넘기고 $h$로 나누면 전진 차분(forward difference) 이 나옵니다.

뒤쪽 점 $x - h$를 쓰면 후진 차분(backward difference) $\dfrac{f(x) - f(x-h)}{h}$이고, 오차는 역시 $O(h)$입니다. 이번에는 $f(x+h)$와 $f(x-h)$의 전개식을 빼 봅시다. 짝수 차 항($f$, $f''$)이 서로 지워집니다.

이것이 중앙 차분(central difference) 입니다. 같은 함수 값 두 개를 쓰는데 오차가 $O(h^2)$입니다. $h$를 반으로 줄이면 전진·후진 차분은 오차가 반, 중앙 차분은 4분의 1이 됩니다. 배터리 모델로 확인해 봅시다. 30분 시점의 참 기울기는 $-4.8\,e^{-0.36} \approx -3.3488$ Wh/분입니다.

>>> E = lambda t: 400 * np.exp(-0.012 * t)
>>> exact = -4.8 * np.exp(-0.36)
>>> round(exact, 6)
np.float64(-3.348846)
>>> for h in [1.0, 0.5, 0.25]:
...     fwd = (E(30 + h) - E(30)) / h
...     bwd = (E(30) - E(30 - h)) / h
...     ctr = (E(30 + h) - E(30 - h)) / (2 * h)
...     print(h, f"{fwd - exact:+.2e} {bwd - exact:+.2e} {ctr - exact:+.2e}")
1.0 +2.00e-02 -2.02e-02 -8.04e-05
0.5 +1.00e-02 -1.01e-02 -2.01e-05
0.25 +5.02e-03 -5.03e-03 -5.02e-06

전진 차분의 오차는 $f''(30)/2 \times h = 0.0201h$로 식과 맞고, 후진 차분은 부호만 반대입니다. 중앙 차분은 1분 간격으로도 0.00008 Wh/분까지 맞고, 걸음을 반으로 줄이면 오차가 정확히 4분의 1이 됩니다. 1분마다 찍힌 기록에서 어느 시각의 기울기를 구한다면, 앞뒤 기록을 함께 쓰는 중앙 차분이 기본입니다.

3-4. h 딜레마: 절단 오차와 반올림 오차의 줄다리기 (10분)

그렇다면 $h$를 $10^{-15}$쯤으로 아주 작게 잡으면 끝일까요? $e^x$의 $x = 1$에서의 기울기(참값 $e$)로 $h$를 10배씩 줄여 봅시다.

>>> f, x0 = np.exp, 1.0
>>> for k in range(1, 16, 2):
...     h = 10.0**-k
...     ef = abs((f(x0 + h) - f(x0)) / h - np.e)
...     ec = abs((f(x0 + h) - f(x0 - h)) / (2 * h) - np.e)
...     print(f"{h:.0e}  {ef:.1e}  {ec:.1e}")
1e-01  1.4e-01  4.5e-03
1e-03  1.4e-03  4.5e-07
1e-05  1.4e-05  5.9e-11
1e-07  1.4e-07  5.9e-11
1e-09  2.2e-07  6.6e-09
1e-11  3.3e-05  1.0e-05
1e-13  4.6e-04  4.6e-04
1e-15  3.9e-01  1.7e-01

오차가 줄다가 어느 순간부터 다시 커집니다. 2교시의 자릿수 상쇄입니다. $h$가 작으면 $f(x+h)$와 $f(x)$가 거의 같은 수가 되고, 둘의 뺄셈에는 각각의 반올림 오차(약 $\varepsilon_{\text{mach}}|f|$)만 남습니다. 그것을 작은 $h$로 나누니 오차가 $\varepsilon_{\text{mach}}|f|/h$로 커집니다. 전진 차분의 전체 오차는 두 오차의 합입니다.

로그-로그 그래프에서 절단 오차는 오른쪽 위로 올라가는 직선(기울기 1), 반올림 오차는 왼쪽 위로 올라가는 직선(기울기 −1)입니다. 합은 U자가 되고, 두 직선이 만나는 곳이 바닥입니다. 합을 $h$로 미분해 0으로 두면 가장 좋은 $h$가 나옵니다. $|f| \approx |f''|$인 보통의 경우 이렇습니다.

>>> eps = np.finfo(float).eps
>>> np.sqrt(eps), eps ** (1 / 3)
(np.float64(1.4901161193847656e-08), np.float64(6.055454452393343e-06))

위 표의 숫자와 맞습니다. 전진 차분은 $h = 10^{-7} \sim 10^{-9}$에서, 중앙 차분은 $10^{-5}$ 근처에서 가장 정확하고, 그보다 작은 $h$는 손해입니다. 실무에서는 $x$의 크기도 반영해 $h = \sqrt{\varepsilon_{\text{mach}}} \cdot \max(1, |x|)$처럼 잡습니다.

아래 실습대는 이 줄다리기를 그대로 보여 줍니다. 가로축이 걸음 h(오른쪽이 큰 h), 세로축이 오차의 크기이고, 둘 다 로그 눈금입니다. 남색이 전진 차분, 청록이 중앙 차분의 오차이며, 호박색 세로선과 점이 슬라이더로 고른 지금의 h입니다. 왼쪽 비탈이 들쭉날쭉한 것은 브라우저가 float64로 실제 계산한 반올림 오차라서입니다. 'h를 줄여 가며 보기'를 누르면 h를 1e-1에서 1e-15까지 줄이며 두 곡선을 그려 나갑니다.

3-5. 기록된 숫자의 기울기: np.diff와 np.gradient (7분)

실제 비행에서는 함수가 아니라 1분마다 찍힌 숫자 배열이 들어옵니다. 배터리 관리 칩은 남은 에너지를 0.1 Wh 단위로 기록합니다. 처음 10분의 기록을 만들어 봅시다.

>>> t = np.arange(0, 11.0)
>>> E_log = np.round(400 * np.exp(-0.012 * t), 1)
>>> E_log
array([400. , 395.2, 390.5, 385.9, 381.3, 376.7, 372.2, 367.8, 363.4,
       359.1, 354.8])
>>> np.diff(E_log)
array([-4.8, -4.7, -4.6, -4.6, -4.6, -4.5, -4.4, -4.4, -4.3, -4.3])
>>> np.diff(E_log) / np.diff(t)
array([-4.8, -4.7, -4.6, -4.6, -4.6, -4.5, -4.4, -4.4, -4.3, -4.3])

np.diff(a)는 이웃한 원소의 차이 a[1:] - a[:-1]입니다. 시각의 차이로 나누면 구간마다의 평균 기울기, 즉 각 시각에서의 전진 차분이 됩니다. 결과가 하나 짧다는 점에 주의하세요. 점이 11개면 구간은 10개입니다.

np.gradient(y, t)는 길이를 그대로 지키면서, 안쪽 점에서는 중앙 차분을, 양 끝에서는 한쪽 차분을 씁니다. 간격이 고르지 않아도 됩니다.

>>> slope = np.gradient(E_log, t)
>>> slope
array([-4.8 , -4.75, -4.65, -4.6 , -4.6 , -4.55, -4.45, -4.4 , -4.35,
       -4.3 , -4.3 ])
>>> true_slope = -4.8 * np.exp(-0.012 * t)
>>> np.round(np.abs(slope - true_slope), 3)
array([0.   , 0.007, 0.036, 0.03 , 0.025, 0.03 , 0.017, 0.013, 0.011,
       0.009, 0.043])

모델의 기울기만 생각하면 1분 간격의 중앙 차분 오차는 0.0001 Wh/분 정도여야 합니다(3-3). 그런데 실제 오차는 0.04 Wh/분을 넘습니다. 기록이 0.1 Wh 단위로 반올림되어 있기 때문입니다. 두 기록의 반올림 오차(각각 최대 0.05 Wh)가 합쳐지면 최대 0.1 Wh이고, 이것이 2분(2h)으로 나뉘니 기울기에는 최대 0.05 Wh/분이 섞입니다. 3-4의 U자와 같은 이야기입니다. 측정값의 오차가 반올림 오차 역할을 합니다. 기록 간격을 1초로 줄이면 절단 오차는 사라지지만 이 오차는 60배로 커집니다. 잡음이 있는 기록에서는 간격을 너무 좁히지 말고, 여러 점을 매끈한 선으로 맞춘 뒤 그 선의 기울기를 쓰는 편이 낫습니다. 9교시에서 배터리 기록에 직선과 곡선을 맞추며 다시 다룹니다.

서윤은 관제 화면의 순간 소모율을 np.gradient로 계산하고, 표시를 소수 첫째 자리(0.1 Wh/분)까지만 하기로 했습니다. 기록의 정밀도보다 더 많은 자릿수를 보여 주는 것은 거짓말이기 때문입니다.

책의 기호 읽는 법

실무에서는

np.gradient(y, t)                     # 기록된 숫자의 기울기 (안쪽은 중앙 차분)
np.gradient(y, t, edge_order=2)       # 양 끝도 O(h²)로
np.diff(y) / np.diff(t)               # 구간마다의 평균 기울기 (길이가 하나 짧다)
np.polynomial.polynomial.polyval(x, c)    # 계수 c로 다항식 계산 (호너 방법)
np.expm1(x), np.log1p(x)              # 작은 x에서 정확한 eˣ - 1, ln(1 + x)
scipy.differentiate.derivative(f, x)  # 걸음을 줄여 가며 기울기와 오차 어림을 함께 (SciPy 1.15+)
scipy.optimize.approx_fprime(x, f)    # 여러 변수 함수의 기울기(그래디언트)를 차분으로
jax.grad(f)  # 또는 torch.autograd   # 자동 미분: 절단 오차 없이 기울기를 계산

머신러닝 라이브러리는 수치미분 대신 자동 미분(automatic differentiation) 을 씁니다. 계산 과정의 연산마다 미분 규칙을 적용해 기울기를 구하므로 h 딜레마가 없습니다. 그래도 자동 미분 코드를 시험할 때는 중앙 차분과 비교하는 것이 표준입니다(그래디언트 검사).

자주 하는 실수 (2분)

h는 작을수록 좋다고 믿기. h = 1e-15로 잡으면 오차가 오히려 커집니다. 전진 차분은 1e-8, 중앙 차분은 1e-5 근처가 기본값이고, 의심스러우면 h를 10배씩 바꿔 보세요.

np.gradient(y)에 시각을 빼먹기. 시각 배열을 주지 않으면 간격을 1로 칩니다. 2분 간격 기록이면 기울기가 두 배로 나옵니다. np.gradient(y, t)처럼 시각을 함께 넘기세요.

np.diff 결과를 원래 시각과 짝짓기. 길이가 하나 짧아서 t와 나란히 그리면 한 칸씩 어긋납니다. 구간의 가운데 시각 (t[1:] + t[:-1]) / 2와 짝지으세요.

근사식을 범위 밖에서 쓰기. 0~60분용으로 맞춘 다항식을 120분 비행에 쓰거나, 수렴 반지름 밖에서 항을 늘리면 오차가 폭발합니다.

연습 문제 (3분)

1. 칩에 넣을 $e^x$ 다항식의 차수를 정하려 합니다. $x \in [-0.72, 0]$에서 나머지 항의 상한 $0.72^{n+1}/(n+1)!$에 400을 곱한 값이 0.001 Wh 아래가 되는 가장 작은 $n$은 얼마일까요?
2. 2계 도함수(휘는 정도)도 차분으로 구할 수 있습니다. $\dfrac{f(x+h) - 2f(x) + f(x-h)}{h^2}$으로 배터리 모델의 $E''(30)$을 $h = 1$, $0.01$, $10^{-5}$에서 계산하고, 참값 $0.012^2 \times E(30)$과 비교하세요. $h = 10^{-5}$에서 무슨 일이 생기나요?
3. 2분마다 찍힌 기록에 np.gradient(E2)와 np.gradient(E2, t2)를 써 보고 차이를 확인하세요.

풀이

>>> for n in range(3, 9):
...     print(n, f"{400 * 0.72**(n + 1) / factorial(n + 1):.1e}")
3 4.5e+00
4 6.4e-01
5 7.7e-02
6 8.0e-03
7 7.2e-04
8 5.7e-05
>>> true2 = 0.012**2 * E(30)
>>> for h in [1.0, 0.01, 1e-5]:
...     d2 = (E(30 + h) - 2 * E(30) + E(30 - h)) / h**2
...     print(h, f"{d2:.6f}", f"{abs(d2 - true2):.1e}")
1.0 0.040187 4.8e-07
0.01 0.040186 5.4e-10
1e-05 0.040359 1.7e-04
>>> t2 = np.arange(0, 11.0, 2)
>>> E2 = 400 * np.exp(-0.012 * t2)
>>> np.round(np.gradient(E2), 3)
array([-9.486, -9.373, -9.151, -8.934, -8.722, -8.617])
>>> np.round(np.gradient(E2, t2), 3)
array([-4.743, -4.687, -4.575, -4.467, -4.361, -4.309])

1번: 상한이 0.001 Wh 아래로 내려가는 것은 $n = 7$입니다. 3-1에서 잰 6차의 실제 최대 오차 0.0073 Wh도 6차의 상한 0.0080 Wh와 비슷하니, 0.001 Wh를 지키려면 7차가 필요합니다. 상한은 실제 오차를 재지 않고도 차수를 정하게 해 줍니다.

2번: $h = 1$과 $0.01$에서는 잘 맞지만 $h = 10^{-5}$에서는 엉뚱한 값이 나옵니다. 반올림 오차가 $\varepsilon_{\text{mach}}|f|/h^2$로 커지기 때문입니다. $h^2$으로 나누니 1계 도함수보다 훨씬 빨리 망가집니다. 2계 차분의 가장 좋은 $h$는 $\varepsilon_{\text{mach}}^{1/4} \approx 10^{-4}$ 근처입니다.

3번: 시각을 주지 않으면 간격을 1분으로 쳐서 기울기가 정확히 두 배로 나옵니다.

퇴근 무렵, 지상 팀장이 1교시의 낙하산 상자 이야기를 다시 꺼냅니다. "상자가 떨어뜨린 지 4초 뒤에 10 m/s를 넘지 않게 낙하산을 다시 고르고 싶어요. 항력 계수가 얼마여야 하죠?" 이번에는 식에 값을 넣어 답을 내는 것이 아니라, 답이 되는 값을 찾아야 합니다. 다음 교시에는 f(x) = 0을 만드는 x를 구간을 반씩 줄여 가며 찾는 이분법과 가위치법을 배웁니다.