[Python]비선형 방정식 - 이분법, 뉴턴법, 할선법

비선형 방정식

• 2차 이상 방정식들을 의미(곡선)

이분법

이분법 Bisection Method

• 연속함수의 속성을 활용

• 각 단계에서 구간 [𝑎, 𝑏]와 𝑢 = 𝑓(𝑎) , 𝑣 = 𝑓(𝑏)의 값이 주어진다.

  • 이때 𝑢와 𝑣는 𝑢𝑣 < 0을 만족(부호가 다르다는 것)

• 이 폐구간의 중점인 \(𝑐 = \frac{1}{2}(a+b)\) 를 정하고, 𝑤 = 𝑓(𝑐)를 구한다.

• 우연히 𝑓(𝑐)=0이라면 알고리즘 종료 후 값 반환

• 보통은 𝑤 ≠ 0, 𝑤𝑢 < 0 || 𝑤𝑣 < 0

• if 𝑤𝑢 < 0 ?

  • 𝑓의 근이 구간 [𝑎, 𝑐]에 존재함

  • 구간을 축소하여 계속 진행

• If 𝑤𝑢 > 0 ?

• 𝑓의 근이 구간 𝑎, 𝑐 사이에 존재하는지의 여부를 알 수 없음.

• 그러나 𝑤𝑣 < 0이므로 [𝑐, 𝑏]에 𝑓의 근이 반드시 존재, 구간을 축소하여 계속 진행

• 구간의 크기가 충분히 줄어들 때까지 반복 (ex) \(|𝑏 − 𝑎| < \frac{1}{2} * 10^{-6}\)

  • 오차 범위 이내여야 한다고 보면 됨

• 이때 가장 합리적인 근의 추정치 : \(\frac{a+b}{2}\)

이분법의 수렴성 분석

• 𝑓가 \([𝑎_0, 𝑏_0]\)의 양 끝에서 서로 다른 부호의 함숫값을 갖는 연속함수.

• \([𝑎_0, 𝑏_0]\) 내에 한 근 𝑟이 존재

• \(c_0 = \frac{𝑎_0+𝑏_0}{2}\)를 근 r의 근삿값으로 잡으면 \(|r - c_0| ≤ \frac{b_0 - a_0}{2}\)

import numpy as np
import sympy as sp

def Bisection(f,a,b,e,nmax):
    fa = f.evalf(subs={x:a})
    fb = f.evalf(subs={x:b})
    if(fa*fb>0):
        print("f(a) and f(b) have same signs.")
        return
    err = b-a
    for i in range(1, nmax+1):
        err = err/2
        c = a + err
        fc = f.evalf(subs={x:c})
        if abs(err)<e:
            print("convergence at " +str(i) +", x = "+str((a+b)/2.0)+", f(x) = "+str(fc))
            return
        if (fa*fc<0):
            b=c
            fb=fc
        else:
            a=c
            fa=fc

뉴턴법

뉴턴법 Newton-Raphson Method

• 미분 가능한 함수 𝑓에 대해 적용

  • 함수 𝑓가 모든 위치에서 하나의 기울기를 가짐

  • 함수 𝑓가 유일한 접선을 가짐

• 𝑓의 그래프 위의 한 점 (𝑥₀, 𝑓(𝑥₀))에 접선이 존재

  • 접선은 접점 근처 함수의 곡선에 대한 좋은 근사

• 미분 조건 : 연속함수, 좌미분계수=우미분계수(기울기가 동일)

• 𝑓의 그래프 위의 한 점 (𝑥0, 𝑓(𝑥₀))에 접선이 존재

  • 접선 : \(𝑙(𝑥) = 𝑓′(𝑥_0)(𝑥 − 𝑥_0) + 𝑓(𝑥_0)\)

• 𝑓의 근에 대한 근사로 𝑙의 근을 선택

• l의 근 : \(x_1 = x_0 - \frac{f(x_0)}{f’(x_0)} \), \(x_2 = x_1 - \frac{f(x_1)}{f’(x_1)} $ , $ x_3 = x_2 - \frac{f(x_2)}{f’(x_2)} …\)

• 이러한 점들의 수열은 𝑓의 근으로 수렴

뉴턴법 - 초기값 설정(나쁜 초기값)

• 발산, 평평한 점, 순환 그래프는 나쁜 초기값을 초래한다.

• 적절한 초기값을 설정해야 뉴턴법이 정상 동작한다.

• 이분법을 적용하여 근이 존재할 범위 파악 후 적절한 초기값을 설정

import numpy as np
import sympy as sp
def Newton(f, fprime, a, nmax, eps, delta):
    fx = f.evalf(subs={x:a})
    for i in range(1, nmax+1):
        fp = fprime.evalf(subs={x:a})
        if(abs(fp)<delta):
            print("small derivative")
            return
        d = fx/fp
        a = a - d
        fx = f.evalf(subs={x:a})
        if(abs(d)<eps):
            print("convergence at "+str(i)+", x = "+str(a)+", f(x) = "+str(fx))
            return

x=sp.symbols("x")
f=21*x**3 - 150*x**2 + 41*x - 220

이분법과 뉴턴법 코드 비교

x = sp.symbols("x")
f = x**3 + 2*x**2 + 10*x + 23
fprime = 3*x**2 + 4*x + 10
print("Bisection") # 이분법
Bisection(f,-10,10,1E-7,100)
print("Newton") # 뉴턴법
Newton(f,fprime,0,100,1E-7,1E-7)
# 뉴턴은 5번, 이분법은 28번 반복했음. 뉴턴이 훨씬 빠름.
# f(x)는 뉴턴이 훨씬 작으니까 정밀도도 훨씬 높음.

Bisection
convergence at 28, x = -2.202034369111061, f(x) = 6.71329234040575e-7
Newton
convergence at 5, x = -2.20203441176566, f(x) = -1.37349444519054e-15

할선법

할선법 Method

• 뉴턴법만큼 빠르게 근으로 수렴하면서, 쉽게 계산하는 방법

  • 뉴턴법 : \(x_{n+1} = x_n - \frac{f(x_n)}{f’(x_n)}\)

  • \(𝑓′(𝑥_𝑛)\) 을 쉽게 계산할 수 있는 근삿값으로 교체

• 할선법 : \(x_{n+1} = x_n - (\frac{x_n-x_{n-1}}{f(x_n)-f(x_{n-1})}) f(x_n) \)

import numpy as np
import sympy as sp

def swap(a, b):
    temp = a
    a = b
    b = temp
    return a, b

def Secant(f, a, b, nmax, eps):
    fa = f.evalf(subs={x:a}) # f(a)
    fb = f.evalf(subs={x:b}) # f(b)
    if abs(fa) > abs(fb):
        a, b = swap(a, b)
        fa, fb = swap(fa, fb)
    for i in range(2, nmax+1):
        if abs(fa) > abs(fb):
            a, b = swap(a, b)
            fa, fb = swap(fa, fb)
        d = (b-a)/(fb-fa)
        b = a
        fb = fa
        d = d*fa # a-d 하기전에 eps와 비교위해(오차범위 확인)
        if abs(d) < eps:
            print("convergence at " +str(i) +", x = "+str(a)+", f(x) = "+str(fa))
            return
        a = a-d # 다음 x값을 구한것임.
        fa = f.evalf(subs={x:a}) # f(a)

이분법, 뉴턴법, 할선법 전체 비교

x = sp.symbols("x")
# 식 1 : f=21*x**3-150*x**2+41*x-220
# 식 2 : f=22*x**3-25*x**2+21*x-100
f1 = 21*x**3-150*x**2+41*x-220
fprime1 = 63*x**2-300*x+41
f2 = 22*x**3-25*x**2+21*x-100
fprime2 = 66*x**2-50*x+21
print("Bisection") # 이분법
Bisection(f1,-10,10,1E-7,100)
Bisection(f2,-10,10,1E-7,100)
print("Newton") # 뉴턴법
Newton(f1,fprime1,0,100,1E-7,1E-7)
Newton(f2,fprime2,0,100,1E-7,1E-7)
print("Secant") # 할선법
Secant(f1, -10, 10, 100, 1E-7)
Secant(f2, 0, 1, 100, 1E-7)

Bisection
convergence at 28, x = 7.07616962492466, f(x) = -1.79966109081769e-5
convergence at 28, x = 1.8966569751501083, f(x) = -4.03194505697596e-6
Newton
convergence at 8, x = 7.07616964170165, f(x) = -1.16819010758179e-13
convergence at 8, x = 1.89665699979684, f(x) = -8.21520068727071e-15
Secant
convergence at 7, x = 7.07616964172657, f(x) = 2.67377372820549e-8
convergence at 12, x = 1.89665702055653, f(x) = 3.39606631760622e-6