비선형방정식의 근

Posted on In , , Reading time ≈ 6 mins.

수치적인 근

방정식의 근은 오차범위 내에서만 파악 가능
ex) 10x=310x = 3의 정확한 근을 표현할 수 없음

따라서 방정식의 근을 찾는다는건 사실 ∣f(x)∣<ε|f(x)| < \varepsilon이 되는 값의 범위을 찾는다고 생각해야 함.
실수를 다룰 때는 정확하게 a==b인지 검사할 수 없다!!

이분법

f(x)f(x)가 연속함수고 f(a)f(b)<0f(a)f(b) < 0이라면 f(x)=0f(x)=0인 x∈(a,b)x \in (a,b)가 존재한다.

매번 c=a+b2c = \frac{a+b}{2}를 계산해서 f(a)f(c)f(a)f(c)나 f(c)f(b)f(c)f(b) 중 음수인 쪽으로 구간을 계속 절반으로 줄이면 된다.
Machine epsilon 때문에 log⁡2∣a−b∣ε\log_2\frac{|a-b|}{\varepsilon}번 이내로 반드시 해가 구해진다는 장점이 있음!
특히 함수가 비해석적이여도 연속적이면 무조건 적용가능!!

하지만 복소함수로 간다면 해를 찾을 수 없음..

뉴턴법

함수 f의 정확한 근이 x∗x^*이고 그 근처의 점 xx가 있을 때, 테일러 전개에 의해

0=f(x∗)=f(x)+(x∗−x)f′(x)+⋯0 = f(x^*) = f(x) + (x^* - x)f'(x) + \cdots

첫 항만 보면 x∗≈x−f(x)f′(x)x^* \approx x - \frac{f(x)}{f'(x)}라고 볼 수 있다!
물론 절단오차가 있지만, 이 방법을 계속해서 적용시키면 해가 구해진다!
대부분 함수는 뉴턴법으로 풀리고, 특히 복소함수도 풀린다.

단점) f′(x)f'(x)가 0에 가까워 값이 급격하게 변하거나 f(x)f′(x)\frac{f(x)}{f'(x)} 부호가 계속 변해 진동만 할 경우 (함수가 대칭이면 자주 발생함) 해를 찾는데 오래 걸리거나 못 구할 수도 있다...

할선법

f′(x)f'(x)를 구하기 어려운 경우 아래 근사를 사용할 수 있다.

f′(x)≈fn−fn−1xn−xn−1∴xn+1=xn−f(xn)f′(xn)=xn−fnxn−xn−1fn−fn−1\begin{align*} f'(x) &\approx \frac{f_n - f_{n-1}}{x_n - x_{n-1}} \\ \therefore x_{n+1} &= x_n - \frac{f(x_n)}{f'(x_n)} \\ &= x_n - f_n \frac{x_n - x_{n-1}}{f_n - f_{n-1}} \end{align*}

뉴턴법과 달리 처음에 두 점이 필요하며, 반드시 ∣f0∣>∣f1∣|f_0| > |f_1|이여야 한다.

다중근의 경우

P(x)=(x−α)nQ(x)P(x) = (x - \alpha)^n Q(x)

다중근을 갖는 다항식의 경우 조그만 변화로도 오차가 급속도로 커지게 된다.
이 경우에는 미분한 값으로 나눠주면 다중근의 영향이 소거된다.

P(x)=(x−α)nQ(x)P′(x)=(x−α)n−1(nQ(x)+(x−α)Q′(x))f(x)=P(x)P′(x)=(x−α)Q(x)nQ(x)+(x−α)Q′(x)\begin{align*} P(x) &= (x - \alpha)^n Q(x) \\ P'(x) &= (x - \alpha)^{n-1} (nQ(x) + (x - \alpha) Q'(x)) \\ f(x) &= \frac{P(x)}{P'(x)} = (x - \alpha) \frac{Q(x)}{nQ(x) + (x - \alpha) Q'(x)} \end{align*}

f(x)f(x)는 x=αx = \alpha에서 단일근을 가지므로, 다중근을 가질 것으로 예상되는 경우 f(x)=0f(x) = 0의 해를 구하면 된다.
사실 QR 반복법으로 해를 구한 다음 다중근의 존재여부를 판단하는 방식이 상위호환이다

이 경우 뉴턴법을 적용시키기엔 f′(x)f'(x)가 너무 더러우므로 할선법을 사용한다.

연속대입법

g(x)=x,∣g′(x)∣<1g(x) = x, |g'(x)| < 1

위 형태로 주어진 방정식의 경우, xn+1=g(xn)x_{n+1} = g(x_n)을 반복하여 근을 구할 수 있다.

사실 뉴턴법과 할선법은 연속대입법의 일종으로 볼 수 있다!

xn+1=xn−f(xn)f′(xn)g(x)=x−f(x)f′(x)\begin{align*} x_{n+1} &= x_n - \frac{f(x_n)}{f'(x_n)} \\ g(x) &= x - \frac{f(x)}{f'(x)} \end{align*}

xn+1=xn−fnxn−xn−1fn−fn−1g(x)=x−f(x)x−x∗f(x)−f(x∗)\begin{align*} x_{n+1} &= x_n - f_n \frac{x_n - x_{n-1}}{f_n - f_{n-1}} \\ g(x) &= x - f(x) \frac{x - x^*}{f(x) - f(x^*)} \end{align*}

오차 및 수렴특성

이론상 반복하면 수렴하지만... 얼마나 반복해야할까?

실제 근을 α\alpha, 오차를 en=α−xne_n = \alpha - x_n이라고 하자.

연속대입법의 경우

α=g(α)\alpha = g(\alpha)이고 xk+1=g(xk)x_{k+1} = g(x_k)이므로, 평균값 정리로부터

∃ξk∈(α,xk) s.t.ek+1=α−xk+1=g(α)−g(xk)=(α−xk)g′(ξk)∴∣ek+1ek∣=∣g′(ξk)∣\begin{align*} \exists \xi_k &\in (\alpha, x_k) \ s.t. \\ e_{k+1} &= \alpha - x_{k+1} \\ &= g(\alpha) - g(x_k) \\ &= (\alpha - x_k)g'(\xi_k) \\ \therefore \left| \frac{e_{k+1}}{e_k} \right| &= |g'(\xi_k)| \end{align*}

따라서 항상 ∣g′(x)∣≤r<1|g'(x)| \leq r < 1이라면 ∣ek+1∣≤rk∣e1∣|e_{k+1}| \leq r^k |e_1|이므로 오차가 0으로 수렴한다.
이때 오차의 비율이 선형적으로 감소하므로, 연속대입법은 선형적 수렴정도를 가진다.

뉴턴법의 경우

f(α)=0f(\alpha) = 0이고 xk+1=xk−f(xk)f′(xk)x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}이므로

ek+1=α−xk+1=α−xk+f(xk)f′(xk)=ek+f(xk)−f(α)f′(xk)\begin{align*} e_{k+1} &= \alpha - x_{k+1} \\ &= \alpha - x_k + \frac{f(x_k)}{f'(x_k)} \\ &= e_k + \frac{f(x_k) - f(\alpha)}{f'(x_k)} \end{align*}

이 때 테일러 전개에 의해

∃ξk∈(α,xk) s.t.f(α)=f(xk)+(α−xk)f′(xk)+12(α−xk)2f′′(ξk)∴f(xk)−f(α)=−ekf′(xk)−12ek2f′′(ξk)∴ek+1=ek+−ekf′(xk)−12ek2f′′(ξk)f′(xk)=−12ek2f′′(ξk)f′(xk)\begin{align*} \exists \xi_k &\in (\alpha, x_k) \ s.t. \\ f(\alpha) &= f(x_k) + (\alpha - x_k)f'(x_k) + \frac{1}{2}(\alpha - x_k)^2f''(\xi_k) \\ \therefore f(x_k) - f(\alpha) &= -e_kf'(x_k) - \frac{1}{2}e_k^2f''(\xi_k) \\ \therefore e_{k+1} &= e_k + \frac{-e_kf'(x_k) - \frac{1}{2}e_k^2f''(\xi_k)}{f'(x_k)} \\ &= -\frac{1}{2}e_k^2\frac{f''(\xi_k)}{f'(x_k)} \end{align*}

이때 오차의 비율이 제곱에 비례하여 감소하므로, 뉴턴법은 2차적 수렴정도를 가진다.
한편 테일러 전개에서 ek2e_k^2 항을 무시하면 ek+1e_{k+1}이 제곱에 비례한다면서 이게 뭔 개소리지 무시할 수 있는거 맞나 ek=−f(xk)f′(xk)e_k = -\frac{f(x_k)}{f'(x_k)}이므로,

∣ek+1ek∣=∣−12ekf′′(ξk)f′(xk)∣=12∣f(xk)f′′(ξk)(f′(xk))2∣\left| \frac{e_{k+1}}{e_k} \right| = \left| -\frac{1}{2}e_k\frac{f''(\xi_k)}{f'(x_k)} \right| = \frac{1}{2} \left| \frac{f(x_k)f''(\xi_k)}{(f'(x_k))^2} \right|

이다. 따라서 12∣f(xk)f′′(ξk)(f′(xk))2∣<1\frac{1}{2} \left| \frac{f(x_k)f''(\xi_k)}{(f'(x_k))^2} \right| < 1이여야 뉴턴법이 수렴하게 된다.
일반적으로 ∣f′(xk)∣|f'(x_k)|가 큰 값을 가지거나 변곡점이 없을 때 뉴턴법이 잘 수렴한다.

또한, ek=−f(xk)f′(xk)e_k = -\frac{f(x_k)}{f'(x_k)}, xk+1=xk−f(xk)f′(xk)x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}에서 ek≈xk+1−xke_k \approx x_{k+1} - x_k다. 오차를 계산 도중에 예측할 수 있다!

가위치법

안 쓰니까 재미로 보기

이분법은 함수값을 사용하지 않는다 -> 함수가 직선에 가깝다면 직선을 이용하는게 더 효과적이지 않을까?
단순 cn=an+bn2c_n = \frac{a_n + b_n}{2} 대신 (an,f(an)),(bn,f(bn))(a_n, f(a_n)), (b_n,f(b_n))을 지나는 직선이 x축과 만나는 점 cn=an−f(an)an−bnf(an)−f(bn)c_n = a_n - f(a_n) \frac{a_n - b_n}{f(a_n) - f(b_n)}을 새로운 경계점으로 사용하는 방법을 가위치법이라고 한다.

문제점) 함수가 직선에 가깝지 않을수가 있음...
함수값의 변화가 완만한 경우 한쪽의 경계점이 계속 고정되는 단점이 있다 (정체점 문제)

수정 가위치기법

이전과 똑같은 경계점은 (an,f(an))(a_n, f(a_n)) 대신 (an,f(an)2)(a_n, \frac{f(a_n)}{2})를 사용한다!
의외로 이렇게 하면 뉴턴법 급으로 빠르게 수렴한다고 한다...

뮬러법

안 쓰니까 재미로 보기 2
진짜 몰라도 됨 시험도 안 냄

왜 두 개만 쓰지? 점 세 개를 가지고 포물선을 만들자!

  1. x1<x2<x3x_1 < x_2 < x_3 세 점을 잡는다. (놀랍게도 함수값 조건은 없다! 대신 함수값 부호가 전부 같으면 [x1,x3][x_1, x_3] 밖에 근이 있을 수도 있으므로 수렴하는데 오래 걸릴 수 있다...)
  2. 세 점 (x1,f1),(x2,f2),(x3,f3)(x_1, f_1), (x_2, f_2), (x_3, f_3)을 지나는 이차방정식을 구하고, x2x_2와 제일 가까운 근 x0x_0을 구한다.
  3. x0,x1,x2,x3x_0, x_1, x_2, x_3 중 x0x_0과 제일 멀리 떨어져있는 점 하나를 버린다.
  4. 남은 세 점을 다시 x1<x2<x3x_1 < x_2 < x_3이 되도록 재배치 한 다음 위 과정을 반복한다.

당연히 이차방정식을 구하고 근을 구하는 과정이 제일 거지같은데, 세 점 (x1,f1),(x2,f2),(x3,f3)(x_1, f_1), (x_2, f_2), (x_3, f_3)을 지나는 이차방정식은

p(x)=f1+f2−f1x2−x1(x−x1)+f3−f2x3−x2−f2−f1x2−x1x3−x1(x−x1)(x−x2)p(x) = f_1 + \frac{f_2 - f_1}{x_2 - x_1} (x - x_1) + \frac{\frac{f_3 - f_2}{x_3 - x_2} - \frac{f_2 - f_1}{x_2 - x_1}}{x_3 - x_1}(x - x_1)(x - x_2)

이고, m1=f2−f1x2−x1,m2=f3−f2x3−x2,a=m2−m1x3−x1,b=m1+a(x2−x1),c=f2m_1 = \frac{f_2 - f_1}{x_2 - x_1}, m_2 = \frac{f_3 - f_2}{x_3 - x_2}, a = \frac{m_2 - m_1}{x_3 - x_1}, b = m_1 + a(x_2 - x_1), c = f_2로 정의하면

p(x)=f2+m1(x−x2)+a(x−x1)(x−x2)=f2+m1(x−x2)+a((x−x2)+(x2−x1))(x−x2)=a(x−x2)2+(m1+a(x2−x1))(x−x2)+f2=a(x−x2)2+b(x−x2)+c\begin{align*} p(x) &= f_2 + m_1(x - x_2) + a(x - x_1)(x - x_2) \\ &= f_2 + m_1(x - x_2) + a((x - x_2) + (x_2 - x_1))(x - x_2) \\ &= a(x - x_2)^2 + (m_1 + a(x_2 - x_1))(x - x_2) + f_2 \\ &= a(x - x_2)^2 + b(x - x_2) + c \end{align*}

가 되므로 위 이차방정식의 근은 x−x2=−b±b2−4ac2ax - x_2 = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a}다.
구간의 크기를 줄여 수렴시키기 위해서 x2x_2에 제일 가까운 근을 선택해야 하므로 분자가 제일 작은 근을 골라야 하는데, 너무 가까운 두 수를 빼면 오차가 커지는 문제점이 있다.
따라서 근의 공식을 변형시켜 루트 부분이 분모로 가게 만들고 분모가 제일 큰 경우를 선택한다.

∴x0={x2−2cb+b2−4acb>0x2−2cb−b2−4acb≤0\therefore x_0 = \begin{cases} x_2 - \frac{2c}{b + \sqrt{b^2 - 4ac}} & b > 0 \\ x_2 - \frac{2c}{b - \sqrt{b^2 - 4ac}} & b \leq 0 \end{cases}

베어스토우법

이것도 QR 반복법에 대체돼서 안 씀

n차 다항식이 사실 어떤 이차다항식으로 나눠진다고 가정

pn(x)=anxn+an−1xn−1+⋯+a1x+a0=(x2+rx+s)(bnxn−2+bn−1xn−3+⋯+b1x+b0)+b1x+b0\begin{align*} p_n(x) =& a_nx^n + a_{n-1}x^{n-1} + \cdots + a_1x + a_0 \\ =& (x^2 + rx + s)(b_nx^{n-2} + b_{n-1}x^{n-3} + \cdots + b_1x + b_0) \\ &+ b_1x + b_0 \end{align*}

만약 b1x+b0=0b_1x + b_0 = 0이라면, x2+rx+s=0x^2+rx+s = 0에서 두 개의 근을 구할 수 있음
r,sr, s를 b1x+b0=0b_1x + b_0 = 0이 되는 참값, r∗,s∗r^*, s^*를 추정값, 참값과 추정값의 오차를 Δr=r−r∗,Δs=s−s∗\Delta r = r - r^*, \Delta s = s - s^*라고 하자.
b0,b1b_0, b_1을 r,sr, s에 대한 함수로 보고 테일러 전개를 하면

0=b0(r,s)≈b0(r∗,s∗)+Δr∂b0∂r+Δs∂b0∂s0=b1(r,s)≈b1(r∗,s∗)+Δr∂b1∂r+Δs∂b1∂s∴ Δr∂b0∂r+Δs∂b0∂s=−b0(r∗,s∗)Δr∂b1∂r+Δs∂b1∂s=−b1(r∗,s∗)∴ [ΔrΔs]=[∂b0∂r∂b0∂s∂b1∂r∂b1∂s]−1[−b0(r∗,s∗)−b1(r∗,s∗)]\begin{align*} 0 = b_0(r,s) &\approx b_0(r^*, s^*) + \Delta r \frac{\partial b_0}{\partial r} + \Delta s \frac{\partial b_0}{\partial s} \\ 0 = b_1(r,s) &\approx b_1(r^*, s^*) + \Delta r \frac{\partial b_1}{\partial r} + \Delta s \frac{\partial b_1}{\partial s} \\ \therefore\ & \Delta r \frac{\partial b_0}{\partial r} + \Delta s \frac{\partial b_0}{\partial s} = -b_0(r^*, s^*) \\ &\Delta r \frac{\partial b_1}{\partial r} + \Delta s \frac{\partial b_1}{\partial s} = -b_1(r^*, s^*) \\ \therefore\ & \begin{bmatrix} \Delta r \\ \Delta s \end{bmatrix} = \begin{bmatrix} \frac{\partial b_0}{\partial r} & \frac{\partial b_0}{\partial s} \\ \frac{\partial b_1}{\partial r} & \frac{\partial b_1}{\partial s} \end{bmatrix}^{-1} \begin{bmatrix} -b_0(r^*, s^*) \\ -b_1(r^*, s^*) \end{bmatrix} \end{align*}

(bi)r≔∂bi∂r(b_i)_r \coloneqq \frac{\partial b_i}{\partial r}, (bi)s≔∂bi∂s(b_i)_s \coloneqq \frac{\partial b_i}{\partial s}라고 정의하면 (b0)r,(b0)s,(b1)r,(b1)s(b_0)_r, (b_0)_s, (b_1)_r, (b_1)_s만 구하면 Δr,Δs\Delta r, \Delta s를 구할 수 있다!!!
왜 아무도 안 쓰는지 알겠다
미분값을 구하는것도 정말 아름다운데, 위의 방정식을 전개한 뒤 계수비교하면

bn=anbn−1=an−1−rbnbn−2=an−2−rbn−1−sbn⋮b1=a1−rb2−sb3b0=a0−sb2\begin{align*} b_n &= a_n \\ b_{n-1} &= a_{n-1} - rb_n \\ b_{n-2} &= a_{n-2} - rb_{n-1} - sb_n \\ &\vdots \\ b_1 &= a_1 - rb_2 - sb_3 \\ b_0 &= a_0 - sb_2 \end{align*}

이므로, 간단한 편미분을 통해

(bn)r=0(bn−1)r=−bn(bn−2)r=−bn−1−r(bn−1)r−s(bn)r⋮(b1)r=−b2−r(b2)r−s(b3)r(b0)r=−s(b2)r(bn)s=0(bn−1)s=0(bn−2)s=−bn−s(bn)s−r(bn−1)s⋮(b1)s=−b3−s(b3)s−r(b2)s(b0)s=−b2−s(b2)s\begin{aligned} (b_n)_r &= 0 \\ (b_{n-1})_r &= -b_n \\ (b_{n-2})_r &= -b_{n-1} - r(b_{n-1})_r - s(b_n)_r \\ &\vdots \\ (b_1)_r &= -b_2 - r(b_2)_r - s(b_3)_r \\ (b_0)_r &= -s(b_2)_r \end{aligned} \quad \begin{aligned} (b_n)_s &= 0 \\ (b_{n-1})_s &= 0 \\ (b_{n-2})_s &= -b_n -s(b_n)_s - r(b_{n-1})_s \\ &\vdots \\ (b_1)_s &= -b_3 - s(b_3)_s - r(b_2)_s \\ (b_0)_s &= -b_2 - s(b_2)_s \end{aligned}

임을 알 수 있다.
따라서 ai,bi,(bi)r,(bi)sa_i, b_i, (b_i)_r, (b_i)_s 배열 4개를 준비한 다음 간단한 DP를 통해 Δr,Δs\Delta r, \Delta s를 구할 수 있다...
Δr,Δs\Delta r, \Delta s가 충분히 작아질 때까지 반복하면 된다...

근데 이따구로 구해도 다중근의 경우 성능이 떨어지는 단점이 있다...

0%