4. 방정식

비선형 방정식

근을 찾아라

구간 $[a, b]$가 있다.
$f(x)$$[a, b]$에서 연속이고 $f(a) \cdot f(b) < 0 \implies$ 구간에 근이 1개 이상 존재

그렇다면 $[a, b]$를 어떻게 알아내는가?

  • 딱히 일반화된 방법은 없다
  • $f(a)$에서 시작해서 일정 간격 숫자를 변화시켜가며 함수값의 부호가 바뀔 때까지 계속
  • 주의:
    • 근이 복수일 때, 두 근의 차가 간격보다 작아서 캐치되지 않을 수 있다

이분양단법(바이섹션)

  • 근을 찾는 가장 간단한 방법
  • 가정:
    • $f(x)$가 범위 $[ a, b ]$에서 근을 가지고 연속이다.
    • $f(x)$가 범위 $[ a, b ]$에서 부호가 바뀐다. $f(a) \cdot f(b) < 0$
    • 정확한 근을 찾을 수는 없지만, 근사값 $\bar{\alpha} \quad \mathrm{s.t.} \quad \left\lvert \alpha - \bar{\alpha} \right\rvert \le \epsilon$ 을 찾는다. $\epsilon > 0$은 오차범위(error tolerance)
  1. $c = (a + b) /2$를 정의한다.
  2. $b - c \le \epsilon \implies c = \alpha,$ 계산 끝
  3. $b - c > \epsilon$이면 부호에 따라 결정한다.
    • $\begin{cases} f(b) \cdot f(c) \le 0 & \implies c \longrightarrow a \\ f(b) \cdot f(c) > 0 & \implies b \longrightarrow c \end{cases}$
    • 1단계로 돌아가 반복한다.

코드

#include <stdio.h>
#include <math.h>
 
float func(float x)        // define equation to solve
{
    float y=pow(x,6)-x-1;
    return y;
}
 
int main(void)
{
    float a=1;
    float b=2;
    float c=(a+b)/2;
    float err=10e-3;
 
    do
      {
        if(func(b) * func(c) <= 0)
        a=c;
        else b=c;
        c=(a+b)/2;
        printf("root = %f \n", c);
      }
 
    while ( b-c > err );
 
    return 0;
}

장점

  • 언제나 수렴한다.
  • 오차범위가 주어져 있고, 루프가 반복될 때마다 그 범위는 줄어든다.
  • 수렴률이 정해져 있다. 루프마다 오차범위는 ½씩 줄어든다.

단점

  • 느리다.
  • 엡실론 값이 컴퓨터의 한계보다 작지 않은지 확인해야 한다 (이것은 계산기 엡실론을 확인해 보면 쉬이 알 수 있다)

이분양단법으로 근처문제(Nearby Problem) 풀기

nearby.png
(1)
\begin{equation} p(x) = f(x_0) + (x - x_0 ) f' (x_0) \end{equation}
  • $x_0$ 가 근의 참값 $\alpha$ 에 매우 가깝다면, $p(x)$의 근도 $\alpha$ 에 가까워야 한다.
  • 근의 근사값 $x_1 = x_0 - {{f(x_0)} \over {f'(x_0)}}$
  • 이하 반복.

뉴턴-랩슨법(NR)

  • 방정식 $f(x) = 0$에 대하여 근 $\alpha$의 최초 예측값 $x_0$가 주어진다.
(2)
\begin{align} x_{n+1} = x_n - {{f(x_n)} \over {f'(x_n)}}, \qquad n = 0, 1, 2, \cdots \end{align}

예: $f(x) \equiv x^6 - x - 1 = 0$의 근사

$n$ $x_n$ $f(x_n)$ $x_n - x_{n-1}$
0 1.5 8.89e+1
1 1.30049088 2.54e+1 -2.00e-1
2 1.18148042 5.38e-1 -1.19e-1
3 1.13945559 4.92e-2 -4.20e-2
4 1.13477763 5.50e-4 -4.68e-3
5 1.13472415 7.11e-8 -5.35e-5
6 1.13472414 1.55e-15 -6.91e-9

코드:

#include <stdio.h>
#include <math.h>
 
float hamsu(float x)
{
    float y=pow(x,6)-x-1;
    return y;
}
 
float dohamsu(float x)
{
    float y=6*pow(x,5)-1;
    return y;
}
 
int main(void)
{
    float x0 = 1.5;
    float x1 = x0 -(hamsu(x0)/dohamsu(x0));
    float ocha=10e-6;
 
    do
    {
        x0 = x1;
        x1 = x0 -(hamsu(x0)/dohamsu(x0));
        printf("%f \n", x1);
    }
    while(sqrt(pow(x1-x0,2)) > sqrt(pow(ocha,2)));
 
    return 0;
}

오차공식
$\alpha, x_n$ 사이의 어떤 $c_n$에 대하여

(3)
\begin{align} f(\alpha) = f(x_n) + (\alpha - x_n ) f' (x_n) + {1 \over 2} ( \alpha - x_n )^2 f'' (c_n) = 0 \end{align}

양변을 $f'(x_n)$으로 나누면

(4)
\begin{align} 0 = {{f(x_n)} \over {f'(x_n)}} + \alpha - x_n + (\alpha - x_n )^2 {{f''(c_n)} \over {2 f'(x_n)}} \end{align}

$f(x_n) / f'(x_n) - x_n = - x_{n+1}$이므로

(5)
\begin{align} \alpha - x_{n+1} = -{{f''(c_n) } \over {2 f'(x_n)}} ( \alpha - x_n )^2 \end{align}

$x_n$$\alpha$에 매우 가까우므로 $c_n$$\alpha$에 매우 가깝다

(6)
\begin{align} \alpha - x_{n+1} \approx - {{f'' (x_n)} \over {2 f'(\alpha)}} ( \alpha - x_n )^2 \end{align}

NR법은 이차수렴하며, $f'(\alpha ) \ne 0$이고 $f(x)$$\alpha$ 근처에서 이계미분가능함수여야 한다.

NR의 실패

fail.png
  • NR은 근의 참값 $\alpha$ 근처에서만 작동한다. 기울기가 0으로 가면 NR은 망한다.
  • 이계도함수의 부호가 근에서 바뀔 경우, NR은 절대 수렴하지 않는다.

장점

  • 빠르다.
  • 공식이 단순해서 적용과 코딩이 쉽다.
  • 직관적이라 그 행동을 이해하기 쉽다.

단점

  • 수렴하지 않을 수 있다.
  • $f'(\alpha) = 0$이라면 골치아파진다.
  • $f(x)$$f'(x)$를 모두 알아야 한다. 이분양단법은 $f(x)$만 알아도 된다.

할선법

원리

  • 근의 추정값을 두 개 선택해서 $x_0, x_1$라 한다.
  • $(x_0, f(x_0)), (x_1, f(x_1))$를 잇는 직선을 그으면 이 직선은 $f(x)$를 관통하는 할선(secant시컨트)이 된다.
(7)
\begin{align} q(x) = {{(x_1 - x) f(x_0) + (x-x_0) f(x_1)} \over {x_1 - x_0}} \end{align}
  • $q(x) = 0$인 점을 찾아 거기를 $(x_2, 0)$으로 삼는다.
(8)
\begin{align} x_2 = x_1 - f(x_1) \div {{f(x_1) - f(x_0)} \over {x_1 - x_0}} \end{align}
  • 이하 반복한다.
(9)
\begin{align} x_{n+1} = x_n - f(x_n) \div {{f(x_n) - f(x_{n-1}) } \over {x_n - x_{n-1} }} \\ n = 1, 2, 3, \cdots \end{align}

할선법의 알고리듬

  1. x1, x2를 정의한다.
  2. f1 = F(x1), f2 = F(x2)를 계산한다.
  3. x3 = x2 - f2 * (x2 - x1)/(f2 - f1) 를 계산한다.
  4. | x3 - x1 | < TOL 이면 멈춘다. x3이 해이다.
  5. | x3 - x1 | > TOL 이면 x1 = x2, x2 = x3 으로 정하고 2단계를 반복한다.

예: $f(x) \equiv x^6 -x -1 = 0$ 풀기

코드:

#include <stdio.h>
#include <math.h>
 
float func(float x)
{
    float y=pow(x,6) -x -1;
    return y;
}
 
int main(void)
{
    double x1=2;
    double x2=1;
    double x3;
    double tol=1e-6;
 
    x3 = x2 - func(x2)*(x2 - x1)/(func(x2) - func(x1));
 
    do
    {
        printf("%e \n", x3);
 
        if (func(x1)*func(x2) < 0) x2 = x3;
        else    x1 = x2;
            x2 = x3;
 
        x3 = x2 - func(x2)*(x2 - x1)/(func(x2) - func(x1));
 
    } while(abs(x3 - x1) > tol);
 
    return 0;
 
}

고정소수점 반복법