x3−2x−5=0에는 근의 공식이 있지만 쓰기 번거롭고, x=cosx나 케플러 방정식 M=E−esinE에는 아예 공식이 없습니다. 공학의 방정식 대부분이 이런 꼴입니다. 이 장은 f(x)=0을 반복으로 푸는 법을 다룹니다. 구간을 반씩 줄이는 가장 안전한 방법(이분법)에서 출발해, 방정식을 x=g(x)로 고쳐 쓰는 고정점 반복, 접선으로 한 번에 다가가는 뉴턴 방법, 도함수 없이 비슷한 속도를 내는 할선법, 그리고 안전과 속도를 함께 얻는 브렌트 방법까지 갑니다. 모든 방법을 같은 잣대, 곧 오차가 한 단계마다 어떤 법칙으로 줄어드는가(수렴 차수)로 비교하고, 언제 멈춰야 하는지, 어떻게 실패하는지도 직접 봅니다.
이분법의 오차 한계로 필요한 반복 수를 미리 계산하고, 부호 판정만으로 근을 가두는 원리를 설명한다.
고정점 반복이 수렴할 조건 ∣g′∣<1을 축소 사상 정리로 이해하고, 거미줄 그림으로 수렴·발산을 읽는다.
수렴 차수와 점근 상수를 정의하고, 반복 표의 오차에서 차수를 실험으로 측정한다.
뉴턴 방법을 유도하고 2차 수렴을 증명하며, 중근에서 1차로 떨어지는 이유와 처방을 안다.
할선법의 차수가 황금비 1.618임을 보이고, 가위치법·일리노이 수정·브렌트 방법의 장단점을 비교한다.
잔차와 간격에 근거한 멈춤 조건을 고르고, 뉴턴 방법의 실패(순환, 발산, 평평한 곳)를 알아보고 피한다.
이 수식이 없었다면제곱근 버튼과 행성의 위치
계산기의 √ 버튼, 나눗셈 회로가 없는 프로세서의 나눗셈, 위성의 궤도 예보, 배관의 압력 손실 계산, 회로 시뮬레이터가 다이오드의 동작점을 찾는 일은 모두 같은 일을 합니다. 공식으로 풀 수 없는 방정식 f(x)=0을 어림값에서 시작해 고쳐 나가는 일입니다. 오늘날 수치 계산 라이브러리의 근 찾기 함수(MATLAB의 fzero, SciPy의 brentq)는 이 장의 이분법·할선법·역2차 보간을 조합한 방법입니다.
고쳐 나가는 생각은 오래되었습니다. 바빌로니아 점토판 YBC 7289에는 2가 60진법으로 1;24,51,10(십진법 1.41421296…)으로 적혀 있습니다. 17세기에 뉴턴(Isaac Newton)과 랩슨(Joseph Raphson)이 다항식 방정식을 푸는 반복법을 내놓았고, 같은 세기 초 케플러(Johannes Kepler)가 행성의 위치를 구하려고 세운 방정식 M=E−esinE는 지금도 궤도 계산 프로그램 안에서 매번 뉴턴 방법으로 풀립니다.
f(x)=x3−2x−5는 f(2)=−1, f(3)=16입니다. 연속함수가 음수에서 양수로 바뀌었으니 2와 3 사이 어딘가에 근이 있습니다. 함숫값의 부호만 물어볼 수 있다면, 근을 소수 열째 자리까지 알려면 몇 번 물어봐야 할까요? 그리고 그 횟수를 계산을 시작하기 전에 알 수 있을까요?
답의 뿌리는 중간값 정리입니다(CalculusBook 2장: 이분법). f가 [a,b]에서 연속이고 f(a)와 f(b)의 부호가 다르면 그 사이에 근이 적어도 하나 있습니다. 이 상태를 근을 가뒀다고 합니다.
정의 3.1 (근 가둠과 이분법) 연속함수 f에 대해 f(a)f(b)<0인 구간 [a,b]를 근의 가둠 구간(bracket)이라 한다. 이분법(bisection method)은 중점 m=(a+b)/2에서 f(m)의 부호를 보고, 부호가 바뀌는 쪽 절반 [a,m] 또는 [m,b]를 새 가둠 구간으로 삼는 일을 되풀이한다.
그림 3-1. 이분법은 함수의 모양을 전혀 쓰지 않는다. 중점의 부호 하나로 근이 있는 절반을 고르므로, 가둠 구간의 길이는 매번 정확히 절반이 된다. 근이 어디 있든 똑같은 속도다.
function bisect(f, a, b, tol)
fa = f(a) // sign(f(a)) != sign(f(b)) 라고 가정while (b - a) / 2 > tol
m = a + (b - a) / 2 // (a+b)/2 보다 넘침·반올림에 안전
fm = f(m)
if fm == 0 then return m
if sign(fm) == sign(fa) then a = m; fa = fm
else b = m
return a + (b - a) / 2
부호 판정을 fa*fm < 0으로 쓰면 두 값이 아주 작을 때 곱이 언더플로로 0이 되어 판정이 틀릴 수 있으므로 sign을 비교합니다. 반복마다 f를 한 번만 계산한다는 것도 기억해 두세요. 비용을 셀 때의 단위입니다.
정리 3.2 (이분법의 오차 한계)f가 [a0,b0]에서 연속이고 f(a0)f(b0)<0이면, n번 나눈 뒤의 가둠 구간 [an,bn]에는 근 r이 있고 그 중점 mn은
∣mn−r∣≤2bn−an=2n+1b0−a0
을 만족한다. 따라서 ∣mn−r∣≤ε을 보장하려면 n≥log2εb0−a0−1이면 된다.
증명. 매 단계 bk+1−ak+1=(bk−ak)/2이고, 부호가 바뀌는 절반을 남기므로 중간값 정리에 의해 근이 늘 구간 안에 있다. 구간 안의 점과 중점의 거리는 반 길이를 넘지 않는다. ∎
정리 3.2는 함수와 무관한 한계라는 점이 특별합니다. 근이 어디 있든, 함수가 얼마나 휘었든 같은 식입니다. 구간 [2, 3]에서 오차 10−10을 보장하려면 log21010−1=32.2, 곧 33번이면 되고, 근 근처에서 자릿수 하나(10배)를 얻는 데 평균 log210≈3.32번이 듭니다. 배정밀도에서는 길이 1인 구간을 약 52번 나누면 더 나눌 수 없을 만큼 좁아집니다(1장의 기계 엡실론).
f(2.5)=5.625>0이므로 근은 [2,2.5]에 있다. 같은 방식으로 f(2.25)>0, f(2.125)>0, f(2.0625)<0이 차례로 구간을 고른다.
오차가 단조롭게 줄지는 않는다. 5단계 중점 2.09375는 우연히 근에 매우 가까웠지만(8.0×10−4), 6단계 중점은 다시 멀어졌다(1.5×10−2). 줄어드는 것은 오차의 상한2−(n+1)이다.
f(2.09375)=−0.0089로 아주 작았는데도 이분법은 이 정보를 쓰지 못한다. 부호만 보기 때문이다. 함숫값의 크기까지 쓰는 방법이 뒤에 나올 할선법·가위치법이다.
SIMULATOR
이분법 한 단계씩: 가둠 구간이 반으로
목표 오차 ε
f(x)가둠 구간 [a, b]중점 m참 근
단계 n—
구간 길이 b − a—
실제 오차 |m − r|—
ε까지 필요한 n (정리 3.2)—
해볼 것: ① '다음 단계'를 누르며 구간 길이가 정확히 반씩 줄고, 실제 오차는 들쭉날쭉하지만 언제나 길이의 반 이하인 것을 확인하세요. ② 축 위의 두 끝점 a, b를 끌어 구간을 바꾸면 필요한 반복 수가 log₂(길이/ε)로 바뀝니다. 부호가 같은 구간은 거부됩니다. ③ tan x를 고르면 이분법은 근이 아니라 극점 π/2로 수렴합니다. 부호가 바뀐 이유가 근이 아니라 불연속이기 때문입니다. 모델: 중점 m = a + (b − a)/2, 부호만 판정.
고정점 반복과 축소 사상
같은 방정식, 다른 운명
x3−2x−5=0은 x=32x+5로도, x=(x3−5)/2로도, x=5/(x2−2)로도 고쳐 쓸 수 있습니다. 셋 다 오른쪽에 어림값을 넣고 나온 값을 다시 넣는 일을 되풀이할 수 있습니다. 같은 근을 가진 식인데 어떤 것은 근으로 다가가고 어떤 것은 무한대로 달아납니다. 무엇이 운명을 가를까요?
정의 3.3 (고정점과 고정점 반복)g(r)=r인 점 r을 g의 고정점(fixed point)이라 한다. 어림값 x0에서 xn+1=g(xn)을 되풀이하는 방법을 고정점 반복(fixed-point iteration)이라 한다. f(x)=0을 x=g(x) 꼴로 고쳐 쓰면 f의 근이 g의 고정점이 된다.
그림으로 보면 고정점은 곡선 y=g(x)와 직선 y=x의 교점입니다. 반복 한 번은 "곡선까지 세로로 가서(g(xn)을 계산), 대각선까지 가로로 가는(그 값을 새 x로)" 일이라서, 궤적이 거미줄이나 계단 모양을 그립니다. 이것을 거미줄 그림(cobweb diagram)이라 합니다.
그림 3-2. 거미줄 그림. 곡선이 대각선보다 완만하면(∣g′∣<1) 계단이 고정점으로 빨려 들어가고, 가파르면(∣g′∣>1) 밀려난다. 기울기가 음수면 계단 대신 고정점 둘레를 도는 거미줄이 된다.
고정점 근처에서 평균값 정리를 쓰면 xn+1−r=g(xn)−g(r)=g′(ξn)(xn−r)입니다. 한 번 반복할 때마다 오차에 g′가 곱해지는 셈입니다. 그래서 ∣g′∣<1이면 오차가 줄고, ∣g′∣>1이면 커집니다.
정리 3.4 (축소 사상 정리, 1차원)g가 [a,b]를 자기 자신으로 보내고(g([a,b])⊆[a,b]), 어떤 L<1에 대해 [a,b]에서 ∣g′(x)∣≤L이면
g는 [a,b]에 고정점 r을 정확히 하나 가진다.
임의의 x0∈[a,b]에서 xn+1=g(xn)은 r로 수렴하고 ∣xn−r∣≤Ln∣x0−r∣이다.
이런 g를 축소 사상(contraction mapping)이라 합니다. 이 정리는 바나흐(Stefan Banach)의 고정점 정리의 1차원 판이고, 4장의 연립방정식, 미분방정식 해의 존재 증명(피카르 반복), 13장의 반복 해법이 모두 같은 뼈대를 씁니다. 사후 추정 3번은 참 근을 모를 때 멈춤 판정에 쓸 수 있다는 점에서 실용적입니다. L=0.5이면 ∣xn−r∣≤∣xn−xn−1∣이지만, L=0.99이면 이웃 반복의 차이보다 오차가 99배까지 클 수 있습니다.
예제 3.2 (세 가지 고쳐 쓰기)x3−2x−5=0의 근 r=2.0945514815를 x0=2에서 세 가지 g로 반복하라.
g1(x)=(2x+5)1/3: g1′(r)=3(2r+5)2/32=0.1520. 반복값 2.0800838, 2.0923507, 2.0942170, 2.0945007, … 오차가 9.46×10−2→1.45×10−2→2.20×10−3→3.35×10−4으로, 매번 정확히 약 0.152배가 된다. 11번이면 10−10 아래다.
같은 근인데 ∣g′(r)∣가 1보다 작은 g1만 수렴한다. 수렴하는 경우에도 속도는 ∣g′(r)∣가 정한다. x=cosx는 ∣g′(r)∣=∣sin0.7390851∣=0.674라서 x0=1에서 10−10까지 55번이 걸린다.
SIMULATOR
거미줄 그림: 고정점 반복의 수렴과 발산
y = g(x)와 y = x, 반복 궤적
|xₙ − r| (로그 눈금)
x₀ (끌기)—
g′(r)—
측정한 |eₙ₊₁/eₙ|—
xₙ—
해볼 것: ① (2x + 5)^(1/3)에서 오른쪽 그래프가 직선(기울기 log 0.152)으로 내려가는 것을 보세요. 선형 수렴입니다. ② (x³ − 5)/2와 5/(x² − 2)는 x₀를 근 가까이 두어도 달아납니다. ③ 완화 반복 x − λ(x² − 2)에서 λ를 움직이면 g′(√2) = 1 − 2√2·λ가 바뀝니다. λ ≈ 0.354(= 1/(2√2))에서 g′ = 0이 되어 오차 그래프가 갑자기 꺾여 내려갑니다. 2차 수렴이고, 뉴턴 방법의 정체입니다. 모델: 축 위의 점 x₀를 끌어 출발점을 바꿉니다.
수렴 차수: 오차가 줄어드는 법칙
이분법은 한 번에 오차 상한이 1/2배, g1 반복은 0.152배, cosx 반복은 0.674배가 되었습니다. 모두 "매번 일정한 비율"로 줄어드는 같은 종류입니다. 거미줄 시뮬레이터의 λ = 0.354처럼 비율 자체가 0으로 가는 더 빠른 종류도 있습니다. 이를 구별하는 언어가 수렴 차수입니다.
정의 3.5 (수렴 차수)xn→r이고 en=xn−r이라 하자. 상수 p≥1과 C>0이 있어
n→∞lim∣en∣p∣en+1∣=C
이면 수열이 수렴 차수(order of convergence)p, 점근 상수 C로 수렴한다고 한다. p=1(이때 C<1이어야 한다)이면 선형 수렴(linear convergence), p=2이면 2차 수렴(quadratic convergence), ∣en+1∣/∣en∣→0이면 초선형 수렴(superlinear convergence)이라 한다.
차이를 숫자로 느껴 봅시다. 오차 10−1에서 출발해 C=0.5로 선형 수렴하면 10−16까지 약 50번이 걸립니다. 2차 수렴(C=1)이면 10−2,10−4,10−8,10−16, 곧 네 번입니다. 2차 수렴은 맞는 자릿수가 매번 두 배가 된다는 뜻입니다.
차수를 실험으로 재려면 ∣en+1∣≈C∣en∣p의 로그를 이웃한 두 단계에서 빼서 C를 지웁니다.
p≈log(∣en∣/∣en−1∣)log(∣en+1∣/∣en∣)
참 근을 모르면 en 대신 이웃 반복의 차이 xn+1−xn을 쓴다. 초선형 수렴에서는 xn+1−xn≈−en이기 때문이다. 오차가 반올림 수준(10−16 근처)에 닿은 단계는 빼야 한다.
또 하나의 그림이 있습니다. log∣en∣을 n에 대해 그리면 선형 수렴은 직선(기울기 logC), 2차 수렴은 아래로 점점 가팔라지는 곡선입니다. log∣en+1∣을 log∣en∣에 대해 그리면 어떤 방법이든 직선이 되고 그 기울기가 바로 p입니다. 이 장의 비교 시뮬레이터가 두 그림을 모두 보여 줍니다.
고정점 반복의 차수는 g의 도함수가 근에서 몇 번째까지 0인가로 정해집니다. 이것이 뉴턴 방법이 빠른 이유를 한 줄로 설명합니다.
정리 3.6 (고정점 반복의 차수)g가 r 근처에서 p번 연속 미분 가능하고 g′(r)=⋯=g(p−1)(r)=0, g(p)(r)=0이면, r에 충분히 가까운 x0에서 고정점 반복은 차수 p로 수렴하고 C=∣g(p)(r)∣/p!이다. p=1이면 0<∣g′(r)∣<1이어야 하고 C=∣g′(r)∣이다.
이분법은 부호만, 고정점 반복은 g만 씁니다. 현재 점에서 함숫값 f(xn)과 기울기 f′(xn)을 안다면 함수를 직선(접선)으로 어림하고 그 직선의 근으로 단숨에 건너갈 수 있습니다. 이것이 얼마나 빠르고, 언제 믿을 수 있을까요?
CalculusBook에서 접선 근사로 뉴턴 방법을 만났습니다(CalculusBook 6장: 뉴턴 방법). 여기서는 오차 분석을 끝까지 합니다. xn에서의 테일러 전개를 1차에서 자르면 f(x)≈f(xn)+f′(xn)(x−xn)이고, 이 직선이 0이 되는 점을 다음 어림값으로 삼습니다.
정의 3.7 (뉴턴 방법)f′(xn)=0일 때
xn+1=xn−f′(xn)f(xn),n=0,1,2,…
로 반복하는 방법을 뉴턴 방법(Newton's method, Newton–Raphson)이라 한다. 한 번에 f와 f′을 하나씩 계산한다. 고정점 반복으로 보면 g(x)=x−f(x)/f′(x)이다.
그림 3-3. 뉴턴 방법의 두 단계. 현재 점의 접선이 x축과 만나는 곳이 다음 어림값이다. 곡선과 접선의 차이는 21f′′(x−xn)2 정도라서, 근에 가까워질수록 접선의 근이 참 근에 제곱으로 가까워진다.
function newton(f, df, x, tol, maxit)
for n = 1 .. maxit
d = df(x)
if d == 0 then fail("평평한 곳: 접선이 x축과 만나지 않음")
s = f(x) / d
x = x - s
if |s| <= tol * (1 + |x|) then return x // 걸음 크기 판정 (멈춤 조건 절 참고)
fail("maxit 안에 수렴하지 않음")
정리 3.8 (뉴턴 방법의 2차 수렴)f가 근 r 근처에서 두 번 연속 미분 가능하고 f(r)=0, f′(r)=0이면, r에 충분히 가까운 x0에서 뉴턴 방법은 r로 수렴하고
증명의 핵심.xn에서 테일러 전개해 r을 넣으면 0=f(r)=f(xn)+f′(xn)(r−xn)+21f′′(ξn)(r−xn)2. 양변을 f′(xn)으로 나누고 뉴턴 식 xn+1=xn−f(xn)/f′(xn)을 쓰면 xn+1−r=2f′(xn)f′′(ξn)(xn−r)2. 근 근처에서 ∣f′′/(2f′)∣≤K이고 K∣e0∣<1이면 ∣en+1∣≤(K∣en∣)∣en∣이라 오차가 줄고, 줄수록 더 빨리 준다. 정리 3.6으로 보면 g′(x)=f(x)f′′(x)/f′(x)2이라 g′(r)=0이다. ∎
정리의 "충분히 가까운"이 중요합니다. 뉴턴 방법은 국소적으로만 수렴이 보장됩니다. 상수 K=∣f′′/(2f′)∣가 크면(근에서 기울기가 작거나 곡률이 크면) 2차 수렴이 시작되는 반경 1/K이 좁습니다. 그 바깥에서는 뒤의 실패 갤러리처럼 무슨 일이든 일어날 수 있습니다.
예제 3.3 (뉴턴 방법 손 계산)f(x)=x3−2x−5, f′(x)=3x2−2, x0=2.
n
xₙ
f(xₙ)
|eₙ|
|eₙ|/|eₙ₋₁|²
0
2
−1
9.46×10⁻²
1
2.1
6.10×10⁻²
5.45×10⁻³
0.610
2
2.0945681211041851
1.86×10⁻⁴
1.66×10⁻⁵
0.560
3
2.0945514816981992
1.74×10⁻⁹
1.56×10⁻¹⁰
0.563
4
2.0945514815423265
−8.9×10⁻¹⁶
0
x1=2−(−1)/10=2.1. 첫 단계에서 이미 오차가 17분의 1이다.
맞는 자릿수: x1 2자리, x2 5자리, x3 10자리, x4는 배정밀도 전부. 매번 거의 두 배다.
보통 뉴턴: 오차 1→0.556→0.298→0.155→0.0796→0.0403→0.0203. 비율이 (m−1)/m=1/2로 간다. 10−6까지 21번이 걸리고, 27번째에 9.7×10−9에서 더 나아가지 못한다.
수정 뉴턴(m=2): 오차 1→0.111→1.95×10−3→6.33×10−7→3.83×10−11. 2차 수렴이 돌아왔지만 3.8×10−11에서 멈춘다.
둘 다 10−16까지 가지 못하는 이유: 이중근 근처에서 f(x)≈3(x−1)2이라 ∣x−1∣=10−8이면 f≈3×10−16, 이미 f를 계산하는 반올림 오차와 같은 크기다. m중근은 대략 ε1/m까지만 구할 수 있다. 근 자체의 조건수가 나쁜 문제다(2장: 조건수).
할선법: 도함수 없이
f′을 식으로 구할 수 없거나 계산이 비싸면(예: f가 시뮬레이션 한 번), 접선 대신 마지막 두 점을 잇는 할선을 씁니다. 도함수를 차분 몫 xn−xn−1f(xn)−f(xn−1)으로 바꾼 뉴턴 방법입니다.
정의 3.10 (할선법) 어림값 두 개 x0,x1에서 출발해
xn+1=xn−f(xn)f(xn)−f(xn−1)xn−xn−1
로 반복하는 방법을 할선법(secant method)이라 한다. 앞 단계의 f(xn−1)을 다시 쓰므로 한 번에 f를 한 번만 새로 계산한다.
정리 3.11 (할선법의 차수는 황금비)f′(r)=0, f′′(r)=0이고 x0,x1이 r에 충분히 가까우면
en+1≈2f′(r)f′′(r)enen−1,p=21+5≈1.618.
증명의 핵심. 할선은 xn−1,xn을 지나는 1차 보간이고, 보간 오차 공식(5장)으로 f(x)−ℓ(x)=21f′′(ξ)(x−xn)(x−xn−1). x=r을 넣고 정리하면 위의 곱 공식이 나온다. ∣en+1∣≈C∣en∣p라 놓고 ∣en−1∣≈(∣en∣/C)1/p를 넣으면 지수 비교에서 p=1+1/p, 곧 p2=p+1. 양의 근이 황금비다. ∎
차수 1.618은 뉴턴의 2보다 작지만 한 단계 비용이 절반입니다. 공정한 비교 잣대는 효율 지수(efficiency index)p1/w(w는 한 단계의 함수 계산 수)입니다. f′의 비용이 f와 같다면 뉴턴은 21/2=1.414, 할선법은 1.6181/1=1.618. 함수 계산 한 번당 할선법이 더 많이 줄입니다. 대신 할선법은 f(xn)−f(xn−1)에서 자릿수 상쇄가 생기고(1장), 근을 가두지 않으므로 뉴턴처럼 달아날 수 있습니다.
예제 3.5 (할선법) 같은 f=x3−2x−5에 x0=2, x1=3으로 할선법을 써라.
n
xₙ
|eₙ|
|eₙ|/(|eₙ₋₁||eₙ₋₂|)
2
2.0588235294117645
3.57×10⁻²
3
2.0812636598450229
1.33×10⁻²
4
2.0948241460940524
2.73×10⁻⁴
0.574
5
2.0945494310352473
2.05×10⁻⁶
0.566
6
2.0945514812275992
3.15×10⁻¹⁰
0.563
7
2.0945514815423270
4.4×10⁻¹⁶
비율 ∣en+1∣/(∣en∣∣en−1∣)이 정리 3.11의 0.563으로 모인다. 오차의 지수 −1.4,−1.9,−3.6,−5.7,−9.5가 앞의 두 지수의 합에 가깝게 자란다(피보나치 수열처럼). f를 2 + 6 = 8번 계산해 배정밀도에 닿았다. 뉴턴 방법은 같은 정확도에 f와 f′을 4번씩, 합해서 8번 썼다.
SIMULATOR
뉴턴 접선 · 할선 · 수정 뉴턴 한 단계씩
방법
f(x)접선 / 할선반복점참 근
단계 n—
xₙ—
|eₙ|—
측정 차수 p—
해볼 것: ① x³ − 2x − 5에서 뉴턴을 4번 누르며 표의 |eₙ|/|eₙ₋₁|²가 0.563으로 모이는 것을 보세요. ② 이중근 함수로 바꾸면 비율 |eₙ|/|eₙ₋₁|이 0.5에 머뭅니다(1차). '수정 (m = 2)'를 고르면 다시 2차입니다. 단순근 함수에서 수정 방법을 쓰면 근을 지나쳐 오히려 1차로 느려집니다. ③ 할선법은 끄는 점이 두 개(x₀, x₁)입니다. 측정 차수가 1.6 근처에 머뭅니다. 모델: 참 근은 미리 고정밀로 구해 둔 값, 측정 차수는 마지막 세 오차로 계산합니다.
가위치법과 일리노이 수정
할선법의 속도와 이분법의 안전을 함께
할선법은 빠르지만 근을 가두지 않아 달아날 수 있고, 이분법은 안전하지만 함숫값의 크기를 버립니다. 가둠 구간의 두 끝을 잇는 할선을 쓰면 둘을 함께 얻을 수 있지 않을까요? 이 자연스러운 생각에는 함정이 하나 있습니다.
가둠 구간 [a,b]의 두 끝점 (a,f(a)), (b,f(b))를 잇는 직선이 x축과 만나는 점
c=f(b)−f(a)af(b)−bf(a)를 중점 대신 쓰고, 부호가 바뀌는 쪽을 남기는 방법이 가위치법(regula falsi, method of false position)입니다. 근은 늘 갇혀 있으니 안전합니다.
그림 3-4. 가위치법의 함정. 근 근처에서 f가 볼록(또는 오목)하면 할선의 근이 늘 같은 쪽에 떨어져, 한 끝점이 처음 자리에 묶인다. 묶인 끝점이 멀리 있으니 할선의 기울기가 접선과 달라 매번 일정한 비율밖에 줄지 않는다.
함정은 그림 3-4입니다. 묶인 끝점 b 때문에 오차가 en+1≈Cen(선형)으로만 줄고, f(b)가 클수록 C가 1에 가까워집니다. 처방은 놀랍도록 간단합니다. 같은 끝점이 두 번 연속 남으면 그 끝점의 함숫값을 반으로 줄여 할선을 그쪽으로 눕힙니다. 이것이 일리노이 수정(Illinois algorithm)입니다.
function illinois(f, a, b, tol)
fa = f(a); fb = f(b); side = 0
while |b - a| > tol
c = (a*fb - b*fa) / (fb - fa) // 할선의 근 (가위치)
fc = f(c)
if sign(fc) == sign(fb) // c가 b를 대신함, a는 남음
b = c; fb = fc
if side == -1 then fa = fa / 2 // a가 두 번 연속 남았다
side = -1
else// c가 a를 대신함, b는 남음
a = c; fa = fc
if side == +1 then fb = fb / 2
side = +1
return c
예제 3.6 (묶인 끝점 풀어 주기)f(x)=x3−2x−5, [2,3].
가위치법: c = 2.0588235, 2.0812637, 2.0896392, 2.0927396, … 모두 근의 왼쪽이고 b=3은 끝까지 그대로다. 오차 비율이 0.3719,0.3697,0.3689,⋯→0.3684로 선형이다. 10−10까지 21번.
일리노이: 처음 두 번은 같다. 3번째에 f(3)=16을 8로 줄인 할선이 c=2.0978634로 근의 오른쪽에 떨어져 b가 드디어 움직인다. 이후 오차 2.49×10−5,4.63×10−8,4.62×10−8,8.9×10−16. 7번이면 끝난다.
이분법은 오차를 보장하려면 33번이 필요하다(이 예에서는 32번째 중점에서 처음 10−10 아래로 내려간다). 더 볼록한 함수에서 차이가 커진다. x10−1을 [0,1.3]에서 풀면 10−10까지 가위치법 95번, 일리노이 14번이고, ex−2를 [0,4]에서 풀면 171번 대 10번이다.
일리노이 방법의 수렴 차수는 약 31/3≈1.44로 알려져 있다. 반으로 줄이는 단계와 할선 단계가 섞여 있기 때문이다.
SIMULATOR
가위치법 vs 일리노이 vs 이분법: 묶인 끝점
선택한 방법의 할선과 가둠 구간
|오차| vs 반복 n (세 방법)
위 그림에 보일 방법
가위치법일리노이이분법
단계 n—
가위치 |e|—
일리노이 |e|—
이분 |e|—
해볼 것: ① 가위치법을 보며 단계를 넘기면 한쪽 끝점(빨간 고리)이 움직이지 않습니다. 오른쪽 그래프에서 빨간 선은 기울기가 일정한 직선(선형 수렴)입니다. ② 일리노이로 바꾸면 반으로 줄인 함숫값(빈 원)이 할선을 눕혀 묶인 끝점을 풀어 줍니다. ③ x¹⁰ − 1은 가위치법이 이분법보다도 훨씬 느린 예입니다. 모델: 반복 최대 60번, 오차는 참 근과의 거리.
브렌트 방법: 안전하게, 그러나 빠르게
실무의 근 찾기 함수는 하나의 공식이 아니라 전략입니다. 브렌트 방법(Brent's method)은 데커(T. J. Dekker)의 아이디어를 브렌트(Richard Brent)가 다듬은 것으로, 늘 가둠 구간을 유지하면서 매 단계 가장 빠른 후보를 시도하고, 후보가 미덥지 않으면 이분법으로 물러납니다.
후보 1 – 역2차 보간: 최근 세 점 (a,fa),(b,fb),(c,fc)을 지나는 이차식을 x를 y의 함수로 보고 세운 뒤(역2차 보간(inverse quadratic interpolation)), y=0을 넣는다. 근을 구할 필요 없이 값 하나가 나온다. 수렴 차수 약 1.839.
후보 2 – 할선: 세 점이 없거나 두 값이 같을 때.
안전장치: 후보가 가둠 구간의 안쪽 3/4 바깥이거나, 걸음이 지지난 걸음의 절반보다 크게 줄지 않으면 이분법 한 걸음을 대신 둔다.
역2차 보간의 라그랑주 형식(5장). 세 y값이 서로 달라야 한다. 이 x가 다음 후보다.
function brent(f, a, b, tol) // f(a), f(b)의 부호가 다름loop
b = 지금까지 |f|가 가장 작은 점, c = b와 부호가 다른 끝점 (가둠 [b, c])
if |c - b| / 2 <= tol then return b
if 직전 걸음이 충분히 컸다
x = 역2차 보간(세 점이 있으면) or 할선
if x가 가둠 구간 안쪽이고 걸음이 직전 걸음의 절반보다 작다
then 받아들임 else x = 이분 중점
else x = 이분 중점
f(x) 계산, 가둠 구간 갱신
브렌트 방법은 매끄러운 함수에서 초선형으로 수렴하고, 최악의 경우에도 이분법의 몇 배를 넘지 않도록 설계되어 있습니다(브렌트는 이분법이 N번 걸리는 문제에서 많아야 약 N2번임을 보였다). 아래 비교 시뮬레이터의 구현으로 x3−2x−5를 [2,3]에서 풀면 f를 8번 계산하고 끝나며, 단계는 할선 → 역2차 → 할선 → 할선 → 역2차 → 할선 순입니다. x10−1에서는 처음 할선이 멀리 떨어지자 이분 단계가 두 번 끼어든 뒤 다시 빠른 단계로 돌아가 10번 만에 끝납니다.
방법 비교: 오차 그래프의 기울기
지금까지의 방법을 한 그림에 모아 봅시다. 세로축은 로그 눈금의 오차, 가로축은 반복 수 또는 함수 계산 횟수입니다. 반복 한 번의 값이 방법마다 다르므로(뉴턴은 f와 f′, 나머지는 f 하나) 공정한 비교는 후자입니다. 두 번째 보기 log∣en+1∣ 대 log∣en∣에서는 각 방법이 직선이 되고 그 기울기가 수렴 차수입니다.
SIMULATOR
여섯 방법의 오차 경주: 수렴 차수를 기울기로 읽기
그래프
이분고정점(기울기 고정)가위치일리노이할선뉴턴브렌트
10⁻¹² 도달 f 계산 수: 이분—
일리노이—
할선—
뉴턴 (f + f′)—
브렌트—
해볼 것: ① 'eₙ₊₁ vs eₙ'을 고르면 이분·가위치·고정점은 기울기 1, 할선은 약 1.6, 뉴턴은 2의 직선이 됩니다(기준선과 비교). ② 'f 계산 수'로 바꾸면 뉴턴과 할선의 차이가 줄거나 뒤집힙니다. 효율 지수 √2 = 1.41 대 1.62입니다. ③ x¹⁰ − 1에서 x₀ = 1.3의 뉴턴은 처음 몇 단계를 거의 일정한 비율(0.9배)로 기어가다가 근 근처에서 갑자기 2차로 바뀝니다. 브렌트 아래 메시지에서 단계 종류를 보세요. 모델: 고정점 반복은 기울기를 처음 값 f′(x₀)로 고정한 "평행선 방법"입니다.
방법
가둠
필요한 것
차수 p
단계당 계산
효율 지수
이분법
예
f의 부호
1 (C = 1/2)
f 1번
1 (확실)
고정점 반복
아니오
g, |g′(r)| < 1
1 (C = |g′(r)|)
g 1번
—
가위치법
예
f
1 (끝점 묶임)
f 1번
—
일리노이
예
f
≈ 1.44
f 1번
1.44
할선법
아니오
f, 시작점 2개
1.618
f 1번
1.618
뉴턴 방법
아니오
f, f′
2 (단순근)
f, f′ 1번씩
1.414
브렌트 방법
예
f
초선형(≈ 1.8)
f 1번
≈ 1.8, 최악에도 안전
규칙은 단순합니다. 가둘 수 있으면 브렌트(fzero, brentq), 도함수가 싸고 좋은 어림값이 있으면 뉴턴(케플러 방정식, 제곱근, 음적 ODE 방법의 매 걸음), 다변수라면 4장의 다변수 뉴턴과 그 변형입니다.
언제 멈출까: 잔차와 간격
|f(x)|가 10⁻¹⁴인데 오차는 0.01
p(x)=(x−1)7을 전개한 다항식 x7−7x6+21x5−⋯−1에서 p(1.01)=8.0×10−15입니다. "∣f(x)∣<10−12이면 멈춘다"는 프로그램은 1.01을 근이라고 보고하겠지만 오차는 0.01입니다. 반대로 f(x)=1010(x−2)에서는 오차가 10−12인 점도 ∣f∣=0.01이라 영원히 멈추지 않습니다. 무엇을 기준으로 멈춰야 할까요?
멈춤 조건에는 세 종류가 있습니다.
잔차(residual)∣f(xn)∣≤εf: f를 얼마나 0에 가깝게 만들었는가. 근 근처에서 f(x)≈f′(r)(x−r)이므로 오차 ≈ 잔차 / |f′(r)|. 기울기가 작으면(평평하면) 잔차가 작아도 오차가 크고, 기울기가 크면 그 반대다. 잔차는 2장의 후방 오차, 근의 오차는 전방 오차다.
걸음∣xn+1−xn∣≤εabs+εrel∣xn+1∣: 초선형 방법에서는 ∣xn+1−xn∣≈∣en∣이라 좋은 오차 추정이다. 절대·상대 허용 오차를 섞는 것은 근이 106이든 10−6이든 의미 있게 하려는 것이다. 선형 수렴에서는 정리 3.4의 L/(1−L) 배를 잊지 말 것.
가둠 구간 길이∣b−a∣≤2ε: 가두는 방법에서만, 하지만 가장 확실하다.
그리고 어떤 조건도 반올림 바닥 아래로는 내려갈 수 없습니다. f를 계산한 값에는 대략 δ≈ε⋅∑∣akxk∣ 크기의 잡음이 섞이고, ∣f(x)∣가 그보다 작아지는 근 근처의 구간에서는 계산한 f의 부호가 무작위가 됩니다. 이 불확실 구간의 반지름은 대략 δ/∣f′(r)∣(단순근), m중근이면 (δm!/∣f(m)(r)∣)1/m입니다. 그래서 허용 오차는 4ε∣x∣보다 작게 잡지 않고, 반복 수 상한(maxit)을 꼭 둡니다.
[0.98,1.02]에서 계산값과 참값 (x−1)7의 차이는 최대 6.0×10−15이다. 참값이 이보다 작아지는 ∣x−1∣<(6×10−15)1/7=0.0093에서는 부호를 믿을 수 없다. 실제로 [0.99,1.01]을 2000칸으로 나눠 보면 계산한 p의 부호가 694번 바뀐다.
이분법은 여전히 "부호가 바뀌는 점"을 찾아 주지만, 그 점은 근이 아니다. [0.5,1.6]에서 시작하면 0.99978(오차 2.2×10−4), [0.7,1.25]에서 시작하면 0.99648(오차 3.5×10−3)을 보고한다.
교훈: 문제의 조건수가 나쁘면(m중근, 평평한 함수) 어떤 알고리즘도 불확실 구간보다 정확할 수 없다. 이때 할 일은 알고리즘이 아니라 문제를 바꾸는 것, 예를 들어 인수 형태 (x−1)7을 그대로 계산하는 것이다.
SIMULATOR
근 근처의 반올림 잡음: 불확실 구간 확대경
참 곡선배정밀도 계산값잡음 크기 ±δ이분법이 보고한 근
잡음 크기 δ (창 안 최대)—
불확실 구간 반지름—
창 안 부호 바뀜—
이분법 결과 오차 범위—
해볼 것: ① (x − 1)⁷에서 창을 10⁻² 근처로 확대하면 매끄럽던 곡선이 톱니 같은 잡음으로 바뀝니다. 잡음 띠 안에서는 부호가 무작위입니다. ② '이분법 20번 실행'은 무작위 처음 구간 20개로 이분법을 돌려 보고된 근을 빨간 점으로 찍습니다. 점들이 불확실 구간 전체에 흩어집니다. ③ 단순근 (x − 1)(x − 2)(x − 3)에서는 10⁻¹⁵까지 확대해야 잡음이 보입니다. 불확실 구간이 δ/|f′(1)| = δ/2 정도로 작습니다. 모델: 계수로 호너 계산, 참값은 곱 형태.
뉴턴 방법이 실패할 때
뉴턴 방법의 수렴은 국소적입니다. 출발점이 나쁘면 다음과 같은 일이 실제로 일어납니다.
순환: f(x)=x3−2x+2에서 x0=0이면 x1=1, x2=0, … 영원히 0과 1을 오간다. 실근은 −1.7692923542인데 근처에도 가지 않는다. x0=0.1이면 1.0142, 0.0797, 1.0091, 0.0522, …로 순환 근처를 맴돌다가 천천히 벗어난다.
발산: f(x)=arctanx는 근 0이 하나뿐이고 단조인데도, ∣x0∣>1.3917452이면 ∣xn∣이 점점 커진다. x0=1.5이면 −1.694,2.321,−5.114,32.30,−1575. 경계값은 2x=(1+x2)arctanx의 근, 곧 2주기 궤도다. x0=1.39는 경계 안쪽이라 천천히 수렴하고, 1.4는 바깥이라 천천히 발산한다.
평평한 곳: f′(xn)≈0이면 접선이 거의 수평이라 아주 멀리 튄다. x3−2x−5에서 x0=0.82(극값 2/3=0.8165 근처)이면 x1=354.8. 그 뒤로는 매번 거리의 약 2/3씩만 돌아온다(멀리서 x3에 뉴턴을 쓰는 셈).
근에서 미분 불가능: f(x)=3x이면 xn+1=xn−3xn=−2xn. 어떤 x0=0에서도 매번 두 배로 멀어진다. 정리 3.8의 가정 f′′ 연속이 깨졌다.
엉뚱한 근: 근이 여럿이면 어느 근으로 갈지 출발점이 정하는데, 그 경계가 매우 복잡하다. 4장의 뉴턴 프랙털이 이 현상을 복소평면에서 보여 준다.
처방은 셋입니다. ① 가둠 구간을 알면 안전장치를 단 뉴턴: 뉴턴 걸음이 가둠 구간 밖으로 나가거나 충분히 줄지 않으면 이분 걸음을 대신 둔다(브렌트와 같은 생각, Numerical Recipes의 rtsafe). ② 걸음을 줄이는 감쇠 뉴턴: ∣f∣가 줄어들 때까지 걸음을 반으로 줄인다(4장). ③ 물리적 직관이나 근사식으로 좋은 출발점을 고른다(케플러 방정식의 E0=π, 콜브룩 방정식의 하알란 근사).
SIMULATOR
뉴턴 방법 실패 갤러리와 안전장치
사례
f(x)접선반복점참 근가둠 구간
x₀ (끌기)—
xₙ—
|f(xₙ)|—
이분 걸음 수—
해볼 것: ① '순환'에서 x₀를 0에서 조금씩 움직여 보세요. 대부분은 결국 순환을 벗어나 근 −1.769로 가지만, 0 바로 근처는 오래 갇힙니다. ② '발산'에서 x₀를 1.39와 1.40 사이로 옮기며 수렴과 발산의 경계 1.3917을 찾으세요. ③ 아무 사례에서나 '안전장치'를 켜면 걸음이 가둠 구간 밖으로 나갈 때마다 이분 걸음(초록 구간이 반으로)으로 바뀌어 언제나 근에 닿습니다. 모델: 가둠 구간은 사례마다 미리 정한 부호 바뀜 구간입니다.
응용: 제곱근, 역수, 케플러 방정식, 콜브룩 방정식
제곱근과 역수: 하드웨어 안의 뉴턴
f(x)=x2−a에 뉴턴 방법을 쓰면 xn+1=xn−2xnxn2−a=21(xn+xna)입니다. "어림값과 a/어림값의 평균"이라는 이 식을 흔히 바빌로니아 방법이라 부릅니다. 2를 x0=1에서 구하면 1.5, 1.4166667, 1.4142157, 1.41421356237469, 1.4142135623730949로 오차가 8.6×10−2,2.5×10−3,2.1×10−6,1.6×10−12,2.2×10−16입니다. 다섯 번이면 배정밀도 끝입니다.
곱셈과 뺄셈만 쓴다. 1−axn+1=(1−axn)2이라 상대 오차가 매번 정확히 제곱된다. 0<x0<2/a이면 수렴한다.
a=7, x0=0.1이면 1−axn이 0.3,0.09,0.0081,6.6×10−5,4.3×10−9이고 다섯 번째에 1/7이 배정밀도로 나옵니다. 곱셈기는 있지만 나눗셈 회로가 느리거나 없는 하드웨어(GPU, DSP 등)가 표에서 읽은 어림값에 이런 반복을 몇 번 더해 역수와 역제곱근 1/a(반복 x(3−ax2)/2)을 만듭니다.
케플러 방정식: 행성은 어디에 있는가
이심률 e인 타원 궤도에서 근일점 통과 후 시간 t의 위치는 평균 근점 이각 M=2πt/T(시간에 비례)로부터 이심 근점 이각 E를 구해 얻습니다. 둘을 잇는 식이 케플러 방정식(Kepler's equation)입니다.
M=E−esinE,x=a(cosE−e),y=a1−e2sinE
E를 알면 태양(초점)을 원점으로 한 위치가 바로 나온다. 0≤e<1이면 f(E)=E−esinE−M은 f′=1−ecosE>0이라 단조 증가, 근이 하나다.
예제 3.8 (케플러 방정식)M=1, e=0.5일 때 E를 구하라.
고정점 반복 E=M+esinE(케플러 자신이 쓴 방식과 같은 꼴): g′(E)=ecosE라 ∣g′∣≤e<1, 늘 수렴한다. E0=1에서 1.4207355, 1.4943810, 1.4985409, 1.4986954, …로 ∣ΔE∣<10−15까지 12번.
이심률이 1에 가까우면 사정이 달라진다. e=0.99, M=0.1에서 E0=M으로 뉴턴을 쓰면 f′(0.1)=1−0.99cos0.1=0.015로 평평해서 E1=6.71, E2=−55.3으로 튀고 21번 만에 겨우 E=0.8316604에 닿는다. E0=π에서 시작하면 3.1416, 1.6132, 1.1102, 0.8899, 0.8351, 0.8317로 8번이면 된다. 고정점 반복은 ∣g′∣≈e라 50번으로도 모자란다.
SIMULATOR
케플러 방정식: 궤도 위 행성의 위치를 뉴턴으로
궤도 (초점 = 태양), 보조원 위의 E
|f(Eₙ)| vs 반복 n
뉴턴 출발점 E₀
뉴턴고정점 E = M + e sin E보조원 · 균일 각 M
E—
뉴턴 반복 수—
고정점 반복 수—
태양 거리 r = a(1 − e cos E)—
해볼 것: ① ▶ 공전을 누르면 M이 일정하게 늘지만 행성은 근일점(오른쪽 끝) 근처에서 빠르고 원일점에서 느립니다. 케플러 제2법칙이고, 매 프레임 뉴턴 방법이 E를 풉니다. ② e = 0.99, M = 0.1로 두고 E₀ = M과 E₀ = π의 반복 수를 비교하세요. ③ 고정점 반복은 e가 클수록 느려집니다(점근 비율 e cos E). 모델: 장반경 a = 1, 멈춤 |f(E)| < 10⁻¹⁴ 또는 50번.
콜브룩 방정식: 배관의 마찰 계수
난류 배관 유동의 압력 손실은 다르시 마찰 계수 f로 계산하고, 표준 공학 관계식인 콜브룩(Colebrook) 방정식은 f에 대한 음함수입니다.
f1=−2log10(3.7ε/D+Ref2.51)
ε/D는 상대 거칠기, Re는 레이놀즈 수. y=1/f로 두면 y=G(y) 꼴의 고정점 문제가 된다.
예제 3.9 (콜브룩 방정식의 고정점 반복)Re=105, ε/D=10−4일 때 f를 구하라.