[numerical computation] Bisection Method (w/ MATLAB)

Bisection Method

수치해석에서 이분법(Bisection Method)은 근을 탐색하는 방법 중 하나입니다.

Bisection Method

 

Bisection Method는 Bracketing methods로 폐구간을 설정하여 근을 구하는 방법.

 

수학에서 이분법은 근이 반드시 존재하는 폐구간을 이분한 후, 이 중 근이 존재하는 하위 폐구간을 선택하는 것을 반복하여서 근을 찾는 알고리즘이다. 

 

간단하고 견고하며, 근의 대략적 위치를 안다면 일정 오차 내에 있는 1개의 해는 무조건 도출이 가능하나, 상대적으로 느린 방법이다. 구간안에 근이 존재한다는 전제하에 도출하는 방법이므로, 근이 존재하는 구간을 미리 파악하여야 한다. 

 

Algorithm

이분법을 통해 근을 계산하는 방법은 다음과 같습니다.

 

1. 탐색 구간 $[x_l, x_u]$과 임계 기준값(tolerance)를 선정합니다.
2. 선정된 구간의 양 끝점에서의 함수 값을 곱하여 음수인지 양수인지 확인합니다.
   $$ f(x_l)\  f(x_u) > 0 $$
3. 근이 존재한다면 추정되는 근 $x_r$을 계산합니다.
   $$ x_r = \frac{x_l+x_u}{2}$$
4. $f(x_u)$혹은 $f(x_l)$중에서 이전 단계에서 추정근 $f(x_r)$의 오차를 계산한다. 
   $$abs(\ f(x_r)-f(x_l)\ ) < tol$$
5. 설정해둔 임계 기준값 안에 오차값이 들어오면 $x_r$이 추정되는 근의 값이 된다.
   $$Approximate \ Solution = x_r$$
6. 만약, 오차값이 임계 기준값 안으로 들어오지 못하면, $x_r$혹은 $x_l$을 $x_r$로 업데이트하여 3번으로 되돌아가 다시 오차값안으로 들어올때까지 iterate한다. 

Result

Bisection Method를 이용한 근 찾기

Iter 횟수 및 결과

MATLAB 코드

clc;
clear all;

x_l = 1; x_u = 2; x_r = 0;
tol = 0.001;

y_l = func(x_l);
y_u = func(x_u);

if(func(x_l) * func(x_u) > 0)
	fprintf(1, "No root exists\n");
else
    iter = 0;
    fprintf(1, '\n\n');
    fprintf(1,' iter|    a     |     b    |     c    |    b-c   | f(b)*f(c)\n'); 
    fprintf(1,'------------------------------------------------------------\n');
    while(1)
        iter = iter + 1;
        x_r = (x_l + x_u) / 2;
        y_r = func(x_r);
        fprintf(1,'  %d  | %f | %f | %f | %f | %f\n',iter, x_l, x_u, x_r, abs(x_u-x_r),y_l*y_u);
        if(abs(x_r - x_l) < tol)
        	fprintf(1,'Approximate solution c = %.7f\n', x_r);
            break;
        elseif(y_l * y_r < 0)
             x_u = x_r;
             y_u = y_r;
        else
             x_l = x_r;
             y_l = y_r;
        end
    end
end

x = linspace(-1, 3);
y = func(x);

figure;
plot([-1,3], [0,0], x,y, '-', x_r, y_r, '*', LineWidth=2);
grid on;
% % set(gca, 'XAxisLocation', 'origin');
set(gca, 'YAxisLocation', 'origin');

function y = func(x)
	y = x.^3 - x.^2 - x - 1;
end