01_root_finding.py

  1"""
  2근 찾기 (Root Finding)
  3Numerical Root Finding Methods
  4
  5f(x) = 0 을 만족하는 x를 수치적으로 찾는 방법들입니다.
  6"""
  7
  8import numpy as np
  9import matplotlib.pyplot as plt
 10from typing import Callable, Tuple, Optional
 11
 12
 13# =============================================================================
 14# 1. 이분법 (Bisection Method)
 15# =============================================================================
 16def bisection(
 17    f: Callable[[float], float],
 18    a: float,
 19    b: float,
 20    tol: float = 1e-10,
 21    max_iter: int = 100
 22) -> Tuple[float, int, list]:
 23    """
 24    이분법으로 f(x) = 0의 근 찾기
 25
 26    조건: f(a)와 f(b)의 부호가 달라야 함 (중간값 정리)
 27    수렴 속도: 선형 (매 반복마다 구간 절반)
 28
 29    Args:
 30        f: 목표 함수
 31        a, b: 초기 구간
 32        tol: 허용 오차
 33        max_iter: 최대 반복 횟수
 34
 35    Returns:
 36        (근, 반복 횟수, 중간값 히스토리)
 37    """
 38    if f(a) * f(b) >= 0:
 39        raise ValueError("f(a)와 f(b)의 부호가 달라야 합니다")
 40
 41    history = []
 42
 43    for i in range(max_iter):
 44        c = (a + b) / 2
 45        history.append(c)
 46
 47        if abs(f(c)) < tol or (b - a) / 2 < tol:
 48            return c, i + 1, history
 49
 50        if f(a) * f(c) < 0:
 51            b = c
 52        else:
 53            a = c
 54
 55    return (a + b) / 2, max_iter, history
 56
 57
 58# =============================================================================
 59# 2. 뉴턴-랩슨 방법 (Newton-Raphson Method)
 60# =============================================================================
 61def newton_raphson(
 62    f: Callable[[float], float],
 63    df: Callable[[float], float],
 64    x0: float,
 65    tol: float = 1e-10,
 66    max_iter: int = 100
 67) -> Tuple[float, int, list]:
 68    """
 69    뉴턴-랩슨 방법으로 f(x) = 0의 근 찾기
 70
 71    x_{n+1} = x_n - f(x_n) / f'(x_n)
 72
 73    수렴 속도: 2차 (제곱 수렴)
 74    단점: 도함수 필요, 초기값에 민감
 75
 76    Args:
 77        f: 목표 함수
 78        df: 도함수
 79        x0: 초기값
 80        tol: 허용 오차
 81        max_iter: 최대 반복 횟수
 82
 83    Returns:
 84        (근, 반복 횟수, 히스토리)
 85    """
 86    x = x0
 87    history = [x]
 88
 89    for i in range(max_iter):
 90        fx = f(x)
 91        dfx = df(x)
 92
 93        if abs(dfx) < 1e-15:
 94            raise ValueError("도함수가 0에 가까움: 발산 위험")
 95
 96        x_new = x - fx / dfx
 97        history.append(x_new)
 98
 99        if abs(x_new - x) < tol:
100            return x_new, i + 1, history
101
102        x = x_new
103
104    return x, max_iter, history
105
106
107# =============================================================================
108# 3. 할선법 (Secant Method)
109# =============================================================================
110def secant(
111    f: Callable[[float], float],
112    x0: float,
113    x1: float,
114    tol: float = 1e-10,
115    max_iter: int = 100
116) -> Tuple[float, int, list]:
117    """
118    할선법으로 f(x) = 0의 근 찾기
119
120    뉴턴법의 도함수를 차분으로 근사:
121    x_{n+1} = x_n - f(x_n) * (x_n - x_{n-1}) / (f(x_n) - f(x_{n-1}))
122
123    수렴 속도: 약 1.618차 (황금비)
124    장점: 도함수 불필요
125
126    Args:
127        f: 목표 함수
128        x0, x1: 두 초기값
129        tol: 허용 오차
130        max_iter: 최대 반복 횟수
131
132    Returns:
133        (근, 반복 횟수, 히스토리)
134    """
135    history = [x0, x1]
136
137    for i in range(max_iter):
138        f0, f1 = f(x0), f(x1)
139
140        if abs(f1 - f0) < 1e-15:
141            raise ValueError("분모가 0에 가까움")
142
143        x2 = x1 - f1 * (x1 - x0) / (f1 - f0)
144        history.append(x2)
145
146        if abs(x2 - x1) < tol:
147            return x2, i + 1, history
148
149        x0, x1 = x1, x2
150
151    return x1, max_iter, history
152
153
154# =============================================================================
155# 4. 고정점 반복법 (Fixed-Point Iteration)
156# =============================================================================
157def fixed_point(
158    g: Callable[[float], float],
159    x0: float,
160    tol: float = 1e-10,
161    max_iter: int = 100
162) -> Tuple[float, int, list]:
163    """
164    고정점 반복: x = g(x)를 만족하는 x 찾기
165
166    f(x) = 0 을 x = g(x) 형태로 변환
167    예: x² - 2 = 0  →  x = 2/x 또는 x = (x + 2/x)/2
168
169    수렴 조건: |g'(x*)| < 1
170
171    Args:
172        g: 반복 함수
173        x0: 초기값
174        tol: 허용 오차
175        max_iter: 최대 반복 횟수
176
177    Returns:
178        (고정점, 반복 횟수, 히스토리)
179    """
180    x = x0
181    history = [x]
182
183    for i in range(max_iter):
184        x_new = g(x)
185        history.append(x_new)
186
187        if abs(x_new - x) < tol:
188            return x_new, i + 1, history
189
190        x = x_new
191
192    return x, max_iter, history
193
194
195# =============================================================================
196# 5. Brent's Method (scipy와 비교)
197# =============================================================================
198def brents_method(
199    f: Callable[[float], float],
200    a: float,
201    b: float,
202    tol: float = 1e-10,
203    max_iter: int = 100
204) -> Tuple[float, int]:
205    """
206    Brent's Method (간소화 버전)
207    이분법, 할선법, 역2차보간을 조합한 방법
208
209    실무에서는 scipy.optimize.brentq 사용 권장
210    """
211    fa, fb = f(a), f(b)
212    if fa * fb >= 0:
213        raise ValueError("f(a)와 f(b)의 부호가 달라야 합니다")
214
215    if abs(fa) < abs(fb):
216        a, b = b, a
217        fa, fb = fb, fa
218
219    c, fc = a, fa
220    d = c  # d 초기화 (첫 iteration에서 사용됨)
221    mflag = True
222
223    for i in range(max_iter):
224        if fa != fc and fb != fc:
225            # 역2차 보간
226            s = (a * fb * fc / ((fa - fb) * (fa - fc)) +
227                 b * fa * fc / ((fb - fa) * (fb - fc)) +
228                 c * fa * fb / ((fc - fa) * (fc - fb)))
229        else:
230            # 할선법
231            s = b - fb * (b - a) / (fb - fa)
232
233        # 조건 체크 후 이분법으로 대체
234        conditions = [
235            not ((3 * a + b) / 4 <= s <= b or b <= s <= (3 * a + b) / 4),
236            mflag and abs(s - b) >= abs(b - c) / 2,
237            not mflag and abs(s - b) >= abs(c - d) / 2,
238        ]
239
240        if any(conditions):
241            s = (a + b) / 2
242            mflag = True
243        else:
244            mflag = False
245
246        fs = f(s)
247        d, c, fc = c, b, fb
248
249        if fa * fs < 0:
250            b, fb = s, fs
251        else:
252            a, fa = s, fs
253
254        if abs(fa) < abs(fb):
255            a, b = b, a
256            fa, fb = fb, fa
257
258        if abs(b - a) < tol or abs(fb) < tol:
259            return b, i + 1
260
261    return b, max_iter
262
263
264# =============================================================================
265# 시각화
266# =============================================================================
267def plot_convergence(f, methods_data, x_range, title="근 찾기 수렴 비교"):
268    """수렴 과정 시각화"""
269    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
270
271    # 왼쪽: 함수와 근
272    x = np.linspace(x_range[0], x_range[1], 500)
273    y = [f(xi) for xi in x]
274    axes[0].plot(x, y, 'b-', label='f(x)')
275    axes[0].axhline(y=0, color='k', linestyle='-', linewidth=0.5)
276
277    colors = plt.cm.tab10.colors
278    for i, (name, root, _, history) in enumerate(methods_data):
279        axes[0].scatter([root], [0], s=100, color=colors[i], zorder=5, label=f'{name}: x={root:.6f}')
280
281    axes[0].set_xlabel('x')
282    axes[0].set_ylabel('f(x)')
283    axes[0].set_title('함수와 근')
284    axes[0].legend()
285    axes[0].grid(True, alpha=0.3)
286
287    # 오른쪽: 수렴 속도
288    for i, (name, root, iters, history) in enumerate(methods_data):
289        if history:
290            errors = [abs(h - root) for h in history]
291            errors = [e if e > 1e-16 else 1e-16 for e in errors]
292            axes[1].semilogy(errors, 'o-', color=colors[i], label=f'{name} ({iters}회)')
293
294    axes[1].set_xlabel('반복 횟수')
295    axes[1].set_ylabel('오차 (log scale)')
296    axes[1].set_title('수렴 속도 비교')
297    axes[1].legend()
298    axes[1].grid(True, alpha=0.3)
299
300    plt.suptitle(title)
301    plt.tight_layout()
302    plt.savefig('/opt/projects/01_Personal/03_Study/Numerical_Simulation/examples/root_finding.png', dpi=150)
303    plt.close()
304    print("    그래프 저장: root_finding.png")
305
306
307# =============================================================================
308# 테스트
309# =============================================================================
310def main():
311    print("=" * 60)
312    print("근 찾기 (Root Finding) 예제")
313    print("=" * 60)
314
315    # 예제 1: f(x) = x³ - x - 2 = 0 (근 ≈ 1.5214)
316    print("\n[예제 1] f(x) = x³ - x - 2 = 0")
317    print("-" * 40)
318
319    f = lambda x: x**3 - x - 2
320    df = lambda x: 3*x**2 - 1
321
322    methods_data = []
323
324    # 이분법
325    root, iters, hist = bisection(f, 1, 2)
326    methods_data.append(("Bisection", root, iters, hist))
327    print(f"이분법:     근 = {root:.10f}, 반복 = {iters}")
328
329    # 뉴턴-랩슨
330    root, iters, hist = newton_raphson(f, df, 1.5)
331    methods_data.append(("Newton", root, iters, hist))
332    print(f"뉴턴-랩슨: 근 = {root:.10f}, 반복 = {iters}")
333
334    # 할선법
335    root, iters, hist = secant(f, 1, 2)
336    methods_data.append(("Secant", root, iters, hist))
337    print(f"할선법:     근 = {root:.10f}, 반복 = {iters}")
338
339    # 시각화
340    try:
341        plot_convergence(f, methods_data, (0, 3), "f(x) = x³ - x - 2")
342    except Exception as e:
343        print(f"    그래프 생성 실패: {e}")
344
345    # 예제 2: √2 구하기 (x² - 2 = 0)
346    print("\n[예제 2] √2 구하기 (x² - 2 = 0)")
347    print("-" * 40)
348
349    f2 = lambda x: x**2 - 2
350    df2 = lambda x: 2*x
351    g = lambda x: (x + 2/x) / 2  # 바빌로니아 방법
352
353    root, iters, _ = newton_raphson(f2, df2, 1.0)
354    print(f"뉴턴-랩슨:   √2 = {root:.15f}, 반복 = {iters}")
355
356    root, iters, _ = fixed_point(g, 1.0)
357    print(f"고정점 반복: √2 = {root:.15f}, 반복 = {iters}")
358
359    print(f"실제 √2:        {np.sqrt(2):.15f}")
360
361    # 예제 3: cos(x) = x 고정점
362    print("\n[예제 3] cos(x) = x (Dottie Number)")
363    print("-" * 40)
364
365    g_cos = lambda x: np.cos(x)
366    root, iters, _ = fixed_point(g_cos, 0.5)
367    print(f"고정점 x = cos(x): {root:.10f}, 반복 = {iters}")
368
369    print("\n" + "=" * 60)
370    print("근 찾기 방법 비교")
371    print("=" * 60)
372    print("""
373    | 방법        | 수렴 속도 | 장점                | 단점                |
374    |------------|----------|---------------------|---------------------|
375    | 이분법      | 선형     | 항상 수렴, 안정적    | 느림, 구간 필요      |
376    | 뉴턴-랩슨   | 2차      | 매우 빠름           | 도함수 필요, 발산 가능|
377    | 할선법      | ~1.618차 | 도함수 불필요        | 뉴턴보다 느림        |
378    | 고정점반복  | 선형~2차 | 간단                | 수렴 조건 확인 필요   |
379    | Brent      | 조합     | 안정 + 빠름         | 구현 복잡           |
380
381    실무 권장:
382    - scipy.optimize.brentq: 안정적인 근 찾기
383    - scipy.optimize.newton: 뉴턴-랩슨/할선법
384    - scipy.optimize.fsolve: 다변수 방정식
385    """)
386
387
388if __name__ == "__main__":
389    main()