구간을 반으로: 이분법과 가위치법

서윤이 입사한 지 나흘째 아침입니다. 민재 팀장이 낙하산 견본 두 개를 들고 관제실로 들어옵니다. "지상 팀이 4초 뒤에 10 m/s 아래로 떨어져야 받을 수 있대요. 지금 낙하산은 15 m/s가 넘고요. 얼마나 큰 걸 주문해야 해요?"

서윤은 1교시에 쓴 공식을 화이트보드에 적습니다. $v(4) = \sqrt{gm/c_d}\,\tanh(\sqrt{g c_d/m}\cdot 4) = 10$. "이걸 $c_d$에 대해 풀면 되겠네요." 10분 뒤, 서윤은 펜을 내려놓습니다. $c_d$가 제곱근 밖에도, $\tanh$ 안에도 들어 있어서 한쪽으로 모아지지 않습니다.

오 박사가 웃습니다. "$c_d$를 꺼낼 수는 없어도, $c_d$를 넣어서 속도를 계산할 수는 있지? 넣어 보고, 너무 빠르면 크게, 너무 느리면 작게. 그 '넣어 보기'를 똑똑하게 하는 게 오늘 배울 근 찾기야."

이번 교시와 다음 교시는 방정식의 근 찾기(root finding) 입니다. 오늘은 근을 품은 구간을 정하고 그 구간을 좁혀 가는 방법 두 가지를 배웁니다. 느리지만 절대 놓치지 않는 방법들입니다.

4-1. 방정식을 f(x) = 0으로 바꾸기 (6분)

"4초 뒤 속도가 10이 되는 $c_d$"를 찾는 문제를 한 가지 모양으로 바꿉니다. 오른쪽의 10을 왼쪽으로 넘기면 됩니다.

$f(c_d)$는 "이 낙하산이면 4초 뒤 속도가 목표보다 얼마나 빠른가"입니다. 양수면 아직 빠르고, 음수면 목표보다 느립니다. $f(x) = 0$을 만드는 $x$를 함수의 근(root) 또는 영점(zero) 이라고 합니다. 어떤 방정식이든 한쪽으로 모으면 이 모양이 되므로, 근 찾기 방법 하나로 온갖 방정식을 풀 수 있습니다.

먼저 할 일은 언제나 그래프로 어림하기입니다. $c_d$를 0.25부터 2.5까지 0.25 간격으로 넣어 봅니다. np.arange(시작, 끝, 간격)은 끝 값을 빼고 만드니 끝을 2.75로 적었습니다.

>>> import numpy as np
>>> g, m = 9.81, 12
>>> def f(cd):
...     return np.sqrt(g * m / cd) * np.tanh(np.sqrt(g * cd / m) * 4) - 10
>>> round(f(0.48), 4)
np.float64(5.4532)
>>> cds = np.arange(0.25, 2.75, 0.25)
>>> cds
array([0.25, 0.5 , 0.75, 1.  , 1.25, 1.5 , 1.75, 2.  , 2.25, 2.5 ])
>>> np.round(f(cds), 2)
array([10.56,  5.16,  2.48,  0.83, -0.3 , -1.14, -1.8 , -2.33, -2.77,
       -3.14])

지금 낙하산(0.48)은 목표보다 5.45 m/s 빠릅니다. 값이 $c_d = 1$에서 0.83, $1.25$에서 −0.3이니 근은 1과 1.25 사이, 1.2보다 조금 작은 곳에 있습니다. 눈으로 얻은 답은 여기까지입니다. 지상 팀이 원하는 소수 넷째 자리까지 가려면 이 어림을 체계적으로 좁혀야 합니다.

4-2. 부호가 바뀌면 근이 있다: 중간값 정리 (5분)

종이에서 펜을 떼지 않고 x축 위의 점에서 x축 아래의 점까지 선을 그으면, 선은 반드시 x축을 한 번은 지납니다. 이것을 수학에서는 중간값 정리(intermediate value theorem) 라고 부릅니다.

부호가 다른지는 $f(a)\,f(b) < 0$으로 확인합니다. 근을 품은 이런 구간을 구간(bracket) 이라 부르고, 오늘 배우는 방법들을 구간법(bracketing method) 이라고 합니다. NumPy의 np.sign은 양수면 1, 음수면 −1, 0이면 0을 돌려줍니다. 부호가 바뀌는 칸을 찾을 때 편합니다.

>>> np.sign(f(cds))
array([ 1.,  1.,  1.,  1., -1., -1., -1., -1., -1., -1.])
>>> np.nonzero(np.diff(np.sign(f(cds))))
(array([3]),)
>>> cds[3], cds[4]
(np.float64(1.0), np.float64(1.25))
>>> f(1.0) * f(1.25) < 0
np.True_

np.diff는 이웃한 값의 차이입니다. 부호가 그대로면 0, 바뀌면 ±2가 되니 0이 아닌 칸이 부호가 바뀐 곳입니다. np.nonzero가 그 칸 번호 3을 찾아 주었습니다. 3번 칸과 4번 칸 사이, 즉 $[1.0,\ 1.25]$가 근을 품은 구간입니다.

4-3. 이분법: 반으로 자르기 (12분)

구간을 얻었으면 좁히기만 하면 됩니다. 가장 단순한 방법은 가운데를 잘라 보는 것입니다.

1. 가운데 점 $x = \dfrac{a + b}{2}$에서 $f(x)$를 계산한다.
2. $f(x)$의 부호가 $f(a)$와 같으면 근은 오른쪽 반에 있으니 $a \leftarrow x$, 다르면 왼쪽 반에 있으니 $b \leftarrow x$.
3. 구간이 충분히 좁아질 때까지 되풀이한다.

이것이 이분법(bisection method) 입니다. 숫자 맞히기 놀이에서 "50보다 커요?"부터 묻는 것과 같은 생각입니다. 한 번 물을 때마다 후보가 절반으로 줄어듭니다.

넓게 $[0.5,\ 2]$에서 시작해 손으로 세 번 잘라 봅시다. $f(0.5) = 5.1608$(양수), $f(2) = -2.3285$(음수)입니다.

세 번 만에 구간 폭이 1.5에서 0.1875로, 정확히 $1/2^3$이 되었습니다. 같은 일을 반복문으로 시킵니다.

>>> a, b = 0.5, 2.0
>>> for k in range(1, 7):
...     x = (a + b) / 2
...     print(k, a, b, x, round(f(x), 4))
...     if np.sign(f(x)) == np.sign(f(a)):
...         a = x
...     else:
...         b = x
1 0.5 2.0 1.25 -0.3015
2 0.5 1.25 0.875 1.5723
3 0.875 1.25 1.0625 0.5138
4 1.0625 1.25 1.15625 0.0817
5 1.15625 1.25 1.203125 -0.1154
6 1.15625 1.203125 1.1796875 -0.0183

부호를 f(x) * f(a) > 0이 아니라 np.sign끼리 비교한 데에는 이유가 있습니다. '자주 하는 실수'에서 봅니다. 여섯 번 잘랐을 때 근은 $[1.15625,\ 1.1796875]$ 안에 있습니다. 폭은 0.023이고, $x$가 1.18 근처에서 왔다 갔다 하며 근에 다가갑니다.

아래 실습대에서 x축 위의 두 손잡이 a, b를 끌어 구간을 정하세요. 남색 곡선이 f, 남색 띠가 지금 구간, 옅은 회색 띠가 바로 앞 단계의 구간입니다. 주황 점은 이번 단계에 새로 계산한 점(이분법은 가운데 점)이고, 축 위의 작은 회색 점들은 그 전에 계산한 점입니다. '재생'이나 '다음'을 누르면 한 단계씩 구간이 줄어들고, 아래 표에서 그 단계의 줄이 강조됩니다. 구간이 너무 좁아 보이지 않으면 '지금 구간 확대'를 켜세요.

4-4. 몇 번 잘라야 하나, 언제 멈추나 (8분)

이분법의 큰 장점은 필요한 반복 횟수를 미리 안다는 것입니다. 처음 폭이 $\Delta x_0 = b - a$이면 $n$번 자른 뒤의 폭은 $\Delta x_0 / 2^n$입니다. 이 폭을 허용 오차 $E$ 아래로 만들려면 이렇게 됩니다.

$[0.5,\ 2]$에서 폭 $10^{-4}$를 원하면 $\log_2(1.5 / 10^{-4}) \approx 13.9$이니 14번입니다. 함수가 아무리 꼬여 있어도, 연속이고 부호만 바뀌면 이 횟수가 지켜집니다. 그 대신 자릿수 하나(오차 1/10)를 더 얻으려면 $\log_2 10 \approx 3.3$번씩 더 잘라야 합니다. 이 속도를 다음 교시의 뉴턴법과 견주게 됩니다.

>>> np.ceil(np.log2((2.0 - 0.5) / 1e-4))
np.float64(14.0)
>>> 1.5 / 2**14
9.1552734375e-05

언제 멈출까? 실제 문제에서는 참값을 모르니 참 오차로 멈출 수 없습니다. 1교시에서 "두 근삿값의 차이를 오차의 어림값으로 쓴다"고 했습니다. 여기서도 같습니다. 이번 $x$와 지난번 $x$의 차이를 이번 $x$로 나눈 상대 오차 어림(approximate relative error) 이 정한 기준 $\varepsilon_s$보다 작아지면 멈춥니다.

이 정지 조건을 넣어 이분법을 함수로 묶습니다. 첫 단계의 $x_{\text{old}}$로는 $a$를 씁니다.

>>> def bisect(f, a, b, es=1e-4, maxit=100):
...     if np.sign(f(a)) == np.sign(f(b)):
...         raise ValueError("f(a)와 f(b)의 부호가 같습니다")
...     x = a
...     for k in range(1, maxit + 1):
...         x_old = x
...         x = (a + b) / 2
...         ea = abs(x - x_old) / abs(x)
...         if np.sign(f(x)) == np.sign(f(a)):
...             a = x
...         else:
...             b = x
...         if ea < es:
...             break
...     return x, k
>>> x, k = bisect(f, 0.5, 2.0)
>>> round(x, 6), k
(1.175385, 14)
>>> root = 1.1753511424325662
>>> abs(x - root) / root
2.8399216713820995e-05

maxit는 무슨 일이 있어도 끝나도록 거는 안전장치입니다. root는 기준을 아주 작게 해서 미리 구해 둔 값으로, 참 오차를 재 보려고만 씁니다. 실습대와 같은 14번에 멈췄고, 참 상대 오차는 $2.8\times10^{-5}$로 기준 $10^{-4}$보다 작습니다. 이분법에서는 어림 오차 $\varepsilon_a$가 참 오차보다 크게 나오기 때문에, $\varepsilon_a$로 멈추면 안전한 쪽으로 멈춥니다.

기준을 100배씩 엄하게 하면 반복 횟수는 6~7번씩 늘어납니다. 자릿수 두 개에 $2 \times 3.3$번입니다.

>>> for es in [1e-2, 1e-4, 1e-6, 1e-8]:
...     x, k = bisect(f, 0.5, 2.0, es)
...     print(es, k, round(x, 9), f"{abs(x - root) / root:.1e}")
0.01 8 1.173828125 1.3e-03
0.0001 14 1.175384521 2.8e-05
1e-06 21 1.175350904 2.0e-07
1e-08 27 1.175351139 2.8e-09

4-5. 가위치법: 직선으로 어림하기 (8분)

이분법은 $f$의 크기를 쓰지 않고 부호만 봅니다. $f(a) = 5.16$, $f(b) = -2.33$이면 근은 $b$ 쪽에 더 가까울 것 같은데, 이분법은 무조건 한가운데를 고릅니다.

크기를 쓰는 방법이 가위치법(false position, 라틴어로 regula falsi) 입니다. 두 끝점 $(a, f(a))$와 $(b, f(b))$를 직선으로 잇고, 그 직선이 x축과 만나는 곳을 새 점으로 씁니다. 곡선을 직선으로 '거짓으로' 바꿔 친다고 해서 이런 이름이 붙었습니다. 두 직각삼각형의 비율을 맞추면 이런 식이 나옵니다.

$[0.5,\ 2]$에 넣으면 $x = 2 - \dfrac{-2.3285 \times (0.5 - 2)}{5.1608 + 2.3285} \approx 1.534$입니다. 첫걸음의 가운데 점 1.25보다 근(1.1754)에서 멀어 보이지만, 남는 구간을 고르는 방법은 이분법과 똑같습니다. 그래서 근은 언제나 구간 안에 남습니다.

>>> def false_position(f, a, b, es=1e-4, maxit=100):
...     if np.sign(f(a)) == np.sign(f(b)):
...         raise ValueError("f(a)와 f(b)의 부호가 같습니다")
...     x = a
...     for k in range(1, maxit + 1):
...         x_old = x
...         x = b - f(b) * (a - b) / (f(a) - f(b))
...         ea = abs(x - x_old) / abs(x)
...         if np.sign(f(x)) == np.sign(f(a)):
...             a = x
...         else:
...             b = x
...         if ea < es:
...             break
...     return x, k
>>> x, k = false_position(f, 0.5, 2.0)
>>> round(x, 6), k
(np.float64(1.175401), 12)

12번으로, 이분법(14번)보다 조금 낫습니다. 생각보다 크게 낫지는 않습니다. 실습대에서 '가위치법'을 누르고 재생해 보면 이유가 보입니다. 낙하산 곡선은 아래로 볼록해서 직선이 늘 근의 오른쪽에 떨어지고, 오른쪽 끝 b만 계속 움직이고 왼쪽 끝 a는 0.5에 그대로 남습니다. 구간 폭이 0으로 줄지 않습니다.

이 약점이 극단적으로 드러나는 함수가 $x^{10} - 1$입니다. 근은 1인데, $[0,\ 1.3]$에서 곡선이 왼쪽은 거의 평평하고 오른쪽은 가파르게 솟습니다.

>>> def p10(x):
...     return x**10 - 1
>>> bisect(p10, 0, 1.3)
(0.9999938964843751, 14)
>>> false_position(p10, 0, 1.3)
(0.9996886513360608, 39)
>>> a, b = 0, 1.3
>>> for k in range(1, 6):
...     x = b - p10(b) * (a - b) / (p10(a) - p10(b))
...     print(k, round(a, 4), b, round(x, 4))
...     if np.sign(p10(x)) == np.sign(p10(a)):
...         a = x
...     else:
...         b = x
1 0 1.3 0.0943
2 0.0943 1.3 0.1818
3 0.1818 1.3 0.2629
4 0.2629 1.3 0.3381
5 0.3381 1.3 0.4079

$b = 1.3$이 끝까지 움직이지 않아 직선이 매번 거의 같은 기울기로 그어지고, $x$는 왼쪽에서 0.07~0.09씩 기어갑니다. 가위치법은 39번이나 걸렸고, 그렇게 멈춘 답도 참값 1과 $3\times10^{-4}$나 다릅니다. 기준(0.01%)의 세 배나 되는 오차입니다. 한 걸음씩 조금만 움직이니 '두 근삿값의 차이'가 작게 나와 오차를 실제보다 작게 어림한 것입니다. 이분법은 이런 함수에서도 14번입니다.

실습대로 돌아가 'x¹⁰ − 1' 칩과 '가위치법' 칩을 눌러 보세요. 청록 점선이 양 끝을 잇는 직선입니다. 재생하면 주황 점이 왼쪽에서 조금씩 기어가고, 상태 문구가 "한쪽에 갇혔다"고 알려 줍니다. 표의 b 열이 1.3에서 바뀌지 않는 것도 확인하세요.

4-6. 구간법이 실패할 때, 그리고 답의 뜻 (6분)

구간법은 "연속이고 부호가 바뀐다"는 조건이 맞을 때만 약속을 지킵니다. 조건이 깨지는 세 가지 경우를 기억해 두세요.

>>> def touch(x):
...     return (x - 1)**2
>>> touch(0) * touch(2) > 0
True
>>> x, k = bisect(np.tan, 1, 2)
>>> round(x, 6), k
(1.570679, 13)
>>> round(np.pi / 2, 6)
1.570796
>>> round(np.tan(1), 2), round(np.tan(x))
(np.float64(1.56), 8502)

$\tan x$는 $\pi/2$에서 $+\infty$에서 $-\infty$로 건너뛰므로 근이 없는데도 이분법이 그곳으로 다가갑니다. $|f(x)|$가 처음 1.56에서 8502로 커졌다는 것이 신호입니다. 끝난 뒤 $f(x)$가 정말 0에 가까운지 꼭 확인하세요. 그래서 4-1처럼 그래프부터 그려 보는 습관이 중요합니다.

찾은 답의 뜻. 기준을 $10^{-8}$로 엄하게 해서 다시 구하면 $c_d \approx 1.1754$ kg/m입니다.

>>> cd_new = bisect(f, 0.5, 2.0, 1e-8)[0]
>>> round(cd_new, 4), round(cd_new / 0.48, 2)
(1.1754, 2.45)
>>> round(np.sqrt(g * m / cd_new), 2)
np.float64(10.01)
>>> rho, CD = 1.2, 1.5
>>> for c in [0.48, cd_new]:
...     A = 2 * c / (rho * CD)
...     print(round(A, 2), round(np.sqrt(4 * A / np.pi), 2))
0.53 0.82
1.31 1.29

항력 계수가 지금의 2.45배여야 합니다. 새 종단 속도는 10.01 m/s라서, 상자는 4초 만에 거의 끝 속도에 닿고 그 뒤로는 10 m/s 남짓으로 내려옵니다. 공학 책에서는 $c_d = \tfrac{1}{2}\rho\,C_D\,A$(공기 밀도 × 모양 계수 × 넓이)로 쪼개 씁니다. 둥근 낙하산의 $C_D$를 1.5쯤으로 잡으면 펼친 넓이가 0.53 m²에서 1.31 m²로, 지름이 0.82 m에서 1.29 m로 커져야 합니다. 민재 팀장이 주문할 숫자입니다.

무게가 바뀌면? 상자 무게별로 다시 풀어 봅니다. lambda로 무게를 고정한 함수를 그때그때 만들어 넘깁니다.

>>> def f_mass(cd, mass):
...     return np.sqrt(g * mass / cd) * np.tanh(np.sqrt(g * cd / mass) * 4) - 10
>>> for mass in [6, 12, 18, 24]:
...     c, k = bisect(lambda cd: f_mass(cd, mass), 0.1, 5, 1e-8)
...     print(mass, round(c, 4), round(c / mass, 4))
6 0.5877 0.0979
12 1.1754 0.0979
18 1.763 0.0979
24 2.3507 0.0979

$c_d / m$이 늘 0.0979입니다. 식을 다시 보면 $m$과 $c_d$가 언제나 $g m / c_d$, $g c_d / m$처럼 비율로만 나오기 때문입니다. 무게가 두 배면 낙하산도 두 배로 크게. 이렇게 답을 여러 조건에서 다시 구해 보면 숫자 뒤의 구조가 보입니다.

책의 기호 읽는 법

실무에서는

SciPy에 구간법이 모두 들어 있습니다. 이 사이트 서버에는 SciPy가 없으니 이름과 쓰임만 소개합니다.

from scipy import optimize
optimize.bisect(f, 0.5, 2.0, xtol=1e-12)            # 이분법
optimize.brentq(f, 0.5, 2.0)                        # 구간을 지키면서 빠른 방법을 섞은 것 (5교시)
optimize.root_scalar(f, bracket=[0.5, 2.0], method="brentq")

실무에서 구간이 있을 때는 거의 언제나 brentq를 씁니다. 이분법처럼 근을 놓치지 않으면서 훨씬 빠릅니다. 어떻게 그럴 수 있는지는 다음 교시에서 봅니다. 구간을 모를 때는 4-2처럼 np.sign과 np.diff로 부호가 바뀌는 칸을 먼저 찾습니다.

자주 하는 실수 (2분)

곱으로 부호 비교하기. f(a) * f(b) < 0은 두 값이 아주 작거나 아주 크면 틀립니다. 곱이 float64가 담을 수 있는 가장 작은 수보다 작아져 0이 되기 때문입니다(2교시의 언더플로).

>>> fa, fb = 1e-200, -1e-200
>>> fa * fb
-0.0
>>> fa * fb < 0
False
>>> np.sign(fa) != np.sign(fb)
np.True_

끝난 뒤 f(x) 확인하지 않기. 4-6의 $\tan x$처럼 근이 아닌 곳으로 수렴할 수 있습니다. 근이 0 근처일 때 상대 오차 쓰기. $x$로 나누므로 근이 0이면 기준이 영영 만족되지 않거나 0으로 나누게 됩니다. 그럴 때는 구간 폭 같은 절대 기준을 함께 둡니다.

연습 문제 (3분)

1. $[0.5,\ 2]$에서 시작한 이분법으로 구간 폭을 $10^{-6}$ 이하로 줄이려면 몇 번 잘라야 할까요?
2. 깨지기 쉬운 약 상자는 4초 뒤 8 m/s를 넘으면 안 됩니다. 12 kg 상자에 필요한 $c_d$를 bisect로 구하세요. 구간은 어떻게 정하면 될까요?
3. 20 kg 상자가 4초 뒤 10 m/s가 되는 $c_d$를 구하고, 4-6에서 찾은 비율로 예상한 값과 비교하세요.
4. 가위치법이 $x^{10} - 1$에서 느린 이유를 '어느 끝이 움직이지 않는가'로 설명해 보세요.

풀이

>>> int(np.ceil(np.log2(1.5 / 1e-6)))
21
>>> def f8(cd):
...     return np.sqrt(g * m / cd) * np.tanh(np.sqrt(g * cd / m) * 4) - 8
>>> f8(0.5) * f8(5) < 0
np.True_
>>> x, k = bisect(f8, 0.5, 5, 1e-6)
>>> round(x, 4), k
(1.839, 22)
>>> round(np.sqrt(g * m / x), 2)
np.float64(8.0)
>>> c20, k = bisect(lambda cd: f_mass(cd, 20), 0.1, 5, 1e-8)
>>> round(c20, 4), round(cd_new * 20 / 12, 4)
(1.9589, 1.9589)

1번은 $\log_2(1.5 / 10^{-6}) \approx 20.5$이니 21번입니다. 4-4의 표에서 기준 $10^{-6}$일 때도 21번이었습니다.

2번은 목표가 8이니 함수의 마지막 숫자만 바꿉니다. 구간은 넉넉히 [0.5, 5]로 잡고, 시작 전에 부호가 다른지 확인했습니다. $c_d \approx 1.839$로, 10 m/s일 때보다 1.56배 큰 낙하산이 필요합니다. 이번에도 종단 속도가 8.0 m/s로 목표와 거의 같습니다.

3번은 1.9589로, $1.1754 \times 20 / 12$와 같습니다. 4번은 곡선이 한쪽으로 휘어 있으면 양 끝을 잇는 직선이 늘 근의 같은 쪽에 떨어지기 때문입니다. 반대쪽 끝($b = 1.3$)이 바뀌지 않으니 직선의 한쪽 발이 묶여 있고, 구간 폭이 0으로 줄지 않습니다. 묶인 끝의 $f$ 값을 일부러 절반으로 줄여 직선을 기울이는 고친 가위치법(일리노이 방법)도 있습니다.

다음 날 아침, 민재 팀장이 오늘 밤 배송 목록을 내밉니다. "상자가 300개예요. 무게도, 받는 곳에서 허용하는 속도도 제각각이고요. 비행 전에 상자마다 낙하산 크기를 정해야 해요." 상자 하나에 20번씩 자르면 6000번입니다. 다음 교시에는 접선을 타고 근으로 뛰어들어, 대여섯 번 만에 같은 답을 얻는 방법을 배웁니다.